Starting to add diffusion

.gitignore

@@ -9,3 +9,4 @@ docs/dev
 .vscode
 .positai
 .claude
+.codex

CLAUDE.md

@@ -0,0 +1,63 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## About mizer
+
+mizer is an R package for dynamic multi-species size-spectrum modelling of fish communities. It models marine ecosystems subject to fishing, tracking individual fish growth from egg size to maximum size and capturing ontogenetic diet shifts.
+
+## Common Commands
+
+```r
+devtools::load_all()        # Load package for development
+devtools::document()        # Regenerate NAMESPACE and man/ pages from roxygen2
+devtools::test()            # Run all tests
+devtools::check()           # Full R CMD check
+lintr::lint_package()       # Lint the package
+
+# Run a single test file
+testthat::test_file("tests/testthat/test-filename.R")
+
+# After editing C++ source
+devtools::clean_dll(); devtools::load_all()
+```
+
+## Architecture
+
+### Core Classes
+
+**`MizerParams`** (S4, `R/MizerParams-class.R`) — the central object passed to nearly all functions. Holds all model configuration: species parameters, size grids (`w`, `w_full`), interaction matrices, gear selectivity, rate function overrides (`@rates_funcs`), and resource dynamics. Validated by `validMizerParams()`. Modified via setter functions that return new copies: `setFishing(params, ...)`, `setInteraction(params, ...)`, etc.
+
+**`MizerSim`** (S4, `R/MizerSim-class.R`) — stores simulation output: a 3D array `n` (time × species × size), `n_pp` (time × size) for resource, `n_other` for additional biomass components, `effort` history, and the `MizerParams` used.
+
+**`MizerRate`** (S3, `R/MizerRate-class.R`) — wraps 2D arrays (species × size) returned by rate functions with metadata (`rate_name`, `units`). Inherits from `matrix`/`array`. Provides enhanced `print()`, `summary()`, `plot()`, and `as.data.frame()`.
+
+### Execution Flow
+
+1. **Setup**: `newMultispeciesParams()` / `newSingleSpeciesParams()` / `newCommunityParams()`
+2. **Configure**: `set*()` functions (`setFishing()`, `setPredKernel()`, `setMetabolicRate()`, …)
+3. **Rates**: `getEncounter()`, `getFeedingLevel()`, `getPredMort()`, `getFMort()`, `getRates()` — each returns a `MizerRate` matrix (species × size)
+4. **Project**: `project(params, t_max = 100, effort = ...)` → `MizerSim`
+5. **Analyse**: `getYield()`, `getBiomass()`, `getSSB()`, `plotSpectra()`, etc.
+
+### Customisable Rate Functions
+
+Users can replace any rate function by storing a custom function name in `params@rates_funcs`. Calls dispatch via `get(params@rates_funcs$FunctionName)(params, ...)`. This is the primary extensibility mechanism.
+
+### C++ Integration
+
+Performance-critical inner projection loop lives in `src/inner_project_loop.cpp` and `src/project_n_loop.cpp`. `RcppExports.R` and `RcppExports.cpp` are auto-generated — never edit them directly.
+
+### Extensibility via "Other" Components
+
+Arbitrary biomass components beyond species can be added (e.g., detritus, zooplankton). Stored in `n_other`; custom rate functions can be registered for them.
+
+## Code Conventions
+
+- **Indentation**: 4 spaces
+- **Naming**: camelCase or snake_case for functions/variables; PascalCase for classes
+- **Language**: British English (en-GB) — "colour", "behaviour", "modelling"
+- **Documentation**: All exported functions require roxygen2 with `@param`, `@return`, `@export`; use `@seealso` for cross-references
+- **Tests**: testthat edition 3; use `expect_doppelganger()` (vdiffr) for plot tests, snapshot tests for complex outputs
+- After adding exports, run `devtools::document()` to regenerate `NAMESPACE`
+- Update `NEWS.md` when adding features or fixing bugs

DESCRIPTION

@@ -58,6 +58,7 @@ Collate:
     'reproduction.R'
     'saveParams.R'
     'species_params.R'
