From bfbb572ba47c2d4a7be553c125454dfe5f99dde8 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Mon, 3 Aug 2026 16:06:01 +0200 Subject: [PATCH 1/2] Fix double-counted prey-bin quadrature in getDiet() and getTrophicLevel() MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Under `second_order_w(params) <- TRUE`, `setPredKernel()` builds `ft_pred_kernel_e` as the predation kernel integrated over the prey bin, dividing by `beta - 1` so that the plain `w * dw` supplied by the prey vector cancels. `mizerEncounter()` accordingly keeps the point weight. `getDiet()` did not: it weighted its prey vectors with the bin-averaged `w`, applying the prey-bin quadrature a second time. The result was uniformly too large by `(1 + beta) / 2` — 9.7% for `NS_params`. Both branches (FFT and custom kernel) were affected. The default `proportion = TRUE` hid it, because a uniform factor divides out in the normalisation. `getTrophicLevel()` had the same mismatch plus a second one: its trophic-level-weighted numerator was built from the point-sampled `getPredKernel()` while its denominator came from `getEncounter()`, which uses the bin-integrated kernel once the flag is on. The two are then no longer the same integral and the trophic levels moved by up to 0.06. Both now use the plain point weight `w * dw`. New internal `encounter_kernel()` returns the kernel `mizerEncounter()` actually uses — recovered from `ft_pred_kernel_e` by inverse FFT, so it cannot drift from whatever quadrature `setPredKernel()` chose — and `getTrophicLevel()` builds its numerator from that. `getPredKernel()` is refactored onto the shared `expand_kernel_offsets()` helper; its behaviour is unchanged. Tests cover the `encounter_kernel()`/`getEncounter()` identity, the diet-sums-to-consumption identity, and a predator whose prey all have trophic level 1 coming out at exactly 2 — each across both settings of the flag and both kernel branches. All fail against the previous code. Fixes #474 Co-Authored-By: Claude Opus 5 --- NEWS.md | 21 ++++++ R/setPredKernel.R | 76 ++++++++++++++++++-- R/summary_methods.R | 65 +++++++++--------- man/encounter_kernel.Rd | 40 +++++++++++ man/expand_kernel_offsets.Rd | 29 ++++++++ man/getDiet.Rd | 5 +- tests/testthat/test-second_order_summary.R | 80 ++++++++++++++++++++++ 7 files changed, 280 insertions(+), 36 deletions(-) create mode 100644 man/encounter_kernel.Rd create mode 100644 man/expand_kernel_offsets.Rd diff --git a/NEWS.md b/NEWS.md index 7e773452f..6b074d2d4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -64,6 +64,27 @@ stability of steady states. `extinction_threshold` fraction (default `1e-6`) of its value at the start of the run. +## Bug fixes + +- `getDiet(proportion = FALSE)` no longer overcounts when second-order + bin-averaging is switched on with `second_order_w()`. It was applying the + prey-bin quadrature twice — once through its bin-averaged prey weight and + again through the bin-integrated predation kernel that `setPredKernel()` + builds under `second_order_w` — so the diet was uniformly too large by + `(1 + beta) / 2`, where `beta` is the grid ratio (9.7% for `NS_params`). + Summing the diet over prey now reproduces + `getEncounter() * (1 - getFeedingLevel())` under both schemes, on both the + FFT and the custom-kernel path. `getDiet(proportion = TRUE)`, the default, + was unaffected because the factor was uniform and divided out (#474). + +- `getTrophicLevel()` had the same quadrature mismatch: its + trophic-level-weighted numerator was built from the point-sampled + `getPredKernel()` and a bin-averaged prey weight, while its denominator came + from `getEncounter()`. Under `second_order_w` the two are no longer the same + integral, so the reported trophic levels were off by up to 0.06. Numerator + and denominator now use the same quadrature, and a predator whose prey all + have trophic level 1 comes out at exactly 2 in both schemes (#474). + # mizer 3.2.1 This patch release fixes how species and gear parameters are handled when they diff --git a/R/setPredKernel.R b/R/setPredKernel.R index baf6645d3..0cd713400 100644 --- a/R/setPredKernel.R +++ b/R/setPredKernel.R @@ -341,20 +341,88 @@ getPredKernel.MizerParams <- function(params) { phis <- get_phi(species_params, ppmr) # Do not allow feeding at own size phis[, 1] <- 0 + expand_kernel_offsets(phis, params, species_params$species) +} + +#' Expand kernel weights indexed by grid offset into a full kernel array +#' +#' The predation kernel depends only on the predator/prey mass ratio, so on the +#' geometric grid it is a function of the offset \eqn{m} between the predator +#' and prey grid indices alone. `phis[i, m + 1]` holds the weight of species +#' \eqn{i} at offset \eqn{m}. This helper writes those weights into the +#' (predator species x predator size x prey size) array that the non-FFT code +#' paths work with. +#' +#' @param phis A species-by-offset matrix of kernel weights, with the offset +#' running from 0 to `length(params@w_full) - 1`. +#' @param params A MizerParams object supplying the grid. +#' @param species A character vector of species names for the dimnames. +#' @return An array (predator species x predator size x prey size). +#' @concept helper +#' @keywords internal +expand_kernel_offsets <- function(phis, params, species) { + no_sp <- length(species) + no_w <- length(params@w) + no_w_full <- length(params@w_full) pred_kernel <- array(0, dim = c(no_sp, no_w, no_w_full), - dimnames = list(sp = species_params$species, + dimnames = list(sp = species, w_pred = signif(params@w, 3), w_prey = signif(params@w_full, 3))) - for (i in 1:no_sp) { - min_w_idx <- no_w_full - no_w + 1 + min_w_idx <- no_w_full - no_w + 1 + for (i in seq_len(no_sp)) { for (k in seq_len(no_w)) { pred_kernel[i, k, (min_w_idx - 1 + k):1] <- phis[i, 1:(min_w_idx - 1 + k)] } } - return(pred_kernel) + pred_kernel +} + +#' The predation kernel as used by the encounter quadrature +#' +#' Returns the kernel array \eqn{\Phi_i(w_k, w_p)} for which +#' \deqn{E_i(w_k) = \gamma_i(w_k) \sum_p \Phi_i(w_k, w_p) N^{eff}_i(w_p) +#' w_p \Delta w_p} +#' reproduces exactly the available energy computed by [mizerEncounter()], +#' where \eqn{N^{eff}} is the interaction-weighted prey density. It is the +#' kernel that any summary function must use if its result is to be consistent +#' with [getEncounter()]. +#' +#' On the default first-order path this is just the point-sampled kernel +#' returned by [getPredKernel()]. When second-order bin-averaging is switched on +#' (see [second_order_w()]) the two differ: `setPredKernel()` then builds the +#' Fourier-transformed kernel from the kernel *integrated over the prey bin*, +#' divided by \eqn{\beta - 1} so that the plain point weight \eqn{w_p \Delta +#' w_p} carried by the prey vector is cancelled. Those bin-integrated weights +#' are recovered here from `params@ft_pred_kernel_e` by an inverse Fourier +#' transform, which costs one FFT and keeps this helper automatically in step +#' with whatever quadrature `setPredKernel()` used. +#' +#' A summary function that instead pairs the point-sampled [getPredKernel()] +#' with a bin-averaged prey weight double-counts the prey-bin quadrature; that +#' was the bug behind issue #474. +#' +#' @param params A MizerParams object. +#' @return An array (predator species x predator size x prey size). +#' @concept helper +#' @keywords internal +encounter_kernel <- function(params) { + # A kernel that is stored explicitly is used as-is by mizerEncounter(), + # which weights it with the plain `w * dw`, so it is already consistent. + if (length(dim(params@pred_kernel)) > 1) { + return(params@pred_kernel) + } + if (!isTRUE(params@second_order_w[["bin_average"]])) { + return(getPredKernel(params)) + } + no_w_full <- length(params@w_full) + # setPredKernel() stores ft_pred_kernel_e[i, ] = fft(phi_e[i, ]), so the + # real-space weights come back from the inverse transform. + phis <- Re(base::t(mvfft(base::t(params@ft_pred_kernel_e), + inverse = TRUE))) / no_w_full + expand_kernel_offsets(phis, params, params@species_params$species) } #' @rdname setPredKernel diff --git a/R/summary_methods.R b/R/summary_methods.R index b8f256d10..9d0ddbc48 100644 --- a/R/summary_methods.R +++ b/R/summary_methods.R @@ -61,7 +61,10 @@ NULL #' This function performs the same integration as [getEncounter()] but does not #' aggregate over prey species, and multiplies by \eqn{1-f_i(w)} to get the #' consumed biomass rather than the available biomass. Outside the range of -#' sizes for a predator species the returned rate is zero. +#' sizes for a predator species the returned rate is zero. Summing the result +#' of `getDiet(proportion = FALSE)` over prey therefore reproduces +#' `getEncounter(params) * (1 - getFeedingLevel(params))`, whichever quadrature +#' scheme the model uses (see [second_order_w()]). #' #' @param object A \linkS4class{MizerParams} or \linkS4class{MizerSim} object. #' @param proportion If TRUE (default) the function returns the diet as a @@ -139,17 +142,16 @@ getDiet.MizerParams <- function(object, # object@w_full[idx_sp] = object@w idx_sp <- (no_w_full - no_w + 1):no_w_full - # Prey-biomass weight factor K = w. On the default path `w_eff` and - # `w_full_eff` are just the grid weights `w` and `w_full`, so the - # quadratures below are byte-identical to previous mizer versions. When - # second-order is enabled they become the trapezoidal bin-averages of `w`. - if (isTRUE(params@second_order_w[["bin_average"]])) { - w_eff <- bin_average_weight(params@w) - w_full_eff <- bin_average_weight(params@w_full) - } else { - w_eff <- params@w - w_full_eff <- params@w_full - } + # The prey vectors below carry the plain point weight `w * dw`, exactly as + # in mizerEncounter(). This is deliberate also when second-order + # bin-averaging is switched on: there the prey-bin integral has already been + # folded into the kernel by setPredKernel(), which builds + # `ft_pred_kernel_e` as the kernel integrated over the prey bin and divides + # out the `w * dw` that the prey vector supplies. Weighting the prey vector + # by the bin-averaged `w` as well would apply that quadrature twice and + # inflate the diet by (1 + beta) / 2 (issue #474). Summed over prey, the + # diet must reproduce `getEncounter() * (1 - getFeedingLevel())` under both + # schemes. # If the user has set a custom kernel we can not use fft. if (!is.null(comment(params@pred_kernel))) { @@ -159,16 +161,16 @@ getDiet.MizerParams <- function(object, # multiplication for this. Then we multiply 1st and 3rd ae <- matrix(params@pred_kernel[, , idx_sp, drop = FALSE], ncol = no_w) %*% - t(sweep(n, 2, w_eff * params@dw, "*")) + t(sweep(n, 2, params@w * params@dw, "*")) diet[, , 1:no_sp] <- ae # Eating the resource diet[, , no_sp + 1] <- rowSums(sweep( - params@pred_kernel, 3, params@dw_full * w_full_eff * n_pp, "*"), + params@pred_kernel, 3, params@dw_full * params@w_full * n_pp, "*"), dims = 2) } else { prey <- matrix(0, nrow = no_sp + 1, ncol = no_w_full) - prey[1:no_sp, idx_sp] <- sweep(n, 2, w_eff * params@dw, "*") - prey[no_sp + 1, ] <- n_pp * w_full_eff * params@dw_full + prey[1:no_sp, idx_sp] <- sweep(n, 2, params@w * params@dw, "*") + prey[no_sp + 1, ] <- n_pp * params@w_full * params@dw_full ft <- array(rep(params@ft_pred_kernel_e, times = no_sp + 1) * rep(mvfft(t(prey)), each = no_sp), dim = c(no_sp, no_w_full, no_sp + 1)) @@ -356,17 +358,21 @@ getTrophicLevel.MizerParams <- function(params, # Total consumption = (1 - f) * E, used for denominator accumulator consumption <- (1 - feeding_level) * encounter # no_sp x no_w - # Full predation kernel array (no_sp x no_w x no_w_full). - # getPredKernel() computes it from species parameters if not explicitly stored. - pred_kernel <- getPredKernel(params) - - # Prey-biomass weight K = w. When second-order, replace the left-edge - # value w_j by its trapezoidal bin-average. - w_ba <- bin_average_summary_weight(params@w, params) + # Full predation kernel array (no_sp x no_w x no_w_full). The numerator + # below is a trophic-level-weighted copy of the encounter integral, and its + # ratio to `consumption` is only a trophic level if the two use the same + # quadrature. encounter_kernel() therefore returns the kernel that + # mizerEncounter() effectively uses: the point-sampled kernel on the default + # path, and the prey-bin-integrated kernel when second-order bin-averaging + # is on (issue #474). + pred_kernel <- encounter_kernel(params) # prey_mass_tl[j, p] = N_j(w_p) * T_j(w_p) * w_p * dw_p + # The prey-biomass weight is the plain point weight `w * dw`, matching + # mizerEncounter(). Under second-order the prey-bin integral lives in + # `pred_kernel`, so bin-averaging this weight as well would apply it twice. # Initialised with T_j = 1 everywhere; updated as trophic levels are computed - prey_mass_tl <- sweep(n, 2, w_ba * params@dw, "*") # no_sp x no_w + prey_mass_tl <- sweep(n, 2, params@w * params@dw, "*") # no_sp x no_w # Size-dependent trophic level of the resource: # T_R(w) = max(1, 1 + log(w / w_R) / log(beta_R)) @@ -376,12 +382,9 @@ getTrophicLevel.MizerParams <- function(params, # to the numerator can be precomputed before the size loop. This mirrors the # `phi_prey_background` term in mizerEncounter(). tl_R <- pmax(1, 1 + log(params@w_full / w_R) / log(beta_R)) # no_w_full - # Full-grid summary weight, consistent with the species weight w_ba and with - # the resource weighting in mizerEncounter() (plain w_full on the default - # path; trapezoidal bin-average when second_order_w is on). - w_full_ba <- bin_average_summary_weight(params@w_full, params) - # TL-weighted resource biomass per bin. - resource_mass_tl <- tl_R * n_pp * w_full_ba * params@dw_full # no_w_full + # TL-weighted resource biomass per bin (length no_w_full), weighted exactly + # as the resource is weighted in mizerEncounter(). + resource_mass_tl <- tl_R * n_pp * params@w_full * params@dw_full # Resource contribution to the TL-weighted encounter for every predator/size: # ae_R[i, k] = sum_p kernel[i, k, p] * resource_mass_tl[p] ae_R <- rowSums(sweep(pred_kernel, 3, resource_mass_tl, "*"), dims = 2) @@ -427,7 +430,7 @@ getTrophicLevel.MizerParams <- function(params, tl_k <- tl[, k] tl_k[is.na(tl_k)] <- 1 prey_mass_tl[active_k, k] <- - n[active_k, k] * tl_k[active_k] * w_ba[k] * params@dw[k] + n[active_k, k] * tl_k[active_k] * params@w[k] * params@dw[k] } return(ArraySpeciesBySize(tl, value_name = "Trophic level", params = params)) diff --git a/man/encounter_kernel.Rd b/man/encounter_kernel.Rd new file mode 100644 index 000000000..aa5d7a3a1 --- /dev/null +++ b/man/encounter_kernel.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/setPredKernel.R +\name{encounter_kernel} +\alias{encounter_kernel} +\title{The predation kernel as used by the encounter quadrature} +\usage{ +encounter_kernel(params) +} +\arguments{ +\item{params}{A MizerParams object.} +} +\value{ +An array (predator species x predator size x prey size). +} +\description{ +Returns the kernel array \eqn{\Phi_i(w_k, w_p)} for which +\deqn{E_i(w_k) = \gamma_i(w_k) \sum_p \Phi_i(w_k, w_p) N^{eff}_i(w_p) + w_p \Delta w_p} +reproduces exactly the available energy computed by \code{\link[=mizerEncounter]{mizerEncounter()}}, +where \eqn{N^{eff}} is the interaction-weighted prey density. It is the +kernel that any summary function must use if its result is to be consistent +with \code{\link[=getEncounter]{getEncounter()}}. +} +\details{ +On the default first-order path this is just the point-sampled kernel +returned by \code{\link[=getPredKernel]{getPredKernel()}}. When second-order bin-averaging is switched on +(see \code{\link[=second_order_w]{second_order_w()}}) the two differ: \code{setPredKernel()} then builds the +Fourier-transformed kernel from the kernel \emph{integrated over the prey bin}, +divided by \eqn{\beta - 1} so that the plain point weight \eqn{w_p \Delta +w_p} carried by the prey vector is cancelled. Those bin-integrated weights +are recovered here from \code{params@ft_pred_kernel_e} by an inverse Fourier +transform, which costs one FFT and keeps this helper automatically in step +with whatever quadrature \code{setPredKernel()} used. + +A summary function that instead pairs the point-sampled \code{\link[=getPredKernel]{getPredKernel()}} +with a bin-averaged prey weight double-counts the prey-bin quadrature; that +was the bug behind issue #474. +} +\concept{helper} +\keyword{internal} diff --git a/man/expand_kernel_offsets.Rd b/man/expand_kernel_offsets.Rd new file mode 100644 index 000000000..564ea8466 --- /dev/null +++ b/man/expand_kernel_offsets.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/setPredKernel.R +\name{expand_kernel_offsets} +\alias{expand_kernel_offsets} +\title{Expand kernel weights indexed by grid offset into a full kernel array} +\usage{ +expand_kernel_offsets(phis, params, species) +} +\arguments{ +\item{phis}{A species-by-offset matrix of kernel weights, with the offset +running from 0 to \code{length(params@w_full) - 1}.} + +\item{params}{A MizerParams object supplying the grid.} + +\item{species}{A character vector of species names for the dimnames.} +} +\value{ +An array (predator species x predator size x prey size). +} +\description{ +The predation kernel depends only on the predator/prey mass ratio, so on the +geometric grid it is a function of the offset \eqn{m} between the predator +and prey grid indices alone. \code{phis[i, m + 1]} holds the weight of species +\eqn{i} at offset \eqn{m}. This helper writes those weights into the +(predator species x predator size x prey size) array that the non-FFT code +paths work with. +} +\concept{helper} +\keyword{internal} diff --git a/man/getDiet.Rd b/man/getDiet.Rd index 673008563..390612acf 100644 --- a/man/getDiet.Rd +++ b/man/getDiet.Rd @@ -71,7 +71,10 @@ the rate of consumption of biomass from these extra components. This function performs the same integration as \code{\link[=getEncounter]{getEncounter()}} but does not aggregate over prey species, and multiplies by \eqn{1-f_i(w)} to get the consumed biomass rather than the available biomass. Outside the range of -sizes for a predator species the returned rate is zero. +sizes for a predator species the returned rate is zero. Summing the result +of \code{getDiet(proportion = FALSE)} over prey therefore reproduces +\code{getEncounter(params) * (1 - getFeedingLevel(params))}, whichever quadrature +scheme the model uses (see \code{\link[=second_order_w]{second_order_w()}}). } \examples{ diet <- getDiet(NS_params) diff --git a/tests/testthat/test-second_order_summary.R b/tests/testthat/test-second_order_summary.R index 75ce2b37d..801075300 100644 --- a/tests/testthat/test-second_order_summary.R +++ b/tests/testthat/test-second_order_summary.R @@ -130,3 +130,83 @@ test_that("second-order biomass converges to default as the grid is refined", { } expect_lt(rel_diff(fine), rel_diff(coarse)) }) + + +# Diet and trophic level must use the encounter quadrature (issue #474) ---- + +# When second-order bin-averaging is on, the prey-bin integral is folded into +# the Fourier-transformed kernel by setPredKernel(). A summary function that +# also bin-averages its prey weight applies that quadrature twice and inflates +# the result by (1 + beta) / 2. These tests pin the identities that catch it. + +custom_kernel_params <- function(params) { + pk <- getPredKernel(params) + comment(pk) <- "set manually" + setPredKernel(params, pred_kernel = pk) +} + +test_that("encounter_kernel reproduces the encounter rate in both modes", { + check <- function(params) { + no_w <- length(params@w) + no_w_full <- length(params@w_full) + idx_sp <- (no_w_full - no_w + 1):no_w_full + n <- initialN(params) + n_pp <- initialNResource(params) + kernel <- encounter_kernel(params) + n_eff <- sweep(params@interaction %*% n, 2, params@w * params@dw, "*") + species <- rowSums(sweep(kernel[, , idx_sp, drop = FALSE], c(1, 3), + n_eff, "*"), dims = 2) + resource <- params@species_params$interaction_resource * + rowSums(sweep(kernel, 3, params@w_full * params@dw_full * n_pp, + "*"), dims = 2) + direct <- params@search_vol * (species + resource) + + params@ext_encounter + expect_equal(as.vector(direct), as.vector(getEncounter(params)), + tolerance = 1e-4) + } + p <- NS_params_small + check(p) + check(custom_kernel_params(p)) + second_order_w(p) <- c(bin_average = TRUE) + check(p) + check(custom_kernel_params(p)) +}) + +test_that("getDiet summed over prey equals the consumption rate in both modes", { + check <- function(params) { + total <- rowSums(getDiet(params, proportion = FALSE), dims = 2) + consumption <- getEncounter(params) * (1 - getFeedingLevel(params)) + # Outside a species' size range the diet is set to zero, so compare + # only where the abundance is positive. + mask <- initialN(params) > 0 + expect_equal(as.vector(total)[mask], as.vector(consumption)[mask]) + } + p <- NS_params_small + check(p) + check(custom_kernel_params(p)) + second_order_w(p) <- c(bin_average = TRUE) + check(p) + check(custom_kernel_params(p)) +}) + +test_that("getTrophicLevel gives 2 for a predator whose prey all have level 1", { + # At the smallest size on the grid no consumer has yet been assigned a + # trophic level above 1, and choosing w_R above the whole grid makes the + # resource trophic level 1 as well. The trophic-level-weighted encounter in + # the numerator then equals the plain encounter in the denominator, so the + # trophic level must come out as exactly 2. It does so only if numerator and + # denominator use the same quadrature. + check <- function(params) { + w_R <- max(params@w_full) * 10 + tl <- getTrophicLevel(params, w_R = w_R) + first <- which(params@w_min_idx == 1) + expect_equal(as.vector(tl[first, 1]), rep(2, length(first)), + tolerance = 1e-4) + } + p <- NS_params_small + check(p) + check(custom_kernel_params(p)) + second_order_w(p) <- c(bin_average = TRUE) + check(p) + check(custom_kernel_params(p)) +}) From 2fc672611356024a8f97ebab44ac0b0601899c48 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Mon, 3 Aug 2026 16:07:24 +0200 Subject: [PATCH 2/2] Document the two size-grid quadrature schemes The `second_order_w` flag was documented in the FFT vignette from the kernel's side, but nothing said where each bin integral is performed across the model as a whole. That gap is what let #474 through: the prey-bin quadrature is folded into the predation kernel at setup, so a summary function that bin-averages its prey weight as well applies it twice, and the resulting uniform (1 + beta) / 2 factor cancels in any proportion. numerical_details.qmd gains: - "Applying each bin integral exactly once", inside the existing "Point values and bin averages" section. Derives the encounter kernel weight, explains why the `w_p dw_p` in the run-time sum is a normalisation rather than a first-order quadrature weight, and tabulates which factor of each integral is bin-integrated and in which setter. - "The `second_order_w` switch", a new section: the two independent entries, a full FALSE-vs-TRUE table (including the rows that stay point values), and a recipe for making your own diagnostic second-order accurate, ending with the identity check and how to read off (1 + beta) / 2 versus 2 / (1 + beta) when it fails. It also drops a paragraph that had been duplicated verbatim. .claude/skills/size-grid-integrals.md turns this into a decision procedure for anything that integrates over the size grid: four cases, the traps (including both #474 failures), the verification steps, and the helper table. AGENTS.md gets the architectural note that the model has two discretisations and a pointer to the skill. Co-Authored-By: Claude Opus 5 --- .claude/skills/size-grid-integrals.md | 138 ++++++++++++++++++++++++++ AGENTS.md | 3 + NEWS.md | 7 ++ vignettes/numerical_details.qmd | 131 +++++++++++++++++++++++- 4 files changed, 277 insertions(+), 2 deletions(-) create mode 100644 .claude/skills/size-grid-integrals.md diff --git a/.claude/skills/size-grid-integrals.md b/.claude/skills/size-grid-integrals.md new file mode 100644 index 000000000..ce530e58b --- /dev/null +++ b/.claude/skills/size-grid-integrals.md @@ -0,0 +1,138 @@ +# Integrating over the size grid + +Use this skill when writing anything that integrates over the size grid: a +summary or indicator function, a diagnostic derived from a rate, a new rate +setter with a size-dependent parameter, or an extension that reaches inside the +predation convolution. + +mizer has **two quadrature schemes**, selected by the `bin_average` entry of the +`second_order_w` slot. `FALSE` is the default, so code that ignores the flag +looks correct, passes its tests, and is silently wrong by ~10% for every user +who has switched second order on. + +## The one rule + +Each bin integral is performed in exactly one place. A size-dependent factor is +bin-averaged **where its integral is performed, and nowhere else**. + +The theory — why, and where each of mizer's integrals is performed — is in +`vignettes/numerical_details.qmd`, sections *Point values and bin averages* and +*The `second_order_w` switch*. Read those rather than re-deriving; the tables +there are the authoritative inventory and are meant to be kept current. + +## Decide which case you are in + +### Case 1 — a plain integral against the abundance, ∫ K(w) N(w) dw + +Discretise as `sum_j Kbar_j * N_j * dw_j`. Only `K` is approximated: + +```r +K <- bin_average_summary_weight(K, params) # gated on the flag +drop(n %*% (K * params@dw)) +``` + +- **Gate it, don't hard-code it.** `bin_average_summary_weight()` returns `K` + untouched on the default path, so the old numbers stay byte-identical. +- **Never bin-average `N` or `dw`.** `N_j` is already a cell average and `dw_j` + is exact. +- **Average the product, not the factors.** SSB averages `psi * w`; yield + averages `F * w`. Averaging separately is a different (wrong) number. +- **If `K` is an exact power law `w^a`,** use `power_law_bin_average(w, dw, a)` + instead of the trapezoid — it is exact, not merely second order. +- If the result is size-resolved, tag it: `ArraySpeciesBySize(..., representation + = "average")` for a bin average, `"point"` for a boundary quantity. The tag + drives the half-bin plotting shift. + +### Case 2 — a quantity built from rates mizer already computes + +Call `getEncounter()`, `getFeedingLevel()`, `getPredRate()`, `getEGrowth()`, … +and do not rebuild them. The rate functions already carry the right quadrature +for whichever scheme the model is in. Re-deriving a rate is how the two known +bugs in this area were introduced. + +### Case 3 — you need to go inside the encounter or predation convolution + +Use **`encounter_kernel(params)`**, not `getPredKernel(params)`, and pair it +with the **plain point weight** `params@w * params@dw`. + +Under `bin_average`, `setPredKernel()` builds `ft_pred_kernel_e` from the kernel +integrated over the prey bin and divides by `beta - 1` precisely so that the +`w * dw` supplied by the prey vector cancels. That `w * dw` is a normalisation, +not a first-order quadrature weight — bin-averaging it applies the prey-bin +integral twice. + +### Case 4 — a new rate setter with a size-dependent parameter + +Follow `setExtMort()` / `setExtDiffusion()` / `setResource()`: gate on the flag +and use `power_law_bin_average()` for power laws, or a composite midpoint rule +(as `setFishing()` does for selectivity) for anything else. Do the integral once, +at setup, so the projection cost is unchanged. + +## Traps + +### Double-counting is a uniform factor, so normalised outputs hide it + +On a geometric grid `bin_average_weight(w) / w` is exactly `(1 + beta) / 2` +(1.0967 for `NS_params`). Applying the prey-bin quadrature twice therefore +scales the result by a constant, which **cancels in any proportion or ratio**. +`getDiet(proportion = FALSE)` was 9.7% too large for a long time while the +default `proportion = TRUE` stayed correct (#474). If a consistency ratio comes +out as a constant, read off its value: `(1 + beta) / 2` means the quadrature was +applied twice, `2 / (1 + beta)` means it is missing. + +### `getPredKernel()` is not the kernel the encounter uses + +It returns the kernel point-sampled on the grid — right for plotting, and the +form you supply a custom kernel in, but not the bin-integrated coefficients the +convolution consumes. Pairing it with `getEncounter()` in a numerator/denominator +is what made `getTrophicLevel()` wrong (#474). + +### Growth-type rates are never bin-averaged + +`g`, `e`, the encounter rate and the feeding level are point values at `w_j` +under **both** settings — they are boundary velocities. What improves them when +the flag is on is the encounter integral behind them, not any averaging of the +rate itself. Bin-averaging them is an error, not an upgrade. + +### Testing only the default path proves nothing + +Both #474 bugs were invisible with `bin_average = FALSE`. Every new integral +needs a test with the flag on. + +## Verifying a change + +1. **Default path unchanged.** Assert byte-identity (or an existing snapshot) + with `bin_average = FALSE`. Any movement there is a regression. +2. **Flag on: assert the identity your quantity should satisfy.** Anything that + decomposes a rate must reassemble into it: + + ```r + params <- NS_params_small + second_order_w(params) <- c(bin_average = TRUE) + total <- rowSums(getDiet(params, proportion = FALSE), dims = 2) + ratio <- total / (getEncounter(params) * (1 - getFeedingLevel(params))) + range(ratio[initialN(params) > 0]) # 1 1 + ``` + + FFT convolution is circular, so allow ~1e-4 when comparing a direct sum + against a rate function; within one code path the agreement is exact. +3. **Convergence.** The gap between the two schemes should shrink under grid + refinement — see the "second-order biomass converges to default" test for the + pattern. + +New tests go in `tests/testthat/test-second_order_summary.R`, using the +`NS_params_small` fixture and toggling with +`second_order_w(p) <- c(bin_average = TRUE)`. + +## Helpers + +| Helper | Use | +|---|---| +| `bin_average_summary_weight(K, params)` | trapezoidal bin average, gated on the flag — the default entry point | +| `bin_average_weight(K)` | ungated trapezoid; averages along the last dimension of an array | +| `power_law_bin_average(w, dw, a, w_max)` | exact bin average of `w^a`, with optional cutoff | +| `encounter_kernel(params)` | the kernel `mizerEncounter()` actually uses, under either scheme | +| `bin_midpoints(params)` | geometric bin centres, for plotting bin averages | + +Setting `second_order_w(params) <- c(bin_average = TRUE)` re-runs `setParams()`, +because every array in the inventory table is precomputed. diff --git a/AGENTS.md b/AGENTS.md index 4344832f7..303bcd9ac 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -14,6 +14,8 @@ mizer is an R package for dynamic multi-species size-spectrum modelling of fish **Customisable rate functions**: users replace rate functions by storing a custom function name in `params@rates_funcs`. Dispatch via `get(params@rates_funcs$FunctionName)(params, ...)`. +**Two quadrature schemes**: the `second_order_w` slot selects how the model is discretised on the size grid — `flux` picks the advective reconstruction, and `bin_average` decides whether size-dependent factors are integrated over their bin or point-sampled at the left bin boundary. Both default to the first-order scheme, so previous mizer versions reproduce byte-for-byte. This is invisible in the code you are likely to be reading: a function that ignores `bin_average` still looks correct and still passes its tests on the default path. Anything that integrates over the size grid must handle both schemes and be tested under both. + **Auto-generated files** — never edit `NAMESPACE`, `man/`, `RcppExports.R`, or `RcppExports.cpp` directly. ## Code Conventions @@ -23,6 +25,7 @@ mizer is an R package for dynamic multi-species size-spectrum modelling of fish - **Language**: British English (en-GB) — "colour", "behaviour", "modelling" - When documenting a mizer S3 generic whose methods share a man page (combined with `@rdname`/`@name`), follow the steps in `.claude/skills/document-s3-generics.md`. - When adding, moving or removing a species parameter default, follow `.claude/skills/species-param-defaults.md`. A default belongs to the rate setter that reads the parameter; only parameters that no single rate setter owns are defaulted centrally. +- When writing anything that integrates over the size grid — a summary or indicator function, a diagnostic derived from a rate, a rate setter with a size-dependent parameter — follow `.claude/skills/size-grid-integrals.md`. Each bin integral is performed in exactly one place, so a size-dependent factor is bin-averaged where its integral is performed and nowhere else; doing it twice is a silent uniform error. ## Testing diff --git a/NEWS.md b/NEWS.md index 6b074d2d4..0fded2fd7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -85,6 +85,13 @@ stability of steady states. and denominator now use the same quadrature, and a predator whose prey all have trophic level 1 comes out at exactly 2 in both schemes (#474). +## Documentation + +- The "Point values and bin averages" section of `vignette("numerical_details")` + now explains where each bin integral is performed and why it must be applied + exactly once, and a new "The `second_order_w` switch" section documents what + the flag changes and how to make your own diagnostic second-order accurate. + # mizer 3.2.1 This patch release fixes how species and gear parameters are handled when they diff --git a/vignettes/numerical_details.qmd b/vignettes/numerical_details.qmd index 0d886b372..2a73c3d22 100644 --- a/vignettes/numerical_details.qmd +++ b/vignettes/numerical_details.qmd @@ -102,12 +102,139 @@ In summary: | $d$, $\mu$; fishing/reproductive investments; predation and encounter integrands | bin properties (a coefficient inside $\partial(dN)/\partial w$, or a rate integrated against $N$ over a bin) | bin average over $[w_j, w_{j+1}]$ | | $e$ (energy income) | both: growth velocity *and* reproduction integrand | point value for growth; bin-averaged product $\psi e$ for reproduction | -**Plotting follows the same distinction.** A bin average $N_j$ does not live at the bin boundary $w_j$ but at the geometric bin centre $w^*_j=\sqrt{w_j\,w_{j+1}}=w_j\sqrt\beta$ (the log-midpoint, exact for the community spectrum $N\propto w^{-2}$). So under second-order bin-averaging mizer draws bin-averaged quantities (the abundance and the mortality/reproduction sinks) at $w^*_j$ — a uniform half-bin shift to the right on the log axis — while point-valued quantities (the encounter and growth-type rates) stay on the nodes $w_j$. The size-resolved array classes carry a `representation` tag recording which a quantity is, and the shift is applied only when `second_order_w[["bin_average"]]` is set, so default plots are unchanged. +### Applying each bin integral exactly once {#sec-one-quadrature} -**Plotting follows the same distinction.** A bin average $N_j$ does not live at the bin boundary $w_j$ but at the geometric bin centre $w^*_j=\sqrt{w_j\,w_{j+1}}=w_j\sqrt\beta$ (the log-midpoint, exact for the community spectrum $N\propto w^{-2}$). So under second-order bin-averaging mizer draws bin-averaged quantities (the abundance and the mortality/reproduction sinks) at $w^*_j$ — a uniform half-bin shift to the right on the log axis — while point-valued quantities (the encounter and growth-type rates) stay on the nodes $w_j$. The size-resolved array classes carry a `representation` tag recording which a quantity is, and the shift is applied only when `second_order_w[["bin_average"]]` is set, so default plots are unchanged. +Knowing *that* a factor must be bin-averaged is only half of the rule; the other half is *where*. Each bin integral is performed in exactly one place, and a quantity assembled out of parts must not repeat an integral that one of its parts has already done. + +The encounter rate is the clearest example. Continuously, +$$ +E_i(w) = \gamma_i(w)\int \tilde\phi_i\!\left(\frac{w}{w_p}\right) N^{\text{eff}}_i(w_p)\, w_p\, dw_p , +$$ {#eq-encounter-integral} +where $N^{\text{eff}}_i=\sum_j\theta_{ij}N_j+\theta_{iR}N_R$ is the interaction-weighted prey density. The smooth weight multiplying that density is $K(w_p)=\tilde\phi_i(w/w_p)\,w_p$ — the kernel **and** the mass factor $w_p$ together — and it is $K$ that has to be integrated over the prey bin. Mizer does this once, at setup, in `setPredKernel()`, which stores the kernel coefficient for the grid offset $m$ as +$$ +\Phi^E_i[m] = \frac{\ln\beta}{\beta-1}\int_0^1 \tilde\phi_i(\beta^{\,m-s})\;\beta^{2s}\,ds . +$$ {#eq-encounter-kernel-weight} +The Jacobian $\beta^{2s}$ carries the $w_p\,dw_p$ of the integrand, and the division by $\beta-1$ *removes* the factor $w_p\,\Delta w_p=(\beta-1)\,w_p^2$ that the prey vector will supply. At run time the rate function therefore evaluates the plain sum +$$ +E_i(w_k) = \gamma_i(w_k) \sum_p \Phi^E_i[k-p]\; N^{\text{eff}}_{i,p}\; w_p\,\Delta w_p , +$$ {#eq-encounter-sum} +with the **point** value $w_p$ at the bin boundary. (The [FFT vignette](fft.html) derives @eq-encounter-kernel-weight and its predation and predation-diffusion analogues.) + +The factor $w_p\,\Delta w_p$ in @eq-encounter-sum is thus not a quadrature weight that has been left at first order; it is a normalisation the kernel has already divided out. Replacing it by the bin-averaged $\bar w_p\,\Delta w_p$ would apply the prey-bin integral a second time and inflate every encounter by +$$ +\frac{\bar w_p}{w_p} = \frac{w_p+w_{p+1}}{2\,w_p} = \frac{1+\beta}{2}, +$$ {#eq-double-count} +uniformly across the grid — 9.7 % for the North Sea model, where $\beta=1.1934$. That is not hypothetical: it is exactly the error `getDiet(proportion = FALSE)` made (issue #474). It went unnoticed for a while precisely because the factor is uniform and so cancels in the default `proportion = TRUE` normalisation. + +The table below records which factor of each integral is bin-integrated, and where: + +| Integral | Factor that is bin-integrated | Performed in | +|---|---|---| +| encounter, predation and predation-diffusion convolutions | the kernel weight, $\tilde\phi\,w_p$, $\tilde\phi$ and $\tilde\phi\,w_p^2$ respectively | `setPredKernel()`, into `ft_pred_kernel_e`/`_p`/`_d` | +| fishing sink $\int Q\,S(w)\,\text{effort}\,N\,dw$ | the selectivity $S$ | `setFishing()`, into `selectivity` | +| external mortality and external diffusion sinks | the power laws $z_{ext}w^{d}$ and $D_{ext}w^{n+1}$ | `setExtMort()`, `setExtDiffusion()` | +| resource semichemostat terms | the power laws $r_{pp}w^{n-1}$ and $\kappa w^{-\lambda}$ | `setResource()`, into `rr_pp` and `cc_pp` | +| reproduction $\int \psi(w)\,e(w)\,N(w)\,dw$ | the product $\psi e$ | `mizerRDI()` | +| summary integrals $\int K(w)\,N(w)\,dw$ | the weight $K$: $w$ for biomass, $\psi w$ for SSB, $F w$ for yield | `getBiomass()`, `getSSB()`, `getYield()`, … | + +Two corollaries are worth stating explicitly. + +**Never bin-average the density or the bin width.** $N_j$ is already a bin average (@eq-bin-average) and $\Delta w_j$ is already exact. It is only the smooth weight multiplying them that is being approximated, so it is the only thing that gets averaged. + +**A diagnostic that decomposes a rate must borrow that rate's quadrature rather than rebuild it.** `getDiet()` is `getEncounter()` resolved by prey species, so it uses the same kernel and the same point prey weight $w_p\,\Delta w_p$; summed over prey it then reproduces $(1-f_i(w))\,E_i(w)$ exactly, under both schemes. `getTrophicLevel()` is a *ratio* of a trophic-level-weighted encounter to the plain encounter, so its numerator and denominator must be built from the same kernel. The internal helper `encounter_kernel()` exists for this: it returns the kernel that `mizerEncounter()` is actually using — the point-sampled kernel by default, the bin-integrated one when `bin_average` is on, and the stored array when a custom kernel has been set. + +Note that `getPredKernel()` is **not** that kernel when `bin_average` is on. It returns $\tilde\phi_i$ point-sampled on the grid, which is the right object for plotting or inspecting a feeding kernel and the form in which you supply a custom kernel, but it is not the bin-integrated coefficient @eq-encounter-kernel-weight that the convolution consumes. + +### Plotting follows the same distinction + +A bin average $N_j$ does not live at the bin boundary $w_j$ but at the geometric bin centre $w^*_j=\sqrt{w_j\,w_{j+1}}=w_j\sqrt\beta$ (the log-midpoint, exact for the community spectrum $N\propto w^{-2}$). So under second-order bin-averaging mizer draws bin-averaged quantities (the abundance and the mortality/reproduction sinks) at $w^*_j$ — a uniform half-bin shift to the right on the log axis — while point-valued quantities (the encounter and growth-type rates) stay on the nodes $w_j$. The size-resolved array classes carry a `representation` tag recording which a quantity is, and the shift is applied only when `second_order_w[["bin_average"]]` is set, so default plots are unchanged. For the `power`-weighted spectrum plots (`plotSpectra()` and friends) the $w^{\text{power}}$ factor must be evaluated where the density value lives, so it too is taken at the bin centre: each marker is the point $\bigl(w^*_j,\,N_j\,(w^*_j)^{\text{power}}\bigr)$ on the continuous $N(w)\,w^{\text{power}}$ curve. (Sampling the weight at the edge would mis-scale it by a factor $\beta^{\text{power}/2}$, largest for the common $\text{power}=2$ Sheldon plot.) A cumulative plot (`plotCDF()`) is the opposite case: a CDF value is cumulative *up to a size*, a boundary quantity, so its increments use the bin-averaged (centre-weighted) density but the cumulative is plotted on the bin **edges**, not the centres. Because the cumulative sum is inclusive — the sum through bin $k$ is the integral over all bins up to and including bin $k$ — each cumulative value is placed on that bin's *upper* edge $w_k+\Delta w_k$ (in both the default and second-order schemes). This makes the inclusive convention explicit and removes a one-bin offset that would otherwise leave the CDF only first-order accurate in its placement. +## The `second_order_w` switch {#sec-second-order-switch} + +Everything in @sec-point-values describes the second-order scheme, but mizer does not use it by default: the historical first-order behaviour is preserved so that existing models reproduce their published results exactly. The choice lives in the `second_order_w` slot, which has two independent entries. + +* **`flux`** selects the reconstruction of the density at the bin boundary in the advective flux: `"upwind"` (first order, the default), `"van_leer"` (second order, limited, keeps abundances non-negative) or `"centred"` (second order, unlimited). This affects only the transport step; see @sec-reducing-spatial-error. +* **`bin_average`** is a logical flag selecting whether the size-dependent factors listed in @sec-one-quadrature are integrated over their bin (`TRUE`) or point-sampled at the left bin boundary $w_j$ (`FALSE`, the default). + +```r +second_order_w(params) <- TRUE # both: van_leer + bin averaging +second_order_w(params) <- "centred" # flux scheme only +second_order_w(params) <- c(bin_average = TRUE) # bin averaging only +second_order_w(params) <- FALSE # back to the mizer defaults +``` + +The two are independent because they correct different errors: `flux` improves the time evolution of the spectrum, `bin_average` improves the rates that drive it, and either alone leaves the other at first order. Changing `bin_average` re-runs `setParams()`, because all the arrays in the first table of @sec-one-quadrature are precomputed and have to be rebuilt. + +### What `bin_average` changes + +Write $\beta=w_{j+1}/w_j$ for the (constant) grid ratio. On a geometric grid the exact bin average of a power law $w^{a}$ is +$$ +\overline{w^{a}}_j = \frac{1}{\Delta w_j}\int_{w_j}^{w_{j+1}} w^{a}\,dw + = \frac{w_{j+1}^{\,a+1}-w_j^{\,a+1}}{(a+1)\,\Delta w_j} + = w_j^{\,a}\,\frac{\beta^{\,a+1}-1}{(a+1)(\beta-1)}, +$$ {#eq-power-law-bin-average} +which mizer computes with the internal helper `power_law_bin_average()` and uses wherever the factor is genuinely a power law. Where it is not, the trapezoidal average $\bar K_j=\tfrac12(K_j+K_{j+1})$ is used, which is second order for any smooth $K$ and exact for $K$ linear in $w$ (`bin_average_weight()`). Where the factor is a kernel the bin integral is done by composite quadrature at setup, @eq-encounter-kernel-weight. + +| Quantity | `bin_average = FALSE` | `bin_average = TRUE` | +|---|---|---| +| kernel coefficients `ft_pred_kernel_e`, `_p`, `_d` | $\tilde\phi_i(\beta^{m})$, point-sampled | bin-integrated, @eq-encounter-kernel-weight | +| gear selectivity | $S(w_j)$ | $\frac{1}{\Delta w_j}\int S\,dw$, by composite midpoint | +| external mortality $z_{ext}w^{d}$, external diffusion $D_{ext}w^{n+1}$ | point value at $w_j$ | exact bin average, @eq-power-law-bin-average | +| resource rate $r_{pp}w^{n-1}$, capacity and initial spectrum $\kappa w^{-\lambda}$ | point value at $w_j$ | exact bin average, @eq-power-law-bin-average | +| reproduction integrand $\psi\,e$ | point value at $w_j$ | trapezoidal bin average | +| summary weights in `getBiomass()`, `getSSB()`, `getYield()`, `getYieldGear()` | $K(w_j)$ | trapezoidal bin average of $K$ | +| plotting position of a bin-averaged quantity | node $w_j$ | bin centre $w_j\sqrt\beta$ | +| growth rate $g$, encounter rate $E$, feeding level $f$ | point value at $w_j$ | **unchanged** — point value at $w_j$ | + +The last two rows are the ones that are easy to get wrong. `getN()` is also unchanged, because its weight is $K\equiv 1$ and the bin average of a constant is the constant. And the growth-type rates stay point values under both settings: they are boundary velocities (@sec-point-values), so bin-averaging them would be an error, not an improvement — what improves them under `bin_average` is that the encounter integral feeding them is now second order, not any averaging of $g$ itself. + +Because one setting is first order and the other second, the difference between them is itself $O(\Delta x)$ and is a usable estimate of the discretisation error of the default scheme: if flipping the flag moves a result by more than you are willing to tolerate, the grid is too coarse for that result. + +### Making your own quantity second-order accurate {#sec-second-order-recipe} + +If you compute a diagnostic of your own — in an extension package, or in analysis code on top of a `MizerParams` object — it falls into one of two cases. + +**Case 1: a plain integral against the abundance,** $\int K(w)\,N(w)\,dw$. Discretise it as $\sum_j \bar K_j\,N_j\,\Delta w_j$: keep $N_j$ and $\Delta w_j$ exactly as they are and replace the point weight $K(w_j)$ by its bin average. + +```r +# Trapezoidal bin average of a weight along the size axis. The top bin has no +# right neighbour, so it is left one-sided; the density there is negligible. +bin_average <- function(K) { + n <- length(K) + c(0.5 * (K[-n] + K[-1]), K[n]) +} + +my_indicator <- function(params, K) { + if (isTRUE(second_order_w(params)$bin_average)) K <- bin_average(K) + drop(initialN(params) %*% (K * dw(params))) +} + +# biomass above 10 g, say +params <- NS_params +K <- w(params) * (w(params) >= 10) +my_indicator(params, K) +``` + +Gate the averaging on the flag, as above, so that your diagnostic follows the model it is given rather than silently disagreeing with `getBiomass()`. If $K$ is an exact power law, use @eq-power-law-bin-average instead of the trapezoid: it is exact rather than merely second order. If $K$ is a product of a size-dependent rate and a mass factor — $\psi(w)\,w$, or $F(w)\,w$ — average the **product**, not the factors separately. + +**Case 2: a quantity built from rates mizer already computes.** Get the rates from the rate functions (`getEncounter()`, `getFeedingLevel()`, `getPredRate()`, `getEGrowth()`, …) and do not re-derive them, because the rate functions already carry the correct quadrature for the current setting. In particular, do not rebuild an encounter or predation convolution out of `getPredKernel()` and a hand-written prey weight: under `bin_average` that kernel is the point-sampled one and will not agree with the rate function. If you genuinely need the resolved convolution — as `getDiet()` and `getTrophicLevel()` do — pair mizer's own encounter kernel with the plain point prey weight $w_p\,\Delta w_p$, and never with a bin-averaged one. + +**Check the result against an identity.** Any diagnostic that decomposes a rate should reassemble into it. For a diet-like decomposition: + +```r +params <- NS_params +second_order_w(params) <- TRUE +total <- rowSums(getDiet(params, proportion = FALSE), dims = 2) +ratio <- total / (getEncounter(params) * (1 - getFeedingLevel(params))) +range(ratio[initialN(params) > 0]) +#> 1 1 +``` + +If such a ratio comes out as a constant instead of 1, read off its value: $(1+\beta)/2$ means the prey-bin quadrature has been applied twice, and $2/(1+\beta)$ means it is missing. The value is easy to recognise, since $\beta$ is just `w_full(params)[2] / w_full(params)[1]`. + + ## Semi-Implicit Time Discretisation With the diffusion term, an explicit time discretisation would require a very small time step for stability ($\Delta t \sim \Delta w^2$). Therefore, we use a semi-implicit scheme where the densities $N$ are evaluated at time $t+1$, but the rates ($g$, $\mu$, $d$) are evaluated at time $t$. Using a fully implicit scheme would require solving a nonlinear system at each time step, which is more computationally expensive.