+    'getRequiredRDD.R'
     'setColours.R'
     'setInteraction.R'
     'setPredKernel.R'
@@ -79,6 +80,8 @@ Collate:
     'resource_dynamics.R'
     'resource_semichemostat.R'
     'resource_logistic.R'
+    'numerical_methods.R'
+    'transport.R'
     'project_n.R'
     'project.R'
     'mizer-package.R'

NAMESPACE

@@ -62,6 +62,7 @@ S3method(getFMortGear,MizerParams)
 S3method(getFMortGear,MizerSim)
 S3method(getFeedingLevel,MizerParams)
 S3method(getFeedingLevel,MizerSim)
+S3method(getFlux,MizerParams)
 S3method(getGrowthCurves,MizerParams)
 S3method(getGrowthCurves,MizerSim)
 S3method(getInitialEffort,MizerParams)
@@ -254,6 +255,7 @@ export(getExtMort)
 export(getFMort)
 export(getFMortGear)
 export(getFeedingLevel)
+export(getFlux)
 export(getGrowthCurves)
 export(getInitialEffort)
 export(getInteraction)
@@ -280,6 +282,7 @@ export(getRateFunction)
 export(getRates)
 export(getReproductionLevel)
 export(getReproductionProportion)
+export(getRequiredRDD)
 export(getResourceCapacity)
 export(getResourceDynamics)
 export(getResourceLevel)

NEWS.md

@@ -1,11 +1,46 @@
-# Development version 2.5.4.9000
+# Development version 2.5.4.9101
 
 - `t_max` and `t_save` arguments in `project()` are now respected even when an
   effort array is supplied. When `t_max` is provided, the simulation extends
   beyond the times in the effort array using the last known effort values. When
   `t_save` is provided, it controls the save frequency with effort values
   interpolated as needed. This allows users to extend simulations without
   specifying dummy effort values for the final time period (#231).
+- The numerical scheme now supports diffusion in the McKendrick-von Foerster
+  equation, allowing individual variability in growth to be modelled. A new
+  `diffusion` slot in `MizerParams` holds the diffusion coefficient (species x
+  size). Use `setDiffusion()` / `diffusion()` / `diffusion<-()` to set and
+  retrieve it.
+- New `getFlux()` function calculates the flux of individuals entering each size
+  class, combining the advective flux from somatic growth and the diffusive flux.
+- `getRequiredRDD()` is now exported. It calculates the recruitment rate needed
+  to maintain a given initial abundance, accounting for both growth and diffusion.
+- `steadySingleSpecies()` now correctly preserves the steady state under
+  `project()`, including when diffusion is non-zero.
+- Growth is now forced to always be non-negative, preventing unphysical shrinkage.
+  No warning is issued when growth stops at or after maturity size.
+- New vignettes: cohort dynamics demonstrating the effect of diffusion in a
+  single-species model; numerical details documenting the finite-volume scheme and
+  its steady-state solution; and a vignette on using FFT for predation kernel
+  calculations.
+- Many functions now have S3 methods so they can be called with either a
+  MizerParams or MizerSim object and users could define their own subclasses
+  and methods to modify mizer behaviour (#330).
+- `getBiomass()` now has a `use_cutoff` argument to restrict the biomass
+  calculation to sizes above the `biomass_cutoff` species parameter.
+- `plotBiomass()` and `plotlyBiomass()` now have a `use_cutoff` argument,
+  passed to `getBiomass()`.
+- `age_mat_vB()` is now exported.
+- `project_n()` is a new exported function that projects the abundance
+  spectrum forward in time, factored out of `project()`.
+
+## Bug fixes
+
+- `getMeanMaxWeight()` now correctly applies the species selector to the
+  denominator.
+- `plotDataFrame()` now correctly applies custom log-scale x breaks.
+- `get_size_range_array()` no longer gives an error when no size brackets are
+  selected.
 
 # mizer 2.5.4
 
@@ -404,7 +439,7 @@ state, so you will need to also call `steady()` after matching the growth rates.
 * Many improvements in the documentation.
 * Many small improvements to code quality and testing.
 * Better social media cards, especially for twitter.
-* mizer can be run on binder, https://mybinder.org/v2/gh/sizespectrum/mizer/HEAD?urlpath=rstudio
+* mizer can be [run on binder](https://mybinder.org/v2/gh/sizespectrum/mizer/HEAD?urlpath=rstudio)
 
 ## Bug fixes
 
@@ -854,7 +889,7 @@ species. The information is set up via a new `gear_params()` data frame. See
   well as `idxFinalT()` to access the values at the final time of a simulation.
 * New function `getCriticalFeedingLevel()` returns the critical feeding level
   for each species at each size.
-* Mizer reexports the `melt()` function from the reshape2 package which allows
+* Mizer re-exports the `melt()` function from the reshape2 package which allows
   users to convert the arrays returned by mizer functions into data frames
   that can be used for example in ggplot2 and plotly.
 * `validSpeciesParams()` checks validity of species parameter data frame and

R/RcppExports.R

@@ -5,3 +5,7 @@ inner_project_loop <- function(no_sp, no_w, n, A, B, S, w_min_idx) {
     .Call('_mizer_inner_project_loop', PACKAGE = 'mizer', no_sp, no_w, n, A, B, S, w_min_idx)
 }
 
+project_n_loop <- function(n, a, b, c, S, w_min_idx) {
+    .Call('_mizer_project_n_loop', PACKAGE = 'mizer', n, a, b, c, S, w_min_idx)
+}
+

R/calibrate.R

@@ -17,7 +17,7 @@
 #' individual species will not match observations yet, with some species
 #' having biomasses that are too high and others too low. So after this
 #' function you may want to use [matchBiomasses()]. This is described in the
-#' blog post at https://bit.ly/2YqXESV.
+#' blog post at \url{https://bit.ly/2YqXESV}.
 #' 
 #' If you have observations of the yearly yield instead of biomasses, you can
 #' use [calibrateYield()] instead of this function.
@@ -84,7 +84,7 @@ calibrateBiomass.MizerParams <- function(params, ...) {
 #' individual species will not match observations yet, with some species
 #' having numbers that are too high and others too low. So after this
 #' function you may want to use [matchNumbers()]. This is described in the
-#' blog post at https://bit.ly/2YqXESV.
+#' blog post at \url{https://bit.ly/2YqXESV}.
 #'
 #' If you have observations of the yearly yield instead of numbers, you can
 #' use [calibrateYield()] instead of this function.

R/getRequiredRDD.R

@@ -0,0 +1,61 @@
+#' Determine reproduction rate needed for initial egg abundance
+#'
+#' @param params A MizerParams object
+#' @return A vector of reproduction rates for all species
+#' @export
+getRequiredRDD <- function(params) {
+    UseMethod("getRequiredRDD")
+}
+
+#' @export
+getRequiredRDD.MizerParams <- function(params) {
+    # Calculate required rdd
+    no_sp <- nrow(params@species_params)
+
+    # Calculate transport coefficients
+    dt <- 1
+    # We pass a dummy recruitment flux of 0 to trigger the boundary condition
+    # corrections for a and b in get_transport_coefs
+    coefs <- get_transport_coefs(params, n = params@initial_n,
+                                 g = getEGrowth(params),
+                                 mu = getMort(params), dt,
+                                 recruitment_flux = numeric(no_sp))
+
+    reproduction <- params@species_params$erepro # vector of correct length
+    names(reproduction) <- params@species_params$species
+
+    for (i in (1:no_sp)) {
+        w_min_idx <- params@w_min_idx[i]
+
+        # Get coefficients for this species at the boundary
+        # The equation for the first node is:
+        # (N_new - N_old)/dt = -(Flux_matrix * N) + R/dw
+        # In steady state N_new = N_old, so:
+        # Flux_matrix * N = R/dw
+        # The rows of coefs correspond to the linear system A*N_{j-1} + B*N_j + C*N_{j+1} = ...
+        # For the first node j=w_min_idx:
+        # A*N_{j-1} + (B-1)/dt * N_j + C/dt * N_{j+1} = R/dw / dt ?
+        # No, let's look at project_n again.
+        # It solves A N_{i-1} + B N_i + C N_{i+1} = N_old + RHS_source
+        # In steady state: A N_{i-1} + B N_i + C N_{i+1} = N_i + R * dt / dw
+        # So R = ( A N_{i-1} + (B-1) N_i + C N_{i+1} ) * dw / dt
+
+        # Extract coefficients
+        a <- coefs$a[i, w_min_idx]
+        b <- coefs$b[i, w_min_idx]
+        c <- coefs$c[i, w_min_idx]
+
+        # Boundary corrections for a and b are now handled in get_transport_coefs
+
+        n_current <- params@initial_n[i, w_min_idx]
+        n_next <- if (w_min_idx < length(params@w)) params@initial_n[i, w_min_idx + 1] else 0
+        n_prev <- if (w_min_idx > 1) params@initial_n[i, w_min_idx - 1] else 0 # Should be irrelevant if A=0 or boundary
+
+        # Calculate R
+        # R = ( A * n_prev + (B - 1) * n_current + C * n_next ) * dw / dt
+
+        total_rate <- a * n_prev + (b - 1) * n_current + c * n_next
+        reproduction[i] <- total_rate * params@dw[w_min_idx] / dt
+    }
+    reproduction
+}

R/get_initial_n.R

@@ -55,21 +55,29 @@ get_initial_n <- function(params, n0_mult = NULL, a = 0.35) {
         p@w_full ^ (-p@resource_params$lambda)
     p@interaction[] <- 0
     income <- getEReproAndGrowth(p) + p@metab
+    mort <- getFMort(p)
+    growth <- getEGrowth(p)
+    N0_vec <- numeric(no_sp)
+
     for (i in seq_len(no_sp)) {
         # At small sizes the income should be A w^n. Determine A
         # Use w_min_idx + 1 in case user has implemented reduced growth
         # for the smallest size class (see e.g. #241)
         iw <- p@w_min_idx[i] + 1
         A <- income[i, iw] / (p@w[iw] ^ p@species_params[[i, "n"]])
 
-        mort <- 0.4 * A * p@w ^ (p@species_params[[i, "n"]] - 1) + getFMort(p)
-        growth <- getEGrowth(p)[i, ]
+        mort[i, ] <- mort[i, ] + 0.4 * A * p@w ^ (p@species_params[[i, "n"]] - 1)
 
-        idxs <- p@w_min_idx[i]:(min(which(c(growth, 0) <= 0)) - 1)
-        idx <- idxs[1:(length(idxs) - 1)]
+        # We start with an arbitrary population at the smallest size class
+        N0_vec[i] <- 1
+    }
+
+    n_exact_matrix <- get_steady_state_n(p, growth, mort, N0_vec)
+
+    for (i in seq_len(no_sp)) {
+        idxs <- p@w_min_idx[i]:(min(which(c(growth[i, ], 0) <= 0)) - 1)
         # Steady state solution of the upwind-difference scheme used in project
-        p@initial_n[i, idxs] <-
-            get_steady_state_n(growth, mort, p@dw, idx)
+        p@initial_n[i, idxs] <- n_exact_matrix[i, idxs]
     }
     p <- matchBiomasses(p)
     return(p@initial_n)

R/helpers.R

@@ -85,20 +85,45 @@ w2l <- function(w, species_params) {
     (w / sp[["a"]])^(1 / sp[["b"]])
 }
 
-#' Helper function to calculate the steady state abundance using the upwind-difference scheme
+#' Calculate steady state abundance
 #'
-#' @param growth A numeric vector of growth rates.
-#' @param mort A numeric vector of mortality rates.
-#' @param dw A numeric vector of the size step.
-#' @param idx A numeric vector of indices.
-#' @param N0 The initial egg density.
-#' @return A numeric vector representing the steady state abundances.
-#' @keywords internal
-get_steady_state_n <- function(growth, mort, dw, idx, N0) {
-    # Steady state solution of the upwind-difference scheme used in project
-    n_exact <- c(1, cumprod(growth[idx] / ((growth + mort * dw)[idx + 1])))
-    if (!missing(N0)) {
-        n_exact <- N0 * n_exact
-    }
-    return(n_exact)
+#' This function calculates the steady state abundance by solving the
+#' transport equation with given growth and mortality rates. It sets up a
+#' tri-diagonal system and solves it.
+#'
+#' @param params A MizerParams object
+#' @param g A matrix of growth rates (species x size)
+#' @param mu A matrix of mortality rates (species x size)
+#' @param N0 A vector with the abundance at the smallest size for each species
+#' @return A matrix with the steady state abundance
+#' @concept helper
+get_steady_state_n <- function(params, g, mu, N0) {
+    no_sp <- nrow(params@species_params)
+    no_w <- length(params@w)
+    n <- matrix(0, nrow = no_sp, ncol = no_w,
+                dimnames = list(params@species_params$species, dimnames(params@initial_n)[[2]]))
+
+    # Use get_transport_coefs to compute the coefficients with dt = 1
+    # and no recruitment flux (since we handle the boundary manually)
+    coefs <- get_transport_coefs(params, n, g, mu, dt = 1,
+                                 recruitment_flux = rep(0, no_sp))
+
+    a <- coefs$a
+    # For steady state, the diagonal term \tilde{B} is B - 1
+    b <- coefs$b - 1
+    c <- coefs$c
+    S <- coefs$S
+
+    # Boundary conditions at the start of the size spectrum:
+    # A_j = 0, B_j = 1, C_j = 0, S_j = N0
+    j_start <- params@w_min_idx
+    idxs_start <- cbind(1:no_sp, j_start)
+    a[idxs_start] <- 0
+    b[idxs_start] <- 1
+    c[idxs_start] <- 0
+    S[idxs_start] <- N0
+
+    n <- project_n_loop(n, a, b, c, S, j_start)
+
+    return(n)
 }

R/manipulate_species.R

@@ -321,6 +321,10 @@ addSpecies.MizerParams <- function(params, species_params, gear_params = data.fr
 #' refer to the selected species. It does not recalculate the steady state for
 #' the remaining species or retune their reproductive efficiency.
 #'
+#' If a gear was targeting only the removed species, then this function will
+#' NOT remove that gear. If you want to also remove that gear then you can do
+#' that by calling [setFishing()].
+#'
 #' @param params A mizer params object for the original system.
 #' @param species The species to be removed. A vector of species names, or a
 #'   numeric vector of species indices, or a logical vector indicating for

R/match.R

@@ -10,7 +10,7 @@
 #' So after using this function you may want to use `steady()` to run the model 
 #' to steady state, after which of course the biomasses will no longer match
 #' exactly. You could then iterate this process. This is described in the
-#' blog post at https://bit.ly/2YqXESV.
+#' blog post at \url{https://bit.ly/2YqXESV}.
 #' 
 #' Before you can use this function you will need to have added a
 #' `biomass_observed` column to your model which gives the observed biomass in
@@ -91,7 +91,7 @@ matchBiomasses.MizerParams <- function(params, species = NULL, ...) {
 #' So after using this function you may want to use `steady()` to run the model
 #' to steady state, after which of course the numbers will no longer match
 #' exactly. You could then iterate this process. This is described in the
-#' blog post at https://bit.ly/2YqXESV.
+#' blog post at \url{https://bit.ly/2YqXESV}.
 #'
 #' Before you can use this function you will need to have added a
 #' `number_observed` column to your model which gives the observed number of
@@ -180,7 +180,7 @@ matchNumbers.MizerParams <- function(params, species = NULL, ...) {
 #' So after using this function you may want to use `steady()` to run the model 
 #' to steady state, after which of course the yields will no longer match
 #' exactly. You could then iterate this process. This is described in the
-#' blog post at https://bit.ly/2YqXESV.
+#' blog post at \url{https://bit.ly/2YqXESV}.
 #' 
 #' Before you can use this function you will need to have added a
 #' `yield_observed` column to your model which gives the observed yields in