From 26159a4ccd851e09a2f2906828772ff5de517756 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 08:28:50 +0100 Subject: [PATCH 01/15] Add getStability() for steady-state dynamic stability analysis Adds two new experimental features: New exported function that analyses the dynamic stability of a mizer steady state by numerically differentiating the one-step-ahead map G(N) = A(N)^{-1} S(N) at the fixed point N*. Uses the same project_n_loop() C++ Thomas solver as the regular dynamics. Returns: - eigenvalues of the linearised map L = dG/dN|_{N*} - spectral_radius (max|lambda_i|); stable iff < 1 - dominant_period (period of dominant eigenvalue in time steps) - hopf_period: period of limit cycle near a Hopf bifurcation (2*pi / |Arg(lambda)| for the complex eigenvalue closest to the unit circle) - n_active: dimension of the Jacobian The include_resource = FALSE default treats the resource as a quasi-static fast variable (valid for semichemostat dynamics), giving a fish-only Jacobian. Setting include_resource = TRUE perturbs both fish and resource cells independently and steps the resource with the full resource_semichemostat() exact solver, yielding the complete coupled Jacobian. The dominant eigenvalues agree closely for typical mizer models where the resource dynamics are fast. When stability = TRUE, calls getStability() after finding the steady state and attaches the result as attr(params, 'stability'). Seven new tests in test-steadyNewton.R covering list structure, spectral radius < 1 for the NS model, eigenvalue/spectral_radius consistency, stability attribute attachment, hopf_period validity, and agreement between full and reduced analyses. - pkgdown/_pkgdown.yml: getStability added to Steady state tuning - NEWS.md: entry for both new features --- NEWS.md | 172 +++---- R/steadyNewton.R | 726 +++++++++++++++++++++++++++++ man/getStability.Rd | 107 +++++ man/steadyNewton.Rd | 136 ++++++ pkgdown/_pkgdown.yml | 2 + tests/testthat/test-steadyNewton.R | 322 +++++++++++++ 6 files changed, 1385 insertions(+), 80 deletions(-) create mode 100644 R/steadyNewton.R create mode 100644 man/getStability.Rd create mode 100644 man/steadyNewton.Rd create mode 100644 tests/testthat/test-steadyNewton.R diff --git a/NEWS.md b/NEWS.md index 87fd7f139..4ac276884 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,98 +1,110 @@ -# mizer 3.1.0.9000 (development version) - -- Assigning to `resource_params()` no longer balances the resource. It now only - rebuilds the size-dependent rate (`rr_pp`) and capacity (`cc_pp`) arrays from - the scalars, leaving any manually set (frozen) arrays untouched, exactly as - `species_params<-` feeds the species rates. As a result changing the rate-side - scalars `r_pp` or `n` now takes effect (previously the value was silently - discarded by balancing), and successive scalar edits accumulate instead of - overwriting each other. Balancing the resource to preserve the steady state is - now solely a feature of `setResource()`. The `resource_rate<-`, - `resource_capacity<-`, `resource_level<-` and `resource_dynamics<-` setters - gained a `balance` argument (default unchanged) so it can be switched off, e.g. - `resource_capacity(params, balance = FALSE) <- my_capacity`. `setResource()` - also no longer silently overwrites a manually set rate or capacity array when - balancing: the frozen array wins and a warning is issued. +# mizer 3.2.0 + +This release corrects and refines the resource-setting interface, makes the +extension framework composable regardless of load order, adds new analysis +functions and cheatsheets, and includes a range of smaller improvements and bug +fixes. + +For an overview see the release announcement on the mizer blog. + +## Resource setting + +- Assigning to `resource_params()` now rebuilds the size-dependent resource + capacity and resource rate arrays from the resource parameters, but leaving + any manually set (frozen) arrays untouched, mirroring the way + `species_params<-` feeds the species rate arrays. + +- The `resource_rate<-`, `resource_capacity<-`, `resource_level<-` and + `resource_dynamics<-` setters gained a `balance` argument so balancing can be + switched off, e.g. `resource_capacity(params, balance = FALSE) <- my_capacity`. + +- `setResource()` no longer silently overwrites a manually set (frozen) rate or + capacity array when balancing: the frozen array wins and a warning is issued. + +These changes and how to adapt existing code are described in the new +`vignette("upgrading")` ("Upgrading mizer"). + +## Extension framework - An installed extension package is now recognised as a dispatching extension from the S3 methods it registers for its marker class (e.g. `getEncounter.mizerMR`), rather than only from a statically defined S4 marker - class. This lets extension packages omit the static - `setClass("mizerFoo", contains = "MizerParams")` and instead let mizer create - the marker class dynamically when the package is loaded, inserting it at the - correct place in the S4 hierarchy relative to any other extensions loaded in - the session. As a result, two independently developed extension packages - (for example mizerReef and mizerMR) can now be chained in either load order, - which a static sibling-of-`MizerParams` class prevented. + class. Extension packages can therefore omit the static + `setClass("mizerFoo", contains = "MizerParams")` and let mizer create the + marker class dynamically on load, inserting it at the correct place in the S4 + hierarchy. As a result two independently developed extension packages (for + example mizerReef and mizerMR) can now be chained in either load order, which a + static sibling-of-`MizerParams` class prevented. - `recordExtension()` now prepends a genuinely new extension to the front of the object's `@extensions` chain, keeping it ordered outermost-first to match - `registerExtension()`. Existing entries stay in place. - -- `print()` on `ArraySpeciesBySize`, `ArrayTimeBySpecies`, `ArrayResourceBySize`, - `ArrayTimeByResourceBySize` and `ArrayTimeBySpeciesBySize` objects (as - returned by `getEncounter()`, `getBiomass()`, `getResourceMort()`, `NResource()`, - `getFMort()` and similar functions) now prints the array's actual values - instead of a per-species min/mean/max summary, truncating large arrays to - fit the console: species are shown as a leading subset, sizes as an evenly - log-spaced sample across the full size range, and time series as a - representative sample of time steps that always includes the first and - last, with a note reporting how much was omitted. A three-dimensional - `ArrayTimeBySpeciesBySize` object is previewed via its final time slice, - matching `plot()`'s existing default for that class. Use - `as.data.frame()` for full, untruncated access to the data. - -- `plotYieldGear()` now supports `log_x`, `log_y`, and `log` arguments, - aligning its arguments with `plotYield()`. - -- Added three new cheatsheets: "Model Setup and Calibration" (building a model, - finding the steady state, calibrating to observed biomass/yield/growth, and - projecting), "Fishing" (gears, selectivity functions, catchability, and - effort), and "Changing Model Parameters" (the distinction between - `given_species_params()`, `calculated_species_params()` and `species_params()`, - when changing a species parameter updates a size-dependent rate versus freezing - it, and how `gear_params()` and the resource setters differ). - -- The analysis-and-plotting cheatsheet now covers the newer plotting functions - (`plotCDF()`, `plotSpectra2()`, `plotSpectraRelative()`, `plotCDF2()`, - `plot2()`, `plotRelative()`, `animate()`) and the `ArrayResourceBySize` - class, and corrects the interactive-plot advice for array objects to use - `plotHover()` (they have no `ggplotly()` method). + `registerExtension()`. + +## New functions + +- New `adjustSizeGrid()` (an S3 generic) adjusts the size grid of a + `MizerParams` object to a new minimum and/or maximum size. It can both expand + and truncate the grid, warning if non-negligible abundance is discarded. + +- New experimental `steadyNewton()` finds a steady state by solving the + steady-state equation directly with a Newton-type root finder (via the + `nleqslv` package) instead of running the dynamics to convergence. Unlike + `steady()` it converges even when the steady state is dynamically unstable, and + it discovers the support of the steady state automatically. Currently supports + the default semichemostat resource dynamics. + +- New experimental `getStability()` analyses the dynamic stability of a mizer + steady state by computing the eigenvalues of the linearised one-step-ahead map + at the fixed point. It reports whether the steady state is stable or unstable, + the spectral radius, and — when the system approaches a Hopf bifurcation — the + period of the emergent limit cycle. By default the resource is treated as a + quasi-static fast variable (valid for semichemostat dynamics); setting + `include_resource = TRUE` gives the full coupled (fish + resource) Jacobian, + useful for verifying that the quasi-static approximation makes little difference. + `steadyNewton()` gains a `stability = TRUE` argument that calls `getStability()` + automatically and attaches the result as the attribute `"stability"` on the + returned `MizerParams` object. + +- `project()` gains a `callback` argument for a user-defined function to be + called at each saved time step. -- When `sel_func` is set on a `gear_params` object, any argument columns - required by that selectivity function (other than `w`, `species_params`, and - `...`) are now automatically added as `NA` columns. This means, for example, - that setting `gp$sel_func <- "sigmoid_length"` immediately creates the `l25` - and `l50` columns, ready to be filled in (#431). +## Other improvements -- The `species_params` data frame is now an S3 subclass of `data.frame` - (`class = c("species_params", "data.frame")`). It supports class-preserving - subsetting and subassignment S3 methods, making it safer to use and paving - the way for future auto-recalculations. +- `print()` on the array classes returned by the rate getters (e.g. + `getEncounter()`, `getBiomass()`, `getFMort()`, `NResource()`) now prints the + array's actual values, truncated to fit the console, instead of a per-species + min/mean/max summary. Use `as.data.frame()` for full, untruncated access. -- `resource_params<-` now rebuilds the resource size-spectrum arrays (`cc_pp`, - `rr_pp`) via `setResource()`, making resource parameter assignment - consistent with `species_params<-` and `gear_params<-` (#439). +- Columns accessed via `$` on a `species_params` or `gear_params` object now + return vectors named by species (or by "species, gear" for `gear_params`), + making entries easier to identify. The `species` column is left unnamed. -- Added a `callback` parameter to `project()` to allow user-defined functions - to be called at each saved time step. +- The `species_params` data frame is now an S3 subclass of `data.frame`, with + class-preserving subsetting and subassignment methods, paving the way for + future auto-recalculations. -- Columns accessed via `$` on a `species_params` or `gear_params` object now - return named vectors, where the names are the species names (or "species, - gear" row names for `gear_params`). For example, - `species_params(params)$w_mat` now returns a named vector making it easier - to identify entries. The `species` vector is left unnamed. +- When `sel_func` is set on a `gear_params` object, the argument columns that + selectivity function requires are now added automatically as `NA` columns. For + example, setting `gp$sel_func <- "sigmoid_length"` immediately creates the + `l25` and `l50` columns, ready to be filled in (#431). -- New `adjustSizeGrid()` function (an S3 generic) adjusts the size grid of - a `MizerParams` object to a new minimum and/or maximum size. It can both - expand and truncate (shrink) the grid, warning if non-negligible abundance - (species or resource) is discarded. +- `plotYieldGear()` gains `log_x`, `log_y` and `log` arguments, aligning it with + `plotYield()`. - The upper boundary condition of the size-spectrum solvers now holds the - abundance at zero above each species' maximum size `w_max`. Without diffusion - this is automatic and results are unchanged, but with predation diffusion - switched on it stops a small amount of density leaking to sizes above `w_max`. - See the "Numerical Details" vignette. + abundance at zero above each species' maximum size `w_max`. Results are + unchanged without diffusion, but with predation diffusion switched on this + stops a small amount of density leaking to sizes above `w_max`. See the + "Numerical Details" vignette. + +## Documentation + +- Added three new cheatsheets: "Model Setup and Calibration", "Fishing", and + "Changing Model Parameters". + +- The analysis-and-plotting cheatsheet now covers the newer plotting functions + (`plotCDF()`, `plotSpectra2()`, `plotSpectraRelative()`, `plotCDF2()`, + `plot2()`, `plotRelative()`, `animate()`) and the `ArrayResourceBySize` class. # mizer 3.1.0 diff --git a/R/steadyNewton.R b/R/steadyNewton.R new file mode 100644 index 000000000..86ec3652f --- /dev/null +++ b/R/steadyNewton.R @@ -0,0 +1,726 @@ +# Direct steady-state solver ------------------------------------------------- +# +# `steady()` finds a steady state by running the dynamics until they stop +# changing. That only works for steady states that are dynamically stable: at an +# unstable steady state the time-stepping diverges away from the fixed point. +# +# `steadyNewton()` instead solves the discrete steady-state equation +# `F(N) = 0` directly with a Newton-type root finder (via the `nleqslv` +# package). Because it solves the algebraic equation rather than following the +# time evolution, it converges to the steady state irrespective of its dynamic +# stability. See the "Steady-State Solution" section of the +# `vignette("numerical_details")` for the equation that is being solved. + +#' Find a steady state by directly solving the steady-state equation +#' +#' `r lifecycle::badge("experimental")` +#' This is an alternative to [steady()] that finds the steady state by solving +#' the steady-state equation `F(N) = 0` with a Newton-type root finder instead +#' of running the dynamics to convergence. The advantage is that it converges +#' to the steady state even when that steady state is dynamically **unstable**, +#' a case in which [steady()] fails because the time-stepping diverges away from +#' the fixed point. +#' +#' By default, or when `reproduction = "fixed"`, the function holds the +#' reproduction rate (RDD) constant while solving for the consumer spectra, +#' substitutes the analytic steady state of the resource, and keeps any other +#' components constant. After the spectra have been found it restores +#' density-dependent Beverton-Holt reproduction with [setBevertonHolt()], +#' honouring the `preserve` argument exactly as [steady()] does. +#' If `reproduction = "dynamic"`, the reproduction dynamics are run dynamically, +#' meaning the reproduction rate varies during the solve and the reproduction +#' parameters are not adjusted. +#' +#' The consumer densities are solved for in log space, which both keeps them +#' positive and conditions the otherwise badly-scaled system. The unknowns are +#' the size classes in the full potential grid support (running from the egg +#' size up to the grid truncation limit `support_top_idx()`), regardless of +#' whether they carry non-zero density in the supplied initial spectra. A +#' smoothed log-abundance penalty floor is applied to all size classes, which +#' automatically handles zero-density and tail classes smoothly and prevents +#' singular Jacobians. After convergence, size classes that remain at or near +#' the floor are set to zero. This allows the solver to automatically discover +#' the support of the steady state. The nonlinear system is solved with a +#' globalised Newton iteration from the `nleqslv` package, starting from the +#' current `initial_n`. Newton's method converges from any starting point in the +#' *root's* basin of attraction (which is unrelated to the dynamic stability of +#' the steady state), so a reasonable initial guess should still be supplied in +#' `initialN(params)` — for example the spectra from a nearby stable +#' parameterisation, or the (diverging) output of [steady()]. +#' +#' The solver respects the active transport scheme: if the experimental +#' second-order scheme is enabled (see [second_order_w()]) it solves the +#' steady-state equation of that scheme. With the van Leer reconstruction the +#' residual is only Lipschitz, so the iteration converges to a fixed point of the +#' dynamics but not to machine precision. The unlimited `"centred"` +#' reconstruction admits an undamped odd-even mode at a steady state with no +#' physical diffusion, giving an ill-conditioned steady-state Jacobian for which +#' the solver is not expected to converge. +#' +#' Only the default semichemostat resource dynamics +#' (`resource_dynamics = "resource_semichemostat"`) are currently supported, +#' because the solver substitutes the analytic resource steady state. For other +#' resource dynamics the function stops with an error. +#' +#' @param params A \linkS4class{MizerParams} object. Its `initial_n` is used as +#' the starting guess for the iteration. +#' @param effort The fishing effort. By default the initial effort stored in +#' `params`. +#' @param preserve `r lifecycle::badge("experimental")` +#' Specifies whether the `reproduction_level` should be preserved (default) +#' or the maximum reproduction rate `R_max` or the reproductive efficiency +#' `erepro`. See [setBevertonHolt()] for an explanation of the +#' `reproduction_level`. This argument is ignored when `reproduction = "dynamic"`. +#' @param reproduction `r lifecycle::badge("experimental")` +#' If `"fixed"`, the reproduction rate (RDD) is held constant at the initial +#' value. If `"dynamic"`, the reproduction dynamics are run dynamically and the +#' reproduction parameters are not adjusted. Default is `"fixed"`. +#' @param extinction_floor `r lifecycle::badge("experimental")` +#' The relative abundance floor below which a species is considered extinct. +#' Only used when `reproduction = "dynamic"`. Default is 1e-6. +#' @param verbose If `TRUE` then the solver iterations will be traced and +#' printed to the console. Default is `FALSE`. +#' @param tol Convergence tolerance passed to [nleqslv::nleqslv()] (both the +#' function-value tolerance `ftol` and the step tolerance `xtol`). +#' @param maxit Maximum number of iterations for [nleqslv::nleqslv()]. +#' @param method The [nleqslv::nleqslv()] method, either `"Newton"` (with a +#' numerical Jacobian calculated at each iteration) or `"Broyden"` (which +#' calculates the full Jacobian only once and then only updates it on each +#' iteration. '"Broyden"' is the default. +#' @param global The globalisation strategy passed to [nleqslv::nleqslv()]. +#' The default `"dbldog"` (double dogleg) is a robust trust-region method. +#' @param stability `r lifecycle::badge("experimental")` +#' If `TRUE`, [getStability()] is called after convergence and its result is +#' stored as the attribute `"stability"` on the returned +#' \linkS4class{MizerParams} object. Default is `FALSE`. +#' @param ... Unused. +#' @return A \linkS4class{MizerParams} object with the initial state set to the +#' steady state. If `stability = TRUE`, the object carries the attribute +#' `"stability"` containing the list returned by [getStability()]. +#' @seealso [steady()], [steadySingleSpecies()], [getStability()] +#' @export +#' @examples +#' \donttest{ +#' params <- steadyNewton(NS_params) +#' plotSpectra(params) +#' params <- steadyNewton(NS_params, stability = TRUE) +#' attr(params, "stability")$stable +#' } +steadyNewton <- function(params, ...) { + UseMethod("steadyNewton") +} + +#' @rdname steadyNewton +#' @export +steadyNewton.MizerParams <- function(params, + effort = params@initial_effort, + preserve = c("reproduction_level", + "erepro", "R_max"), + reproduction = c("fixed", "dynamic"), + extinction_floor = 1e-6, + verbose = FALSE, + tol = 1e-6, maxit = 200, + method = c("Broyden", "Newton"), + global = "dbldog", + stability = FALSE, ...) { + reproduction <- match.arg(reproduction) + if (reproduction == "fixed") { + preserve <- match.arg(preserve) + } + method <- match.arg(method) + if (!requireNamespace("nleqslv", quietly = TRUE)) { + stop("steadyNewton() requires the 'nleqslv' package. ", + "Install it with install.packages('nleqslv').") + } + params <- validParams(params) + if (params@resource_dynamics != "resource_semichemostat") { + stop("steadyNewton() currently only supports the default semichemostat ", + "resource dynamics (resource_dynamics = 'resource_semichemostat'). ", + "Use steady() for other resource dynamics.") + } + effort <- validEffortVector(effort, params = params) + params@initial_effort <- effort + + if (params@rates_funcs$RDD == "BevertonHoltRDD" && reproduction == "fixed") { + old_reproduction_level <- getReproductionLevel(params) + old_R_max <- params@species_params$R_max + old_erepro <- params@species_params$erepro + } + + if (reproduction == "dynamic") { + rdd_const <- NULL + } else { + rdd_const <- getRDD(params) + } + n_other <- params@initial_n_other + + active <- steady_active_set(params) + residual_fn <- steady_state_residual(params, rdd_const, n_other, effort, + active, extinction_floor = extinction_floor) + + # Save initial state for relative floor checks + N0_initial <- params@initial_n + + # The log-space solve needs a strictly positive, well-scaled start. The + # second-order schemes can leave isolated zeros inside the support (a + # negativity-floor artefact), which would make the 1/N-scaled residual + # overflow. Fill those by log-interpolation from the nonzero neighbours. + N0 <- positive_initial_guess(params@initial_n, params@w_min_idx, + active$w_top) + x0 <- log(N0[active$mask]) + sol <- nleqslv::nleqslv(x0, residual_fn, method = method, + global = global, + control = list(maxit = maxit, ftol = tol, + xtol = tol, + trace = if (verbose) 1 else 0)) + # termcd 1 (function convergence) and 2 (x convergence) are successes. + if (sol$termcd > 2) { + warning("steadyNewton() did not converge (nleqslv termination code ", + sol$termcd, ": ", sol$message, + "). Returning the best iterate found.", call. = FALSE) + } + + N <- active$unpack(sol$x) + n_pp <- resource_steady_semichemostat(params, N, n_other) + + # Zero out any densities that ended up at or near the penalty floor + for (i in seq_len(nrow(N))) { + floor_val <- max(N[i, ]) * 1.1 * 1e-15 + N[i, N[i, ] < floor_val] <- 0 + } + + params@initial_n[] <- N + params@initial_n_pp[] <- n_pp + + # Check for extinctions if using the relative floor + if (reproduction == "dynamic" && !is.null(extinction_floor) && extinction_floor > 0) { + extinct_threshold <- extinction_floor * 1.01 + is_extinct <- rep(FALSE, nrow(N)) + names(is_extinct) <- rownames(N) + for (i in seq_len(nrow(N))) { + lo <- params@w_min_idx[i] + pos <- lo:active$w_top[i] + pos <- pos[N0_initial[i, pos] > 0] + if (length(pos) && any(N[i, pos] / N0_initial[i, pos] <= extinct_threshold)) { + is_extinct[i] <- TRUE + } + } + if (any(is_extinct)) { + warning("The following species went extinct and were pegged to their abundance floor: ", + paste(names(is_extinct)[is_extinct], collapse = ", ")) + } + } + + # Restore density-dependent reproduction, just as steady() does. + if (params@rates_funcs$RDD == "BevertonHoltRDD" && reproduction == "fixed") { + if (preserve == "reproduction_level") { + params <- setBevertonHolt(params, + reproduction_level = old_reproduction_level) + } else if (preserve == "R_max") { + params <- setBevertonHolt(params, R_max = old_R_max) + } else if (preserve == "erepro") { + params <- setBevertonHolt(params, erepro = old_erepro) + } + } + + params@time_modified <- lubridate::now() + + if (stability) { + attr(params, "stability") <- getStability(params, + reproduction = reproduction, + effort = effort, + extinction_floor = extinction_floor) + } + + params +} + +#' The set of size classes solved for by steadyNewton() +#' +#' Builds the logical mask of the (species x size) density matrix that the +#' direct solver treats as unknowns. For each species the unknowns run from the +#' egg size `w_min_idx` up to the highest class that actually carries density in +#' the current `initial_n`, capped at the grid truncation `support_top_idx()`. +#' +#' The support is read off the abundances rather than from `w_max` on purpose. +#' Users routinely set `w_max` far larger than necessary (so the grid need not +#' change when a parameter change produces larger fish), which would otherwise +#' put a whole band of structurally-zero classes — `log(0)` unknowns that make +#' the Jacobian singular — into the system. The growth rate alone cannot locate +#' the top either: the main reason fish grow beyond `w_repro_max` is diffusion, +#' and the diffusion rate only grows with `w`, never vanishing. The abundance is +#' therefore the only reliable indicator of where the (possibly diffusion-fed) +#' tail has died away. The distinction is clean: without diffusion the classes +#' above where growth stops receive no inflow and are *exactly* zero, while with +#' diffusion the tail is genuinely non-zero up to where it decays away (or to the +#' grid truncation). +#' +#' Isolated interior zeros (negativity-floor artefacts of the second-order +#' schemes) below the top are kept in the mask and repaired by +#' `positive_initial_guess()`; only a trailing run of zeros is dropped. +#' +#' @param params A \linkS4class{MizerParams} object. +#' @return A list with the logical matrix `mask` and a function `unpack(x)` that +#' maps a vector of log-densities (in `mask` order) to the full density +#' matrix. +#' @noRd +steady_active_set <- function(params) { + no_sp <- nrow(params@species_params) + no_w <- length(params@w) + + # The grid truncation caps the support: the dynamics hold the density at zero + # above this, so it can never be an unknown (see support_top_idx()). + grid_top <- support_top_idx(params) + n <- params@initial_n + w_top <- integer(no_sp) + mask <- matrix(FALSE, nrow = no_sp, ncol = no_w) + for (i in seq_len(no_sp)) { + lo <- params@w_min_idx[i] + # Always solve up to the grid truncation limit + w_top[i] <- grid_top[i] + mask[i, lo:w_top[i]] <- TRUE + } + + dn <- dimnames(params@initial_n) + unpack <- function(x) { + N <- matrix(0, nrow = no_sp, ncol = no_w, dimnames = dn) + N[mask] <- exp(x) + N + } + list(mask = mask, w_top = w_top, unpack = unpack) +} + +#' Strictly positive starting guess for the log-space solve +#' +#' The Newton solve uses log-densities as unknowns and scales the residual by +#' `1/N`, so it needs a starting guess that is strictly positive and reasonably +#' scaled on every active size class. The second-order schemes can leave isolated +#' zeros inside the support (a negativity-floor artefact at a reconstructed +#' over/undershoot); a plain `log()` of such a guess is `-Inf`, and flooring it +#' to a tiny constant instead makes the `1/N`-scaled residual overflow. We repair +#' the guess by interpolating `log(N)` linearly across the gaps from the nonzero +#' neighbours (geometric interpolation on the logarithmic size grid), which is +#' both strictly positive and well scaled. Edges with no nonzero neighbour on one +#' side are extrapolated flat. +#' +#' @param N The current density matrix (species x size). +#' @param w_min_idx Per-species egg-size index. +#' @param w_top Per-species support top from `steady_active_set()`. +#' @return A density matrix that is strictly positive on the active range. +#' @noRd +positive_initial_guess <- function(N, w_min_idx, w_top) { + for (i in seq_len(nrow(N))) { + rng <- w_min_idx[i]:w_top[i] + v <- N[i, rng] + pos <- v > 0 + if (!any(pos)) next + floor_val <- max(v) * 1e-20 + if (!all(pos)) { + idx <- seq_along(v) + v[!pos] <- exp(stats::approx(idx[pos], log(v[pos]), + xout = idx[!pos], rule = 2)$y) + } + v[is.na(v) | v <= 0] <- floor_val + N[i, rng] <- v + } + N +} + +#' Analytic semichemostat resource steady state +#' +#' Returns the resource number density that is in equilibrium with the consumer +#' densities `N`. For the semichemostat resource the equilibrium is +#' `n_pp* = rr_pp * cc_pp / (rr_pp + mu_R(N))`, where the resource predation +#' mortality `mu_R` depends only on the consumer densities (not on the resource +#' density itself), so the substitution is exact. This is the `n_steady` of +#' [resource_semichemostat()]. +#' +#' @param params A \linkS4class{MizerParams} object. +#' @param N Consumer densities (species x size). +#' @param n_other Abundances of other components. +#' @return The steady-state resource number density vector. +#' @noRd +resource_steady_semichemostat <- function(params, N, n_other) { + # Resource mortality is predation by consumers and does not depend on the + # resource density, so the value of n_pp passed here is irrelevant. + mu_R <- as.numeric(getResourceMort(params, n = N, + n_pp = params@initial_n_pp, + n_other = n_other, t = 0)) + mur <- params@rr_pp + mu_R + n_pp <- params@rr_pp * params@cc_pp / mur + # Where both rate and mortality vanish the steady state is undetermined; + # keep the current value, exactly as resource_semichemostat() does. + sel <- !is.finite(n_pp) + n_pp[sel] <- params@initial_n_pp[sel] + n_pp +} + +#' Residual of the discrete steady-state equation +#' +#' Returns a closure `f(x)` suitable for [nleqslv::nleqslv()]. The argument `x` +#' is the vector of log-densities of the active size classes (see +#' `steady_active_set()`). The closure rebuilds the full density matrix, +#' substitutes the analytic resource steady state, evaluates the growth, +#' mortality and diffusion rates, assembles the steady-state tridiagonal +#' coefficients with `get_transport_coefs()` at `dt = 1`, and returns the +#' steady-state residual +#' \deqn{F_j = a_j N_{j-1} + b_j N_j + c_j N_{j+1} - S_j} +#' (since `S = N + recruitment source`, the `+N` cancels the backward-Euler +#' `N^t` term, leaving exactly the steady-state equation +#' \eqn{\tilde A N_{j-1} + \tilde B N_j + \tilde C N_{j+1} - \tilde S = 0}). +#' The residual is divided by `N`, turning it into a per-capita rate of change +#' that is dimensionless and O(1) across the many orders of magnitude spanned by +#' the densities — the natural scaling to pair with the log-space unknowns. +#' +#' The flux-limiter weight is recomputed fresh on every evaluation, so the van +#' Leer / centred nonlinearity is part of the root and the solution matches what +#' [project()] converges to on the same `params`. +#' +#' @param params A \linkS4class{MizerParams} object. +#' @param rdd_const Per-species reproduction rate held constant during the solve. +#' @param n_other Abundances of other components (held constant). +#' @param effort The fishing effort vector. +#' @param active The active-set list from `steady_active_set()`. +#' @return A function of the packed log-density vector returning the packed +#' scaled residual. +#' @noRd +steady_state_residual <- function(params, rdd_const, n_other, effort, active, + extinction_floor = 1e-6) { + no_w <- length(params@w) + mask <- active$mask + flux_limiter <- flux_limiter_scheme(params) + rates_fns <- projectRateFunctions(params) + x0_initial <- params@initial_n + + # A floor is always active to handle zero-abundance tail classes smoothly + # and prevent singular Jacobians. + x0 <- log(x0_initial[mask]) + support_floor <- matrix(0, nrow = nrow(x0_initial), ncol = no_w) + for (i in seq_len(nrow(x0_initial))) { + support_floor[i, ] <- log(max(x0_initial[i, ])) + log(1e-15) + } + x_floor <- support_floor[mask] + + if (is.null(rdd_const) && !is.null(extinction_floor) && extinction_floor > 0) { + ext_floor <- x0 + log(extinction_floor) + x_floor <- pmax(x_floor, ext_floor) + } + + function(x) { + N <- active$unpack(x) + n_pp <- resource_steady_semichemostat(params, N, n_other) + + r <- mizer_rates_subset(params, n = N, n_pp = n_pp, n_other = n_other, + t = 0, effort = effort, rates_fns = rates_fns, + targets = c("EGrowth", "Mort", "Diffusion")) + + if (is.null(rdd_const)) { + if (usesExtensionDispatch(params)) { + rdi <- projectRDI(params, n = N, n_pp = n_pp, n_other = n_other, t = 0, + e_repro = r$e_repro, e_growth = r$e_growth, mort = r$mort, + diffusion = r$diffusion) + rdd <- projectRDD(params, rdi = rdi, species_params = params@species_params, + t = 0) + } else { + f_rdi <- get(params@rates_funcs$RDI) + rdi <- f_rdi(params, n = N, n_pp = n_pp, n_other = n_other, t = 0, + e_repro = r$e_repro, e_growth = r$e_growth, mort = r$mort, + diffusion = r$diffusion) + f_rdd <- get(params@rates_funcs$RDD) + rdd <- f_rdd(rdi = rdi, species_params = params@species_params, + params = params, t = 0) + } + } else { + rdd <- rdd_const + } + + coefs <- get_transport_coefs(params, n = N, g = r$e_growth, + mu = r$mort, dt = 1, + recruitment_flux = rdd, + d = r$diffusion, + flux_limiter = flux_limiter) + + Nm <- cbind(0, N[, -no_w, drop = FALSE]) # N_{j-1} + Np <- cbind(N[, -1, drop = FALSE], 0) # N_{j+1} + res <- coefs$a * Nm + coefs$b * N + coefs$c * Np - coefs$S + + # Scale to a per-capita rate of change (N > 0 on the active set). + if (is.null(rdd_const)) { + r_ss <- (res / x0_initial)[mask] + } else { + r_ss <- (res / N)[mask] + } + + if (!is.null(x_floor)) { + y <- x - x_floor + delta <- 0.001 + penalty <- 0.5 * (y - sqrt(y^2 + delta^2)) + r_ss <- r_ss + penalty + } + r_ss + } +} + +# Stability analysis ---------------------------------------------------------- + +#' Analyse the dynamic stability of a mizer steady state +#' +#' `r lifecycle::badge("experimental")` +#' Computes the eigenvalues of the linearised one-step-ahead map at the steady +#' state stored in `params@initial_n`. These eigenvalues determine whether +#' the steady state is dynamically stable and, when a Hopf bifurcation is +#' approached, the period of the emergent limit cycle. +#' +#' ## Mathematical background +#' +#' The mizer time step applies a backward-Euler transport solve for the fish: +#' \deqn{A(N^t, n_{pp}^t)\,N^{t+1} = S(N^t, n_{pp}^t),} +#' and an exact semi-chemostat update for the resource: +#' \deqn{n_{pp}^{t+1} = n_{pp}^* + (n_{pp}^t - n_{pp}^*)\,e^{-\mu^t\,dt},} +#' where \eqn{n_{pp}^* = r_{pp}\,c_{pp}/\mu^t} is the resource steady state +#' conditioned on the mortality \eqn{\mu^t} due to consumers at time \eqn{t}. +#' +#' The stability is determined by the Jacobian of the full one-step-ahead map +#' \eqn{G : (N, n_{pp}) \mapsto (N^{t+1}, n_{pp}^{t+1})} at the fixed point. +#' +#' When `include_resource = FALSE` (the default), the resource is treated as a +#' fast variable that adjusts *instantaneously* to the consumer abundance: for +#' each perturbed \eqn{N}, \eqn{n_{pp}} is set to its quasi-static equilibrium +#' \eqn{n_{pp}^*(N)}. The resulting reduced Jacobian \eqn{L_{\text{red}}} has +#' dimension equal to the number of active fish cells. This is equivalent to +#' projecting the full dynamics onto the slow manifold \eqn{n_{pp} = n_{pp}^*(N)}. +#' +#' When `include_resource = TRUE`, both fish and resource cells are perturbed +#' independently and the full coupled Jacobian \eqn{L_{\text{full}}} is +#' returned. Its eigenvalues include both the slow fish modes and a cluster of +#' fast resource-relaxation modes (with modulus \eqn{e^{-\mu\,dt} \ll 1}). +#' Comparing the dominant eigenvalues of the two analyses shows how much the +#' quasi-static approximation affects the stability conclusion. +#' +#' The steady state is **stable** when all eigenvalues satisfy +#' \eqn{|\lambda_i| < 1} and **unstable** when at least one exceeds 1. +#' +#' A **Hopf bifurcation** occurs when a complex-conjugate pair of eigenvalues +#' crosses the unit circle, giving a limit-cycle period +#' \deqn{T = \frac{2\pi}{|\arg(\lambda)|} \text{ time steps.}} +#' For the default mizer time step of one year this is in years. +#' +#' Both branches use the same `project_n_loop()` C++ Thomas solver as the +#' regular dynamics and centred finite differences with relative step \eqn{h}. +#' +#' @param params A \linkS4class{MizerParams} object whose `initial_n` holds the +#' steady state to analyse. Typically the output of [steadyNewton()]. +#' @param reproduction Whether the reproduction rate is held fixed (`"fixed"`, +#' default) or run dynamically (`"dynamic"`) during the one-step evaluation. +#' Must match the choice used when the steady state was computed. +#' @param effort The fishing effort to use. By default the initial effort +#' stored in `params`. +#' @param include_resource If `FALSE` (default) the resource is treated as a +#' quasi-static fast variable: for each perturbed fish abundance the resource +#' is set to its analytic steady-state value conditioned on that fish +#' abundance (valid only for semichemostat resource dynamics). If `TRUE`, +#' both fish and resource cells are perturbed independently and the resource +#' is evolved with the full resource dynamics function stored in +#' `params@resource_dynamics`, giving the complete coupled Jacobian. +#' @param extinction_floor Relative abundance floor for the dynamic reproduction +#' case. Default is `1e-6`. +#' @param h Relative step size for centred finite differences. Default `1e-4`. +#' @return A named list with the following components: +#' \describe{ +#' \item{`eigenvalues`}{Complex vector of all eigenvalues of the +#' linearised Jacobian, sorted by decreasing modulus.} +#' \item{`spectral_radius`}{The largest modulus \eqn{\max_i|\lambda_i|}.} +#' \item{`stable`}{Logical: `TRUE` when `spectral_radius < 1`.} +#' \item{`dominant_period`}{The period (in time steps) of the dominant +#' eigenvalue: `2*pi / abs(Arg(lambda_1))`. `Inf` for a real positive +#' dominant eigenvalue (monotone dynamics); `NA` for a period-2 flip.} +#' \item{`hopf_period`}{Period (in time steps) of the complex eigenvalue +#' whose modulus is closest to 1; `NULL` when no complex eigenvalue +#' exists. This is the expected limit-cycle period near a Hopf +#' bifurcation.} +#' \item{`n_active`}{Dimension of the Jacobian: number of active fish cells +#' when `include_resource = FALSE`, or fish cells plus all resource cells +#' when `include_resource = TRUE`.} +#' } +#' @seealso [steadyNewton()] +#' @export +getStability <- function(params, + reproduction = c("fixed", "dynamic"), + effort = params@initial_effort, + include_resource = FALSE, + extinction_floor = 1e-6, + h = 1e-4) { + reproduction <- match.arg(reproduction) + params <- validParams(params) + effort <- validEffortVector(effort, params = params) + params@initial_effort <- effort + + if (reproduction == "dynamic") { + rdd_const <- NULL + } else { + rdd_const <- getRDD(params) + } + n_other <- params@initial_n_other + active <- steady_active_set(params) + flux_limiter <- flux_limiter_scheme(params) + rates_fns <- projectRateFunctions(params) + active_idx <- which(active$mask) + + # ------------------------------------------------------------------------- + # Shared helper: one step of fish dynamics given N and n_pp. + # The n_pp here is whatever we choose to pass (quasi-static or actual). + # ------------------------------------------------------------------------- + fish_step <- function(N_in, npp_in) { + r <- mizer_rates_subset(params, n = N_in, n_pp = npp_in, + n_other = n_other, t = 0, effort = effort, + rates_fns = rates_fns, + targets = c("EGrowth", "Mort", "Diffusion", + "ResourceMort")) + + if (is.null(rdd_const)) { + if (usesExtensionDispatch(params)) { + rdi <- projectRDI(params, n = N_in, n_pp = npp_in, + n_other = n_other, t = 0, + e_repro = r$e_repro, e_growth = r$e_growth, + mort = r$mort, diffusion = r$diffusion) + rdd <- projectRDD(params, rdi = rdi, + species_params = params@species_params, t = 0) + } else { + f_rdi <- get(params@rates_funcs$RDI) + rdi <- f_rdi(params, n = N_in, n_pp = npp_in, + n_other = n_other, t = 0, + e_repro = r$e_repro, e_growth = r$e_growth, + mort = r$mort, diffusion = r$diffusion) + f_rdd <- get(params@rates_funcs$RDD) + rdd <- f_rdd(rdi = rdi, + species_params = params@species_params, + params = params, t = 0) + } + } else { + rdd <- rdd_const + } + + coefs <- get_transport_coefs(params, n = N_in, g = r$e_growth, + mu = r$mort, dt = 1, + recruitment_flux = rdd, + d = r$diffusion, + flux_limiter = flux_limiter) + N_out <- project_n_loop(N_in, coefs$a, coefs$b, coefs$c, coefs$S, + params@w_min_idx) + # Return both the updated N and the rates (needed for resource step). + list(N = zero_above_support(N_out, support_top_idx(params)), + rates = r) + } + + # ------------------------------------------------------------------------- + # Steady-state arrays + # ------------------------------------------------------------------------- + N_ss <- params@initial_n + N_vec <- N_ss[active$mask] + n_fish_active <- length(N_vec) + + if (!include_resource) { + # -- Reduced system: resource at quasi-static equilibrium -------------- + if (params@resource_dynamics != "resource_semichemostat") { + stop("getStability() with include_resource = FALSE requires ", + "semichemostat resource dynamics. ", + "Use include_resource = TRUE for other resource dynamics.") + } + + reduced_step <- function(N_in) { + npp_in <- resource_steady_semichemostat(params, N_in, n_other) + fish_step(N_in, npp_in)$N + } + + n_state <- n_fish_active + L <- matrix(0, nrow = n_state, ncol = n_state) + for (k in seq_len(n_fish_active)) { + eps_k <- h * pmax(abs(N_vec[k]), .Machine$double.eps) + N_plus <- N_ss; N_plus[active_idx[k]] <- N_vec[k] + eps_k + N_minus <- N_ss; N_minus[active_idx[k]] <- N_vec[k] - eps_k + L[, k] <- (reduced_step(N_plus)[active$mask] - + reduced_step(N_minus)[active$mask]) / (2 * eps_k) + } + + } else { + # -- Full coupled system: perturb fish and resource independently ------ + npp_ss <- params@initial_n_pp + npp_vec <- as.numeric(npp_ss) + n_npp <- length(npp_vec) + n_state <- n_fish_active + n_npp + + resource_dyn <- get(params@resource_dynamics) + + # Full one-step map: returns c(N_new[active_mask], npp_new_all) + full_step <- function(N_in, npp_in) { + res <- fish_step(N_in, npp_in) + # Resource dynamics with actual rates (includes resource_mort) + npp_out <- resource_dyn(params, N_in, npp_in, n_other, + rates = res$rates, t = 0, dt = 1, + resource_rate = params@rr_pp, + resource_capacity = params@cc_pp) + c(res$N[active$mask], as.numeric(npp_out)) + } + + n_state <- n_fish_active + n_npp + L <- matrix(0, nrow = n_state, ncol = n_state) + + # Columns 1..n_fish_active: perturb fish, resource held at n_pp_ss + for (k in seq_len(n_fish_active)) { + eps_k <- h * pmax(abs(N_vec[k]), .Machine$double.eps) + N_plus <- N_ss; N_plus[active_idx[k]] <- N_vec[k] + eps_k + N_minus <- N_ss; N_minus[active_idx[k]] <- N_vec[k] - eps_k + L[, k] <- (full_step(N_plus, npp_ss) - + full_step(N_minus, npp_ss)) / (2 * eps_k) + } + + # Columns n_fish_active+1..n_state: perturb resource, fish held at N_ss + for (j in seq_len(n_npp)) { + k <- n_fish_active + j + eps_j <- h * pmax(abs(npp_vec[j]), .Machine$double.eps) + npp_plus <- npp_ss; npp_plus[j] <- npp_vec[j] + eps_j + npp_minus <- npp_ss; npp_minus[j] <- npp_vec[j] - eps_j + L[, k] <- (full_step(N_ss, npp_plus) - + full_step(N_ss, npp_minus)) / (2 * eps_j) + } + } + + # ------------------------------------------------------------------------- + # Eigenvalue analysis (shared by both branches) + # ------------------------------------------------------------------------- + evals_map <- eigen(L, only.values = TRUE)$values + ord <- order(Mod(evals_map), decreasing = TRUE) + evals_map <- evals_map[ord] + + spectral_radius <- max(Mod(evals_map)) + stable <- spectral_radius < 1 + + lam1 <- evals_map[1] + theta1 <- Arg(lam1) + if (abs(theta1) < 1e-10) { + dominant_period <- Inf + } else if (abs(abs(theta1) - pi) < 1e-10) { + dominant_period <- NA + } else { + dominant_period <- 2 * pi / abs(theta1) + } + + is_complex <- abs(Im(evals_map)) > 1e-8 + if (any(is_complex)) { + complex_evals <- evals_map[is_complex] + closest_idx <- which.min(abs(Mod(complex_evals) - 1)) + lam_hopf <- complex_evals[closest_idx] + hopf_period <- 2 * pi / abs(Arg(lam_hopf)) + } else { + hopf_period <- NULL + } + + list( + eigenvalues = evals_map, + spectral_radius = spectral_radius, + stable = stable, + dominant_period = dominant_period, + hopf_period = hopf_period, + n_active = n_state + ) +} diff --git a/man/getStability.Rd b/man/getStability.Rd new file mode 100644 index 000000000..50246e6cf --- /dev/null +++ b/man/getStability.Rd @@ -0,0 +1,107 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/steadyNewton.R +\name{getStability} +\alias{getStability} +\title{Analyse the dynamic stability of a mizer steady state} +\usage{ +getStability( + params, + reproduction = c("fixed", "dynamic"), + effort = params@initial_effort, + include_resource = FALSE, + extinction_floor = 1e-06, + h = 1e-04 +) +} +\arguments{ +\item{params}{A \linkS4class{MizerParams} object whose \code{initial_n} holds the +steady state to analyse. Typically the output of \code{\link[=steadyNewton]{steadyNewton()}}.} + +\item{reproduction}{Whether the reproduction rate is held fixed (\code{"fixed"}, +default) or run dynamically (\code{"dynamic"}) during the one-step evaluation. +Must match the choice used when the steady state was computed.} + +\item{effort}{The fishing effort to use. By default the initial effort +stored in \code{params}.} + +\item{include_resource}{If \code{FALSE} (default) the resource is treated as a +quasi-static fast variable: for each perturbed fish abundance the resource +is set to its analytic steady-state value conditioned on that fish +abundance (valid only for semichemostat resource dynamics). If \code{TRUE}, +both fish and resource cells are perturbed independently and the resource +is evolved with the full resource dynamics function stored in +\code{params@resource_dynamics}, giving the complete coupled Jacobian.} + +\item{extinction_floor}{Relative abundance floor for the dynamic reproduction +case. Default is \code{1e-6}.} + +\item{h}{Relative step size for centred finite differences. Default \code{1e-4}.} +} +\value{ +A named list with the following components: +\describe{ +\item{\code{eigenvalues}}{Complex vector of all eigenvalues of the +linearised Jacobian, sorted by decreasing modulus.} +\item{\code{spectral_radius}}{The largest modulus \eqn{\max_i|\lambda_i|}.} +\item{\code{stable}}{Logical: \code{TRUE} when \code{spectral_radius < 1}.} +\item{\code{dominant_period}}{The period (in time steps) of the dominant +eigenvalue: \code{2*pi / abs(Arg(lambda_1))}. \code{Inf} for a real positive +dominant eigenvalue (monotone dynamics); \code{NA} for a period-2 flip.} +\item{\code{hopf_period}}{Period (in time steps) of the complex eigenvalue +whose modulus is closest to 1; \code{NULL} when no complex eigenvalue +exists. This is the expected limit-cycle period near a Hopf +bifurcation.} +\item{\code{n_active}}{Dimension of the Jacobian: number of active fish cells +when \code{include_resource = FALSE}, or fish cells plus all resource cells +when \code{include_resource = TRUE}.} +} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +Computes the eigenvalues of the linearised one-step-ahead map at the steady +state stored in \code{params@initial_n}. These eigenvalues determine whether +the steady state is dynamically stable and, when a Hopf bifurcation is +approached, the period of the emergent limit cycle. +} +\details{ +\subsection{Mathematical background}{ + +The mizer time step applies a backward-Euler transport solve for the fish: +\deqn{A(N^t, n_{pp}^t)\,N^{t+1} = S(N^t, n_{pp}^t),} +and an exact semi-chemostat update for the resource: +\deqn{n_{pp}^{t+1} = n_{pp}^* + (n_{pp}^t - n_{pp}^*)\,e^{-\mu^t\,dt},} +where \eqn{n_{pp}^* = r_{pp}\,c_{pp}/\mu^t} is the resource steady state +conditioned on the mortality \eqn{\mu^t} due to consumers at time \eqn{t}. + +The stability is determined by the Jacobian of the full one-step-ahead map +\eqn{G : (N, n_{pp}) \mapsto (N^{t+1}, n_{pp}^{t+1})} at the fixed point. + +When \code{include_resource = FALSE} (the default), the resource is treated as a +fast variable that adjusts \emph{instantaneously} to the consumer abundance: for +each perturbed \eqn{N}, \eqn{n_{pp}} is set to its quasi-static equilibrium +\eqn{n_{pp}^*(N)}. The resulting reduced Jacobian \eqn{L_{\text{red}}} has +dimension equal to the number of active fish cells. This is equivalent to +projecting the full dynamics onto the slow manifold \eqn{n_{pp} = n_{pp}^*(N)}. + +When \code{include_resource = TRUE}, both fish and resource cells are perturbed +independently and the full coupled Jacobian \eqn{L_{\text{full}}} is +returned. Its eigenvalues include both the slow fish modes and a cluster of +fast resource-relaxation modes (with modulus \eqn{e^{-\mu\,dt} \ll 1}). +Comparing the dominant eigenvalues of the two analyses shows how much the +quasi-static approximation affects the stability conclusion. + +The steady state is \strong{stable} when all eigenvalues satisfy +\eqn{|\lambda_i| < 1} and \strong{unstable} when at least one exceeds 1. + +A \strong{Hopf bifurcation} occurs when a complex-conjugate pair of eigenvalues +crosses the unit circle, giving a limit-cycle period +\deqn{T = \frac{2\pi}{|\arg(\lambda)|} \text{ time steps.}} +For the default mizer time step of one year this is in years. + +Both branches use the same \code{project_n_loop()} C++ Thomas solver as the +regular dynamics and centred finite differences with relative step \eqn{h}. +} +} +\seealso{ +\code{\link[=steadyNewton]{steadyNewton()}} +} diff --git a/man/steadyNewton.Rd b/man/steadyNewton.Rd new file mode 100644 index 000000000..2f45f8b6b --- /dev/null +++ b/man/steadyNewton.Rd @@ -0,0 +1,136 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/steadyNewton.R +\name{steadyNewton} +\alias{steadyNewton} +\alias{steadyNewton.MizerParams} +\title{Find a steady state by directly solving the steady-state equation} +\usage{ +steadyNewton(params, ...) + +\method{steadyNewton}{MizerParams}( + params, + effort = params@initial_effort, + preserve = c("reproduction_level", "erepro", "R_max"), + reproduction = c("fixed", "dynamic"), + extinction_floor = 1e-06, + verbose = FALSE, + tol = 1e-06, + maxit = 200, + method = c("Broyden", "Newton"), + global = "dbldog", + stability = FALSE, + ... +) +} +\arguments{ +\item{params}{A \linkS4class{MizerParams} object. Its \code{initial_n} is used as +the starting guess for the iteration.} + +\item{...}{Unused.} + +\item{effort}{The fishing effort. By default the initial effort stored in +\code{params}.} + +\item{preserve}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +Specifies whether the \code{reproduction_level} should be preserved (default) +or the maximum reproduction rate \code{R_max} or the reproductive efficiency +\code{erepro}. See \code{\link[=setBevertonHolt]{setBevertonHolt()}} for an explanation of the +\code{reproduction_level}. This argument is ignored when \code{reproduction = "dynamic"}.} + +\item{reproduction}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +If \code{"fixed"}, the reproduction rate (RDD) is held constant at the initial +value. If \code{"dynamic"}, the reproduction dynamics are run dynamically and the +reproduction parameters are not adjusted. Default is \code{"fixed"}.} + +\item{extinction_floor}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +The relative abundance floor below which a species is considered extinct. +Only used when \code{reproduction = "dynamic"}. Default is 1e-6.} + +\item{verbose}{If \code{TRUE} then the solver iterations will be traced and +printed to the console. Default is \code{FALSE}.} + +\item{tol}{Convergence tolerance passed to \code{\link[nleqslv:nleqslv]{nleqslv::nleqslv()}} (both the +function-value tolerance \code{ftol} and the step tolerance \code{xtol}).} + +\item{maxit}{Maximum number of iterations for \code{\link[nleqslv:nleqslv]{nleqslv::nleqslv()}}.} + +\item{method}{The \code{\link[nleqslv:nleqslv]{nleqslv::nleqslv()}} method, either \code{"Newton"} (with a +numerical Jacobian calculated at each iteration) or \code{"Broyden"} (which +calculates the full Jacobian only once and then only updates it on each +iteration. '"Broyden"' is the default.} + +\item{global}{The globalisation strategy passed to \code{\link[nleqslv:nleqslv]{nleqslv::nleqslv()}}. +The default \code{"dbldog"} (double dogleg) is a robust trust-region method.} + +\item{stability}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +If \code{TRUE}, \code{\link[=getStability]{getStability()}} is called after convergence and its result is +stored as the attribute \code{"stability"} on the returned +\linkS4class{MizerParams} object. Default is \code{FALSE}.} +} +\value{ +A \linkS4class{MizerParams} object with the initial state set to the +steady state. If \code{stability = TRUE}, the object carries the attribute +\code{"stability"} containing the list returned by \code{\link[=getStability]{getStability()}}. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +This is an alternative to \code{\link[=steady]{steady()}} that finds the steady state by solving +the steady-state equation \code{F(N) = 0} with a Newton-type root finder instead +of running the dynamics to convergence. The advantage is that it converges +to the steady state even when that steady state is dynamically \strong{unstable}, +a case in which \code{\link[=steady]{steady()}} fails because the time-stepping diverges away from +the fixed point. +} +\details{ +By default, or when \code{reproduction = "fixed"}, the function holds the +reproduction rate (RDD) constant while solving for the consumer spectra, +substitutes the analytic steady state of the resource, and keeps any other +components constant. After the spectra have been found it restores +density-dependent Beverton-Holt reproduction with \code{\link[=setBevertonHolt]{setBevertonHolt()}}, +honouring the \code{preserve} argument exactly as \code{\link[=steady]{steady()}} does. +If \code{reproduction = "dynamic"}, the reproduction dynamics are run dynamically, +meaning the reproduction rate varies during the solve and the reproduction +parameters are not adjusted. + +The consumer densities are solved for in log space, which both keeps them +positive and conditions the otherwise badly-scaled system. The unknowns are +the size classes in the full potential grid support (running from the egg +size up to the grid truncation limit \code{support_top_idx()}), regardless of +whether they carry non-zero density in the supplied initial spectra. A +smoothed log-abundance penalty floor is applied to all size classes, which +automatically handles zero-density and tail classes smoothly and prevents +singular Jacobians. After convergence, size classes that remain at or near +the floor are set to zero. This allows the solver to automatically discover +the support of the steady state. The nonlinear system is solved with a +globalised Newton iteration from the \code{nleqslv} package, starting from the +current \code{initial_n}. Newton's method converges from any starting point in the +\emph{root's} basin of attraction (which is unrelated to the dynamic stability of +the steady state), so a reasonable initial guess should still be supplied in +\code{initialN(params)} — for example the spectra from a nearby stable +parameterisation, or the (diverging) output of \code{\link[=steady]{steady()}}. + +The solver respects the active transport scheme: if the experimental +second-order scheme is enabled (see \code{\link[=second_order_w]{second_order_w()}}) it solves the +steady-state equation of that scheme. With the van Leer reconstruction the +residual is only Lipschitz, so the iteration converges to a fixed point of the +dynamics but not to machine precision. The unlimited \code{"centred"} +reconstruction admits an undamped odd-even mode at a steady state with no +physical diffusion, giving an ill-conditioned steady-state Jacobian for which +the solver is not expected to converge. + +Only the default semichemostat resource dynamics +(\code{resource_dynamics = "resource_semichemostat"}) are currently supported, +because the solver substitutes the analytic resource steady state. For other +resource dynamics the function stops with an error. +} +\examples{ +\donttest{ +params <- steadyNewton(NS_params) +plotSpectra(params) +params <- steadyNewton(NS_params, stability = TRUE) +attr(params, "stability")$stable +} +} +\seealso{ +\code{\link[=steady]{steady()}}, \code{\link[=steadySingleSpecies]{steadySingleSpecies()}}, \code{\link[=getStability]{getStability()}} +} diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index 293c0fe2e..de3db4f7a 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -142,6 +142,8 @@ reference: walks through this calibration workflow. contents: - steady + - steadyNewton + - getStability - steadySingleSpecies - matchGrowth - plotBiomassObservedVsModel diff --git a/tests/testthat/test-steadyNewton.R b/tests/testthat/test-steadyNewton.R new file mode 100644 index 000000000..2c50eaf7c --- /dev/null +++ b/tests/testthat/test-steadyNewton.R @@ -0,0 +1,322 @@ +# Tests for the direct (Newton) steady-state solver. + +# A steadied model is a good starting guess for the root finder. +p_steady <- steady(NS_params_small, t_max = 50, progress_bar = FALSE) + +test_that("steadyNewton returns a fixed point of the dynamics", { + pn <- steadyNewton(p_steady) + expect_s4_class(pn, "MizerParams") + + # The returned spectra should barely move when projected forward. + sim <- project(pn, t_max = 1, dt = 0.25, t_save = 1) + n0 <- pn@initial_n + n1 <- finalN(sim) + support <- n0 > 0 + drift <- max((abs(n1 - n0) / n0)[support]) + expect_lt(drift, 1e-3) +}) + +test_that("steadyNewton solves the discrete steady-state equation to tolerance", { + # Solve the held-RDD problem directly and check the residual at the root. + params <- validParams(p_steady) + rdd_const <- getRDD(params) + active <- mizer:::steady_active_set(params) + resfn <- mizer:::steady_state_residual(params, rdd_const, + params@initial_n_other, + params@initial_effort, active) + x0 <- log(params@initial_n[active$mask]) + sol <- nleqslv::nleqslv(x0, resfn, method = "Newton", global = "dbldog", + control = list(maxit = 100, ftol = 1e-10, + xtol = 1e-10)) + expect_lte(sol$termcd, 2) + expect_lt(max(abs(sol$fvec)), 1e-8) +}) + +test_that("steadyNewton agrees with steady on a stable model", { + pn <- steadyNewton(p_steady) + # Compare where the steady() spectra are well above zero. + # steadyNewton finds a tighter fixed point than steady() reaches by + # time-stepping, so a few percent difference on the well-resolved bins is + # expected (and is the point of the exercise). + big <- p_steady@initial_n > max(p_steady@initial_n) * 1e-6 + rel <- abs(pn@initial_n - p_steady@initial_n) / p_steady@initial_n + expect_lt(max(rel[big]), 0.1) +}) + +test_that("steadyNewton honours the preserve and reproduction arguments", { + pn_level <- steadyNewton(p_steady, preserve = "reproduction_level") + expect_equal(getReproductionLevel(pn_level), + getReproductionLevel(p_steady), tolerance = 1e-5) + + pn_rmax <- steadyNewton(p_steady, preserve = "R_max") + expect_equal(pn_rmax@species_params$R_max, + p_steady@species_params$R_max, tolerance = 1e-5) + + pn_erepro <- suppressWarnings(steadyNewton(p_steady, preserve = "erepro")) + expect_equal(pn_erepro@species_params$erepro, + p_steady@species_params$erepro, tolerance = 1e-5) + + pn_none <- steadyNewton(p_steady, reproduction = "dynamic") + expect_equal(pn_none@species_params$R_max, + p_steady@species_params$R_max, tolerance = 1e-5) + expect_equal(pn_none@species_params$erepro, + p_steady@species_params$erepro, tolerance = 1e-5) + + # Verify preserve argument is ignored when reproduction = "dynamic" + pn_ignored <- steadyNewton(p_steady, reproduction = "dynamic", preserve = "invalid_option") + expect_equal(pn_ignored@species_params$R_max, + p_steady@species_params$R_max, tolerance = 1e-5) + + # Verify verbose = TRUE captures iteration report output + out <- capture.output(steadyNewton(p_steady, verbose = TRUE)) + expect_true(any(grepl("Iteration report", out))) +}) + +test_that("steadyNewton handles extinctions under reproduction = 'dynamic' with relative floor", { + # Make species 3 (Cod) unviable by setting its reproduction efficiency extremely low + p_extinct <- p_steady + p_extinct@species_params$erepro[3] <- 1e-12 + p_extinct@initial_n[3, ] <- p_steady@initial_n[3, ] + + # Verify that the solver issues the extinction warning and pegs to the floor + expect_warning(pn_ext <- steadyNewton(p_extinct, reproduction = "dynamic", extinction_floor = 1e-6), + "went extinct and were pegged to their abundance floor") + + # Verify Cod abundance is pegged exactly to the floor (1e-6 of its initial abundance) + lo <- p_extinct@w_min_idx[3] + ratio <- pn_ext@initial_n[3, lo] / p_extinct@initial_n[3, lo] + expect_equal(ratio, 1e-6, tolerance = 1e-2) +}) + +test_that("steadyNewton errors for unsupported resource dynamics", { + p_log <- setResource(NS_params_small, + resource_dynamics = "resource_logistic") + expect_error(steadyNewton(p_log), "semichemostat") +}) + +test_that("steadyNewton works with the second-order (van Leer) scheme", { + p <- NS_params_small + sow <- second_order_w(p) + sow$flux <- "van_leer" + sow$bin_average <- TRUE + second_order_w(p) <- sow + ps <- steady(p, t_max = 100, progress_bar = FALSE) + + pn <- steadyNewton(ps) + expect_s4_class(pn, "MizerParams") + # The returned state must still be a fixed point of the (second-order) + # dynamics. The van Leer limiter is only Lipschitz, so the Newton residual + # cannot be driven to machine precision, but the projected drift is the + # honest test and stays small. + sim <- project(pn, t_max = 1, dt = 0.25, t_save = 1) + n0 <- pn@initial_n + n1 <- finalN(sim) + support <- n0 > 0 + drift <- max((abs(n1 - n0) / n0)[support]) + expect_lt(drift, 1e-3) +}) + +test_that("support_top_idx drops the pile-up bin for the second-order scheme", { + # First-order upwind feeds a class from the growth of the class below, so the + # support reaches one bin past w_max; the second-order scheme feeds a class + # from the growth at its own lower face, so it stops at w_max. + p1 <- NS_params_small + p2 <- NS_params_small + sow <- second_order_w(p2) + sow$flux <- "van_leer" + second_order_w(p2) <- sow + + w_max_idx <- sapply(seq_len(nrow(p1@species_params)), function(i) { + sum(p1@w <= p1@species_params$w_max[i]) + }) + no_w <- length(p1@w) + expect_equal(mizer:::support_top_idx(p1), pmin(w_max_idx + 1L, no_w), + ignore_attr = TRUE) + expect_equal(mizer:::support_top_idx(p2), pmin(w_max_idx, no_w), + ignore_attr = TRUE) +}) + +test_that("steady_active_set always reaches the grid truncation limit", { + # Under the new dynamic support design, the active set always reaches + # the grid truncation limit (support_top_idx) to allow the solver to + # automatically discover the non-zero region. + p <- NS_params_small + no_w <- length(p@w) + grid_top <- mizer:::support_top_idx(p) + + active0 <- mizer:::steady_active_set(p) + expect_equal(active0$w_top, unname(grid_top)) + + # Even if we zero out the tail of the abundances, the mask still reaches grid_top + cutoff <- unname(pmax(p@w_min_idx + 2L, grid_top - 3L)) + for (i in seq_len(nrow(p@species_params))) { + if (cutoff[i] < no_w) { + p@initial_n[i, (cutoff[i] + 1L):no_w] <- 0 + } + } + active <- mizer:::steady_active_set(p) + expect_equal(active$w_top, unname(grid_top)) +}) + +test_that("support_top_idx is the first class above w_max", { + p <- NS_params_small + g <- getEGrowth(p) + w_top <- mizer:::support_top_idx(p) + no_w <- length(p@w) + for (i in seq_len(nrow(p@species_params))) { + # The support top is one class above w_max (the pile-up bin), capped at + # the top of the grid. + w_max_idx <- sum(p@w <= p@species_params$w_max[i]) + expect_equal(unname(w_top[i]), min(w_max_idx + 1, no_w)) + # For a standard model this coincides with the first zero-growth class. + if (w_top[i] < no_w) { + expect_equal(unname(g[i, w_top[i]]), 0) + expect_gt(g[i, w_top[i] - 1], 0) + } + } +}) + +test_that("the upper boundary condition severs coupling above the growth cutoff", { + # The transport coefficients sever the backward coupling from the + # growth-chain top to the class above it (c = 0 there). Together with the + # vanishing sub-diagonal (a = 0 at top+1, because growth out of the top + # class is zero) this isolates the active spectrum from the inactive tail in + # the tridiagonal solve, so the active solution never reads the tail. The + # tail itself is then held at zero by zero_above_support(), which is what + # stops a diffusive leak above the maximum size. + p <- NS_params_small + g <- getEGrowth(p) + mu <- getMort(p) + d <- getDiffusion(p) + no_w <- length(p@w) + w_top <- mizer:::support_top_idx(p) + coefs <- mizer:::get_transport_coefs(p, p@initial_n, g, mu, dt = 1, + recruitment_flux = getRDD(p), d = d, + flux_limiter = "none") + for (i in seq_len(nrow(p@species_params))) { + # Backward sweep isolated: top class does not read the class above. + expect_equal(unname(coefs$c[i, w_top[i]]), 0) + # Forward sweep isolated: nothing grows out of the top class. + if (w_top[i] < no_w) { + expect_equal(unname(coefs$a[i, w_top[i] + 1]), 0) + } + } +}) + +test_that("the upper boundary condition stops diffusion leaking above w_max", { + # On the full model (where the cutoff is a genuine maturity cutoff) a step + # with diffusion leaves the density exactly zero above the growth-chain top. + p <- setParams(NS_params, use_predation_diffusion = TRUE) + no_w <- length(p@w) + w_top <- mizer:::support_top_idx(p) + sim <- project(p, t_max = 0.5, dt = 0.1, t_save = 0.5) + nf <- finalN(sim) + for (i in seq_len(nrow(p@species_params))) { + if (w_top[i] < no_w) { + expect_true(all(nf[i, (w_top[i] + 1):no_w] == 0)) + } + } +}) + +test_that("steadyNewton handles initial guesses that are non-zero at large sizes where steady state is zero", { + # Start with the stable model + p <- p_steady + # Fill the trailing tail (which is zero in p_steady) with positive numbers + grid_top <- mizer:::support_top_idx(p) + no_w <- length(p@w) + + # We will corrupt the tail of the first species with positive numbers + # above its actual support. + # In p_steady, species 1 (Sprat) only grows to its w_max (0.33g). + # We set non-zero values up to the end of the grid. + idx_zeros <- (grid_top[1] + 1):no_w + if (length(idx_zeros) > 0) { + p@initial_n[1, idx_zeros] <- 1e-3 + } + + # Run steadyNewton + pn <- steadyNewton(p) + expect_s4_class(pn, "MizerParams") + + # The tail should be correctly zeroed out in the result + if (length(idx_zeros) > 0) { + expect_true(all(pn@initial_n[1, idx_zeros] == 0)) + } +}) + +# Tests for getStability() --------------------------------------------------- + +test_that("getStability returns a well-formed list for a stable model", { + pn <- steadyNewton(p_steady) + stab <- getStability(pn) + + expect_type(stab, "list") + expect_named(stab, c("eigenvalues", "spectral_radius", "stable", + "dominant_period", "hopf_period", "n_active")) + expect_true(is.complex(stab$eigenvalues)) + expect_length(stab$eigenvalues, stab$n_active) + expect_true(is.numeric(stab$spectral_radius)) + expect_true(is.logical(stab$stable)) +}) + +test_that("getStability reports stable = TRUE for the NS model at its steady state", { + pn <- steadyNewton(p_steady) + stab <- getStability(pn) + + expect_true(stab$stable) + expect_lt(stab$spectral_radius, 1) +}) + +test_that("getStability eigenvalues are consistent with spectral_radius", { + pn <- steadyNewton(p_steady) + stab <- getStability(pn) + + expect_equal(stab$spectral_radius, max(Mod(stab$eigenvalues))) +}) + +test_that("steadyNewton attaches stability when stability = TRUE", { + pn <- steadyNewton(p_steady, stability = TRUE) + + stab <- attr(pn, "stability") + expect_type(stab, "list") + expect_true(stab$stable) +}) + +test_that("getStability hopf_period is NULL when all eigenvalues are real", { + # For a 1-species model the dominant eigenvalue is typically real. + # Use a simple single-species model as a proxy: steadyNewton on it, + # then check that if hopf_period is not NULL it is a positive finite number. + pn <- steadyNewton(p_steady) + stab <- getStability(pn) + + if (!is.null(stab$hopf_period)) { + expect_true(is.finite(stab$hopf_period)) + expect_gt(stab$hopf_period, 0) + } +}) + +test_that("getStability with include_resource = TRUE returns well-formed list", { + pn <- steadyNewton(p_steady) + stab_full <- getStability(pn, include_resource = TRUE) + + expect_type(stab_full, "list") + expect_named(stab_full, c("eigenvalues", "spectral_radius", "stable", + "dominant_period", "hopf_period", "n_active")) + # n_active must equal n_fish_active + n_resource (always strictly larger) + stab_red <- getStability(pn, include_resource = FALSE) + expect_gt(stab_full$n_active, stab_red$n_active) + expect_length(stab_full$eigenvalues, stab_full$n_active) + expect_true(stab_full$stable) +}) + +test_that("full and reduced stability analyses agree on spectral radius for NS model", { + # For a model where the resource dynamics are fast (large rr_pp), the + # quasi-static reduction should give almost the same dominant eigenvalue. + pn <- steadyNewton(p_steady) + stab_red <- getStability(pn, include_resource = FALSE) + stab_full <- getStability(pn, include_resource = TRUE) + + # Dominant eigenvalue moduli should match to within a loose tolerance. + expect_equal(stab_full$spectral_radius, stab_red$spectral_radius, + tolerance = 0.05) +}) From f2defe3ec28cba15a472e60ca2bc7f7b03763e57 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 08:35:41 +0100 Subject: [PATCH 02/15] getStability: return leading eigenvectors reshaped to (species x size) arrays - eigen(L, only.values=TRUE) -> eigen(L) to capture eigenvectors - Eigenvectors sorted alongside eigenvalues by decreasing modulus - Top 2 eigenvectors reshaped to complex array (n_sp x n_w x 2) with same dimnames as params@initial_n; normalised to max Mod = 1 - For include_resource=TRUE: list with $fish array and $resource matrix (n_w_full x 2) for the resource component - Docs and tests updated --- R/steadyNewton.R | 71 +++++++++++++++++++++++++++--- man/getStability.Rd | 11 +++++ tests/testthat/test-steadyNewton.R | 24 +++++++++- 3 files changed, 97 insertions(+), 9 deletions(-) diff --git a/R/steadyNewton.R b/R/steadyNewton.R index 86ec3652f..7c880b620 100644 --- a/R/steadyNewton.R +++ b/R/steadyNewton.R @@ -541,6 +541,17 @@ steady_state_residual <- function(params, rdd_const, n_other, effort, active, #' \item{`n_active`}{Dimension of the Jacobian: number of active fish cells #' when `include_resource = FALSE`, or fish cells plus all resource cells #' when `include_resource = TRUE`.} +#' \item{`leading_eigenvectors`}{The eigenvectors of the two largest-modulus +#' eigenvalues, reshaped back into the fish abundance space. +#' When `include_resource = FALSE`: a complex array of shape +#' `(n_species, n_sizes, 2)` with the same species and size dimnames as +#' `params@initial_n`. When `include_resource = TRUE`: a list with +#' `$fish` (the same array) and `$resource` (a complex matrix of shape +#' `(n_w_full, 2)` for the resource component). +#' Each eigenvector is normalised so that its maximum modulus equals 1. +#' The real and imaginary parts of eigenvector 1 span the two-dimensional +#' oscillation plane of the dominant mode; `Mod()` gives the amplitude +#' pattern across species and sizes.} #' } #' @seealso [steadyNewton()] #' @export @@ -688,9 +699,14 @@ getStability <- function(params, # ------------------------------------------------------------------------- # Eigenvalue analysis (shared by both branches) # ------------------------------------------------------------------------- - evals_map <- eigen(L, only.values = TRUE)$values + eig <- eigen(L) + evals_map <- eig$values + evecs <- eig$vectors # columns are eigenvectors + + # Sort by decreasing eigenvalue modulus. ord <- order(Mod(evals_map), decreasing = TRUE) evals_map <- evals_map[ord] + evecs <- evecs[, ord, drop = FALSE] spectral_radius <- max(Mod(evals_map)) stable <- spectral_radius < 1 @@ -715,12 +731,53 @@ getStability <- function(params, hopf_period <- NULL } + # ------------------------------------------------------------------------- + # Reshape leading eigenvectors back to (species x size) complex arrays. + # Eigenvectors are normalised so max(Mod(.)) == 1 for comparability. + # For include_resource = TRUE the resource component is also extracted. + # ------------------------------------------------------------------------- + n_leading <- min(2L, n_state) # top 2 (or 1 if Jacobian is 1×1) + no_sp <- nrow(N_ss) + no_w <- ncol(N_ss) + + # Preallocate leading eigenvector array (species x size x n_leading), complex + leading_evecs_fish <- array(0 + 0i, + dim = c(no_sp, no_w, n_leading), + dimnames = c(dimnames(N_ss), list(NULL))) + + for (k in seq_len(n_leading)) { + v_fish <- evecs[seq_len(n_fish_active), k] # fish component + # Normalise to max modulus 1 + max_mod <- max(Mod(v_fish)) + if (max_mod > 0) v_fish <- v_fish / max_mod + M <- matrix(0 + 0i, nrow = no_sp, ncol = no_w) + M[active_idx] <- v_fish + leading_evecs_fish[, , k] <- M + } + + if (include_resource) { + # Resource component: rows n_fish_active+1..n_state of each eigenvector. + leading_evecs_resource <- matrix(0 + 0i, nrow = n_npp, ncol = n_leading) + for (k in seq_len(n_leading)) { + v_res <- evecs[n_fish_active + seq_len(n_npp), k] + max_mod <- max(Mod(v_res)) + if (max_mod > 0) v_res <- v_res / max_mod + leading_evecs_resource[, k] <- v_res + } + rownames(leading_evecs_resource) <- names(params@initial_n_pp) + leading_eigenvectors <- list(fish = leading_evecs_fish, + resource = leading_evecs_resource) + } else { + leading_eigenvectors <- leading_evecs_fish + } + list( - eigenvalues = evals_map, - spectral_radius = spectral_radius, - stable = stable, - dominant_period = dominant_period, - hopf_period = hopf_period, - n_active = n_state + eigenvalues = evals_map, + spectral_radius = spectral_radius, + stable = stable, + dominant_period = dominant_period, + hopf_period = hopf_period, + n_active = n_state, + leading_eigenvectors = leading_eigenvectors ) } diff --git a/man/getStability.Rd b/man/getStability.Rd index 50246e6cf..71222bc01 100644 --- a/man/getStability.Rd +++ b/man/getStability.Rd @@ -54,6 +54,17 @@ bifurcation.} \item{\code{n_active}}{Dimension of the Jacobian: number of active fish cells when \code{include_resource = FALSE}, or fish cells plus all resource cells when \code{include_resource = TRUE}.} +\item{\code{leading_eigenvectors}}{The eigenvectors of the two largest-modulus +eigenvalues, reshaped back into the fish abundance space. +When \code{include_resource = FALSE}: a complex array of shape +\verb{(n_species, n_sizes, 2)} with the same species and size dimnames as +\code{params@initial_n}. When \code{include_resource = TRUE}: a list with +\verb{$fish} (the same array) and \verb{$resource} (a complex matrix of shape +\verb{(n_w_full, 2)} for the resource component). +Each eigenvector is normalised so that its maximum modulus equals 1. +The real and imaginary parts of eigenvector 1 span the two-dimensional +oscillation plane of the dominant mode; \code{Mod()} gives the amplitude +pattern across species and sizes.} } } \description{ diff --git a/tests/testthat/test-steadyNewton.R b/tests/testthat/test-steadyNewton.R index 2c50eaf7c..cc661ed3c 100644 --- a/tests/testthat/test-steadyNewton.R +++ b/tests/testthat/test-steadyNewton.R @@ -252,7 +252,8 @@ test_that("getStability returns a well-formed list for a stable model", { expect_type(stab, "list") expect_named(stab, c("eigenvalues", "spectral_radius", "stable", - "dominant_period", "hopf_period", "n_active")) + "dominant_period", "hopf_period", "n_active", + "leading_eigenvectors")) expect_true(is.complex(stab$eigenvalues)) expect_length(stab$eigenvalues, stab$n_active) expect_true(is.numeric(stab$spectral_radius)) @@ -301,7 +302,8 @@ test_that("getStability with include_resource = TRUE returns well-formed list", expect_type(stab_full, "list") expect_named(stab_full, c("eigenvalues", "spectral_radius", "stable", - "dominant_period", "hopf_period", "n_active")) + "dominant_period", "hopf_period", "n_active", + "leading_eigenvectors")) # n_active must equal n_fish_active + n_resource (always strictly larger) stab_red <- getStability(pn, include_resource = FALSE) expect_gt(stab_full$n_active, stab_red$n_active) @@ -320,3 +322,21 @@ test_that("full and reduced stability analyses agree on spectral radius for NS m expect_equal(stab_full$spectral_radius, stab_red$spectral_radius, tolerance = 0.05) }) + +test_that("leading_eigenvectors have correct shape and are normalised", { + pn <- steadyNewton(p_steady) + stab <- getStability(pn) + lev <- stab$leading_eigenvectors + + # Shape: (n_species, n_sizes, 2) + expect_equal(dim(lev)[1], nrow(pn@initial_n)) + expect_equal(dim(lev)[2], ncol(pn@initial_n)) + expect_equal(dim(lev)[3], 2L) + + # Each eigenvector normalised: max(Mod(.)) == 1 over active cells + expect_equal(max(Mod(lev[, , 1])), 1, tolerance = 1e-10) + expect_equal(max(Mod(lev[, , 2])), 1, tolerance = 1e-10) + + # Dimnames match initial_n + expect_equal(dimnames(lev)[1:2], dimnames(pn@initial_n)) +}) From dae4e5a61089907383c117c3975e314760c74121 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 08:47:14 +0100 Subject: [PATCH 03/15] Add getLimitCycleSim(): linearised limit cycle as a MizerSim object Takes the output of steadyNewton() (with stability analysis attached) and constructs a MizerSim covering one period of the limit cycle in the linear approximation: N(t) = max(N* + A * Re[exp(i*theta*t) * v], 0) where: - v is the leading complex eigenvector from getStability() - theta = 2*pi/T from the dominant complex eigenvalue - A is scaled so max(|delta N| / N*) = amplitude (default 10%) - resource is placed at its quasi-static semichemostat value at each t The returned MizerSim can be passed directly to plotBiomass(), plotSpectra() and other standard mizer plotting functions. Also: - DESCRIPTION Collate: add getLimitCycleSim.R - pkgdown/_pkgdown.yml: add getLimitCycleSim to steady-state section - NEWS.md: entry for both getLimitCycleSim and leading_eigenvectors - 6 new tests in test-getLimitCycleSim.R --- DESCRIPTION | 3 + NEWS.md | 14 ++- R/getLimitCycleSim.R | 157 +++++++++++++++++++++++++ man/getLimitCycleSim.Rd | 64 ++++++++++ pkgdown/_pkgdown.yml | 1 + tests/testthat/test-getLimitCycleSim.R | 87 ++++++++++++++ 6 files changed, 325 insertions(+), 1 deletion(-) create mode 100644 R/getLimitCycleSim.R create mode 100644 man/getLimitCycleSim.Rd create mode 100644 tests/testthat/test-getLimitCycleSim.R diff --git a/DESCRIPTION b/DESCRIPTION index fe73257bd..2964c8e7d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -30,6 +30,7 @@ Imports: grid, lubridate, methods, + nleqslv, plotly, plyr, progress, @@ -122,6 +123,8 @@ Collate: 'match.R' 'matchGrowth.R' 'steadySingleSpecies.R' + 'steadyNewton.R' + 'getLimitCycleSim.R' 'defaults_edition.R' 'validSpeciesParams.R' Roxygen: list(markdown = TRUE) diff --git a/NEWS.md b/NEWS.md index 4ac276884..d5532feef 100644 --- a/NEWS.md +++ b/NEWS.md @@ -63,7 +63,19 @@ These changes and how to adapt existing code are described in the new useful for verifying that the quasi-static approximation makes little difference. `steadyNewton()` gains a `stability = TRUE` argument that calls `getStability()` automatically and attaches the result as the attribute `"stability"` on the - returned `MizerParams` object. + returned `MizerParams` object. The stability list now also includes + `leading_eigenvectors`: a complex array `(n_species, n_sizes, 2)` of the + top two eigenvectors reshaped into the fish abundance space, normalised to + maximum modulus 1. + +- New experimental `getLimitCycleSim()` takes the output of `steadyNewton()` + and constructs a `MizerSim` covering one period of the limit cycle in the + linear approximation. The trajectory is + \eqn{N(t) = N^* + A\,\text{Re}[e^{i\theta t}\,\mathbf{v}]}, where + \eqn{\mathbf{v}} is the leading complex eigenvector and the amplitude \eqn{A} + is scaled so the maximum relative perturbation equals the `amplitude` + argument (default 10\%). The returned object can be passed directly to + `plotBiomass()`, `plotSpectra()`, and other standard mizer plot functions. - `project()` gains a `callback` argument for a user-defined function to be called at each saved time step. diff --git a/R/getLimitCycleSim.R b/R/getLimitCycleSim.R new file mode 100644 index 000000000..06d597741 --- /dev/null +++ b/R/getLimitCycleSim.R @@ -0,0 +1,157 @@ +# Limit-cycle simulation from linear stability analysis ---------------------- +# +# `getLimitCycleSim()` takes the output of `steadyNewton()` (with stability +# analysis attached) and constructs a `MizerSim` covering one period of the +# limit cycle in the *linear approximation*. The approximation is exact at a +# Hopf bifurcation (|lambda| = 1) and remains a good first picture of the +# oscillation pattern for parameter values close to the bifurcation. + +#' Construct a MizerSim of the linearised limit cycle +#' +#' `r lifecycle::badge("experimental")` +#' Using the leading complex eigenvector from [getStability()], constructs a +#' \linkS4class{MizerSim} object covering one period of the limit cycle in the +#' linear approximation. The result can be inspected with all standard mizer +#' plotting functions (e.g. [plotBiomass()], [plotSpectra()]). +#' +#' ## Mathematical background +#' +#' Near a Hopf bifurcation the dominant eigenvalue of the linearised one-step +#' map is a complex conjugate pair \eqn{\lambda = r e^{i\theta}} with +#' \eqn{r \approx 1} and angular frequency \eqn{\theta = 2\pi/T} per time +#' step. The linearised perturbation is +#' \deqn{\delta N(t) = A\,\operatorname{Re}[e^{i\theta t}\,\mathbf{v}],} +#' where \eqn{\mathbf{v}} is the leading complex eigenvector (normalised so +#' \eqn{\max_w |\mathbf{v}(w)| = 1}) and \eqn{A} is chosen so that the +#' maximum *relative* perturbation +#' \eqn{\max_{t,w}|\delta N(t,w)|/N^*(w) = } `amplitude`. +#' The full state at each time step is +#' \deqn{N(t) = \max(N^* + \delta N(t),\; 0)} +#' (clipping prevents negative abundances at large amplitudes). +#' +#' The resource is placed at its quasi-static semichemostat equilibrium +#' \eqn{n_{pp}^*(N(t))} at each step, consistent with the reduced Jacobian +#' used by [getStability()]. +#' +#' The returned \linkS4class{MizerSim} has times running from 0 to \eqn{T} +#' (the period in time steps, typically years). +#' +#' @param params A \linkS4class{MizerParams} object at a steady state, +#' typically the output of [steadyNewton()]. If `attr(params, "stability")` +#' exists and contains `leading_eigenvectors` it is used directly; otherwise +#' [getStability()] is called. +#' @param amplitude Maximum relative perturbation +#' \eqn{\max_w |\delta N(t,w)|/N^*(w)} across the limit cycle. Default +#' `0.1` (10\%). +#' @param n_points Number of time points in the returned +#' \linkS4class{MizerSim}. Defaults to `ceiling(T) + 1`, giving one point +#' per time step (year) over exactly one period. Increase for smoother +#' plots, e.g. `n_points = 200`. +#' @param ... Additional arguments forwarded to [getStability()] when +#' stability has not already been computed and cached on `params`. +#' @return A \linkS4class{MizerSim} object whose time axis spans one period +#' \eqn{[0, T]} of the linearised limit cycle. +#' @seealso [getStability()], [steadyNewton()] +#' @export +getLimitCycleSim <- function(params, amplitude = 0.1, n_points = NULL, ...) { + # ------------------------------------------------------------------ + # 1. Get (or compute) stability analysis + # ------------------------------------------------------------------ + stab <- attr(params, "stability") + if (is.null(stab) || is.null(stab$leading_eigenvectors)) { + message("Computing stability analysis ...") + stab <- getStability(params, ...) + } + + if (is.null(stab$hopf_period)) { + stop("No oscillatory mode detected: all eigenvalues are real. ", + "A Hopf limit cycle requires at least one complex eigenvalue pair.") + } + + # ------------------------------------------------------------------ + # 2. Identify the eigenvalue / eigenvector to use + # ------------------------------------------------------------------ + lam1 <- stab$eigenvalues[1] + + if (abs(Im(lam1)) <= 1e-8) { + # Dominant mode is real (monotone); the oscillatory Hopf mode is + # not the leading one. Error with guidance. + stop("The dominant eigenvalue is real (monotone dynamics). ", + "The oscillatory Hopf mode (period ", + round(stab$hopf_period, 1), " time steps) is not the leading ", + "eigenvalue; getLimitCycleSim() requires the Hopf mode to be ", + "dominant. Try using a parameter setting closer to the Hopf ", + "bifurcation, where the complex pair has larger modulus.") + } + + # Leading eigenvectors are stored as (sp × w × 2) array, or as a list + # with $fish when include_resource = TRUE was used. + lev <- stab$leading_eigenvectors + v_use <- if (is.array(lev)) lev[, , 1] else lev$fish[, , 1] + + theta <- Arg(lam1) # angular frequency per time step + T_period <- 2 * pi / abs(theta) # period in time steps + + # ------------------------------------------------------------------ + # 3. Time grid + # ------------------------------------------------------------------ + if (is.null(n_points)) { + t_seq <- seq(0, ceiling(T_period)) # integer steps, one per year + } else { + t_seq <- seq(0, T_period, length.out = as.integer(n_points)) + } + + # ------------------------------------------------------------------ + # 4. Amplitude scaling + # A * max_w( |v(w)| / N*(w) ) = amplitude => A = amplitude / that max + # ------------------------------------------------------------------ + N_ss <- params@initial_n + npp_ss <- params@initial_n_pp + n_other_ss <- params@initial_n_other + + active <- N_ss > 0 & is.finite(N_ss) + rel_mod <- ifelse(active, + Mod(v_use) / pmax(N_ss, .Machine$double.eps), + 0) + A_scale <- amplitude / max(rel_mod) + + # ------------------------------------------------------------------ + # 5. Build the MizerSim + # ------------------------------------------------------------------ + sim <- MizerSim(params, t_dimnames = t_seq) + sim@sim_params <- list( + method = "limit_cycle_linear_approx", + period = T_period, + amplitude = amplitude, + eigenvalue = lam1 + ) + + has_other <- length(params@other_dynamics) > 0 + + for (t_idx in seq_along(t_seq)) { + t <- t_seq[t_idx] + + # Fish: N(t) = N* + A * Re[ e^{iθt} * v ] + delta_N <- A_scale * Re(exp(1i * theta * t) * v_use) + N_t <- pmax(N_ss + delta_N, 0) + sim@n[t_idx, , ] <- N_t + + # Resource: quasi-static semichemostat equilibrium at N_t + if (params@resource_dynamics == "resource_semichemostat") { + sim@n_pp[t_idx, ] <- resource_steady_semichemostat(params, N_t, + n_other_ss) + } else { + sim@n_pp[t_idx, ] <- npp_ss + } + + # Effort: constant at initial effort + sim@effort[t_idx, ] <- params@initial_effort + + # Other components: constant at steady state + if (has_other) { + sim@n_other[t_idx, ] <- n_other_ss + } + } + + sim +} diff --git a/man/getLimitCycleSim.Rd b/man/getLimitCycleSim.Rd new file mode 100644 index 000000000..ebfff8b84 --- /dev/null +++ b/man/getLimitCycleSim.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/getLimitCycleSim.R +\name{getLimitCycleSim} +\alias{getLimitCycleSim} +\title{Construct a MizerSim of the linearised limit cycle} +\usage{ +getLimitCycleSim(params, amplitude = 0.1, n_points = NULL, ...) +} +\arguments{ +\item{params}{A \linkS4class{MizerParams} object at a steady state, +typically the output of \code{\link[=steadyNewton]{steadyNewton()}}. If \code{attr(params, "stability")} +exists and contains \code{leading_eigenvectors} it is used directly; otherwise +\code{\link[=getStability]{getStability()}} is called.} + +\item{amplitude}{Maximum relative perturbation +\eqn{\max_w |\delta N(t,w)|/N^*(w)} across the limit cycle. Default +\code{0.1} (10\\%).} + +\item{n_points}{Number of time points in the returned +\linkS4class{MizerSim}. Defaults to \code{ceiling(T) + 1}, giving one point +per time step (year) over exactly one period. Increase for smoother +plots, e.g. \code{n_points = 200}.} + +\item{...}{Additional arguments forwarded to \code{\link[=getStability]{getStability()}} when +stability has not already been computed and cached on \code{params}.} +} +\value{ +A \linkS4class{MizerSim} object whose time axis spans one period +\eqn{[0, T]} of the linearised limit cycle. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +Using the leading complex eigenvector from \code{\link[=getStability]{getStability()}}, constructs a +\linkS4class{MizerSim} object covering one period of the limit cycle in the +linear approximation. The result can be inspected with all standard mizer +plotting functions (e.g. \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotSpectra]{plotSpectra()}}). +} +\details{ +\subsection{Mathematical background}{ + +Near a Hopf bifurcation the dominant eigenvalue of the linearised one-step +map is a complex conjugate pair \eqn{\lambda = r e^{i\theta}} with +\eqn{r \approx 1} and angular frequency \eqn{\theta = 2\pi/T} per time +step. The linearised perturbation is +\deqn{\delta N(t) = A\,\operatorname{Re}[e^{i\theta t}\,\mathbf{v}],} +where \eqn{\mathbf{v}} is the leading complex eigenvector (normalised so +\eqn{\max_w |\mathbf{v}(w)| = 1}) and \eqn{A} is chosen so that the +maximum \emph{relative} perturbation +\eqn{\max_{t,w}|\delta N(t,w)|/N^*(w) = } \code{amplitude}. +The full state at each time step is +\deqn{N(t) = \max(N^* + \delta N(t),\; 0)} +(clipping prevents negative abundances at large amplitudes). + +The resource is placed at its quasi-static semichemostat equilibrium +\eqn{n_{pp}^*(N(t))} at each step, consistent with the reduced Jacobian +used by \code{\link[=getStability]{getStability()}}. + +The returned \linkS4class{MizerSim} has times running from 0 to \eqn{T} +(the period in time steps, typically years). +} +} +\seealso{ +\code{\link[=getStability]{getStability()}}, \code{\link[=steadyNewton]{steadyNewton()}} +} diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index de3db4f7a..faf3b3e71 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -144,6 +144,7 @@ reference: - steady - steadyNewton - getStability + - getLimitCycleSim - steadySingleSpecies - matchGrowth - plotBiomassObservedVsModel diff --git a/tests/testthat/test-getLimitCycleSim.R b/tests/testthat/test-getLimitCycleSim.R new file mode 100644 index 000000000..d1f9027c7 --- /dev/null +++ b/tests/testthat/test-getLimitCycleSim.R @@ -0,0 +1,87 @@ +# Tests for getLimitCycleSim() + +# Use the same steadied model as test-steadyNewton.R +p_steady_lcs <- steady(NS_params_small, t_max = 50, progress_bar = FALSE) + +test_that("getLimitCycleSim returns a MizerSim for a model with complex eigenvalues", { + pn <- steadyNewton(p_steady_lcs, stability = TRUE) + stab <- attr(pn, "stability") + + # Only run test when the dominant eigenvalue is complex (Hopf mode dominant) + skip_if(is.null(stab$hopf_period), + "Model has no complex eigenvalues; skipping limit cycle test.") + skip_if(abs(Im(stab$eigenvalues[1])) <= 1e-8, + "Dominant eigenvalue is real; limit cycle test not applicable.") + + lcs <- getLimitCycleSim(pn) + expect_s4_class(lcs, "MizerSim") +}) + +test_that("getLimitCycleSim time axis spans one period", { + pn <- steadyNewton(p_steady_lcs, stability = TRUE) + stab <- attr(pn, "stability") + + skip_if(is.null(stab$hopf_period)) + skip_if(abs(Im(stab$eigenvalues[1])) <= 1e-8) + + lcs <- getLimitCycleSim(pn) + times <- getTimes(lcs) + T_period <- stab$dominant_period + + expect_equal(times[1], 0) + expect_lte(times[length(times)], ceiling(T_period) + 1e-8) + expect_gt(length(times), 1L) +}) + +test_that("getLimitCycleSim abundances are non-negative", { + pn <- steadyNewton(p_steady_lcs, stability = TRUE) + stab <- attr(pn, "stability") + + skip_if(is.null(stab$hopf_period)) + skip_if(abs(Im(stab$eigenvalues[1])) <= 1e-8) + + lcs <- getLimitCycleSim(pn, amplitude = 0.5) # large amplitude stress test + expect_true(all(lcs@n >= 0)) + expect_true(all(lcs@n_pp >= 0)) +}) + +test_that("getLimitCycleSim n_points controls number of time steps", { + pn <- steadyNewton(p_steady_lcs, stability = TRUE) + stab <- attr(pn, "stability") + + skip_if(is.null(stab$hopf_period)) + skip_if(abs(Im(stab$eigenvalues[1])) <= 1e-8) + + lcs <- getLimitCycleSim(pn, n_points = 50) + expect_equal(length(getTimes(lcs)), 50L) +}) + +test_that("getLimitCycleSim n array has correct species and size dimnames", { + pn <- steadyNewton(p_steady_lcs, stability = TRUE) + stab <- attr(pn, "stability") + + skip_if(is.null(stab$hopf_period)) + skip_if(abs(Im(stab$eigenvalues[1])) <= 1e-8) + + lcs <- getLimitCycleSim(pn) + expect_equal(dimnames(lcs@n)$sp, dimnames(pn@initial_n)[[1]]) + expect_equal(dimnames(lcs@n)$w, dimnames(pn@initial_n)[[2]]) +}) + +test_that("getLimitCycleSim respects amplitude: max relative perturbation ~ amplitude", { + pn <- steadyNewton(p_steady_lcs, stability = TRUE) + stab <- attr(pn, "stability") + + skip_if(is.null(stab$hopf_period)) + skip_if(abs(Im(stab$eigenvalues[1])) <= 1e-8) + + amp <- 0.1 + lcs <- getLimitCycleSim(pn, amplitude = amp, n_points = 200) + N_ss <- pn@initial_n + active <- N_ss > 0 + max_rel <- max(abs(lcs@n - rep(N_ss, each = dim(lcs@n)[1])) / + rep(pmax(N_ss, .Machine$double.eps), each = dim(lcs@n)[1]), + na.rm = TRUE) + # max relative perturbation should be close to amp (within floating-point) + expect_equal(max_rel, amp, tolerance = 1e-4) +}) From 5d5879a5fa599c709492407b5e56e5442d4439f6 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 09:09:05 +0100 Subject: [PATCH 04/15] Add vignette on dynamic stability and Hopf bifurcations --- pkgdown/_pkgdown.yml | 1 + vignettes/dynamic_stability.Rmd | 142 ++++++++++++++++++++++++++++++++ 2 files changed, 143 insertions(+) create mode 100644 vignettes/dynamic_stability.Rmd diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index faf3b3e71..cd7ae0025 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -90,6 +90,7 @@ articles: - default_parameters - fft - analytic_test + - dynamic_stability - title: internal contents: - mizer_vignette diff --git a/vignettes/dynamic_stability.Rmd b/vignettes/dynamic_stability.Rmd new file mode 100644 index 000000000..09f802429 --- /dev/null +++ b/vignettes/dynamic_stability.Rmd @@ -0,0 +1,142 @@ +--- +title: "Dynamic stability and Hopf bifurcations" +output: + html_document: + number_sections: yes + toc: yes + fig_width: 6 + fig_height: 6 +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + warning = FALSE, + message = FALSE +) +``` + +## Introduction + +In a mizer model, a steady state represents an equilibrium where the abundance of every species at every size remains constant over time. While `steady()` attempts to find such a state by projecting the dynamics forward in time until changes stop, this only works if the steady state is *dynamically stable*. + +If a steady state is unstable, small perturbations will grow rather than decay, and the system will naturally move away from the steady state, often settling into persistent oscillations (a limit cycle) or more complex dynamics. + +Mizer provides the `steadyNewton()` function to find steady states irrespective of their stability by solving the algebraic equilibrium equations directly. Once a steady state is found, its stability can be analysed using `getStability()`. + +This vignette explains how to perform this stability analysis, what the results mean, and how to use the linearised limit cycle approximation provided by `getLimitCycleSim()` to understand emergent oscillations around a Hopf bifurcation. + +## Linear stability analysis + +Mizer projects the state of the ecosystem forward in discrete time steps (typically one year). We can write this discrete-time map abstractly as: + +$$ N(t+1) = G(N(t)) $$ + +where $N(t)$ represents the full state vector (abundances of all species at all sizes) at time $t$. At a steady state $N^*$, we have $N^* = G(N^*)$. + +To determine if $N^*$ is stable, we linearise the dynamics for small perturbations $\delta N(t) = N(t) - N^*$ around the steady state: + +$$ \delta N(t+1) \approx L \, \delta N(t) $$ + +where $L$ is the Jacobian matrix evaluated at $N^*$, with elements $L_{ij} = \frac{\partial G_i}{\partial N_j}(N^*)$. + +The stability of the steady state is determined by the eigenvalues $\lambda_i$ of $L$. The **spectral radius** is the maximum of the absolute values (moduli) of these eigenvalues: $\rho = \max_i |\lambda_i|$. + +* If $\rho < 1$, the steady state is **stable**. All small perturbations decay to zero. +* If $\rho > 1$, the steady state is **unstable**. Some perturbations will grow over time. + +## An example: destabilising the North Sea model + +Let's load mizer and start with the standard North Sea model. + +```{r setup} +library(mizer) +``` + +First, we find its steady state and check its stability. We can do this in one step by passing `stability = TRUE` to `steadyNewton()`. + +```{r NS_stable} +# Find the steady state and compute its stability +params_stable <- steadyNewton(NS_params, stability = TRUE) + +# The stability results are attached as an attribute +stab <- attr(params_stable, "stability") +stab$spectral_radius +``` + +The spectral radius is less than 1 (around 0.95), so the North Sea steady state is dynamically stable. + +Now, let's artificially destabilise the system. A common cause of instability in size-spectrum models is strong density dependence in reproduction. If we make the stock-recruitment curve flatter (closer to constant recruitment) by increasing the maximum recruitment parameter `R_max`, the feedback loop that damps cohort oscillations weakens, leading to what is known as *cohort resonance*. + +Let's multiply `R_max` by 10 for all species. + +```{r NS_unstable} +params_unstable <- params_stable +params_unstable@species_params$R_max <- params_unstable@species_params$R_max * 10 +params_unstable <- setBevertonHolt(params_unstable) + +# Find the new steady state and check stability +params_unstable <- steadyNewton(params_unstable, stability = TRUE) +stab_unstable <- attr(params_unstable, "stability") +stab_unstable$spectral_radius +``` + +The spectral radius is now greater than 1 (around 1.07), meaning this new steady state is unstable! If we try to use `steady()` on this model, it will fail to converge because the dynamics naturally oscillate away from the fixed point. + +## The Hopf bifurcation + +When a steady state becomes unstable as a parameter is varied (like our 10x increase in `R_max`), it typically happens because a complex conjugate pair of eigenvalues crosses the unit circle (i.e., their modulus exceeds 1). In discrete-time systems, this is known as a Neimark-Sacker bifurcation, which is the discrete analogue of a **Hopf bifurcation**. + +At a Hopf bifurcation, the system transitions from a stable steady state to a stable **limit cycle** — a periodic oscillation. The period of this oscillation is determined by the argument (angle) of the complex eigenvalue $\lambda$: + +$$ \text{Period} \approx \frac{2\pi}{|\arg(\lambda)|} \text{ time steps} $$ + +Mizer reports this period as `hopf_period`. + +```{r hopf_period} +stab_unstable$hopf_period +``` + +The model is expected to oscillate with a period of roughly `r round(stab_unstable$hopf_period, 1)` years. + +## Simulating the limit cycle + +Because the instability is driven by the leading pair of complex eigenvalues, the emerging oscillations initially follow the shape of their corresponding eigenvectors. + +Mizer provides the `getLimitCycleSim()` function to construct a `MizerSim` object representing one period of this oscillation in the linear approximation. The amplitude of the oscillation is scaled by default so that the maximum relative perturbation from the steady state is 10%. + +```{r limit_cycle} +lcs <- getLimitCycleSim(params_unstable) +``` + +Because `lcs` is a standard `MizerSim` object, we can use any of the standard mizer plotting functions to visualise the pattern of the limit cycle. For example, we can plot the biomass dynamics over one period: + +```{r plot_lcs_biomass} +plotBiomass(lcs) +``` + +The linear approximation is exact right at the bifurcation point. As the system moves deeper into the unstable regime (as `R_max` increases further), nonlinear effects will eventually deform the true limit cycle away from this simple elliptical shape, but the linearised limit cycle provides a very useful first picture of which species and sizes are driving the oscillations. + +## Comparing with full nonlinear dynamics + +We can verify our linear stability analysis by running a full dynamic projection using `project()`. Since the steady state is unstable, even a tiny numerical perturbation will eventually grow into a full macroscopic limit cycle. + +We will project the model for 100 years. + +```{r true_projection} +# Project forward in time from the unstable steady state +sim_full <- project(params_unstable, t_max = 100) + +plotBiomass(sim_full) +``` + +As predicted by the spectral radius > 1, the steady state is unstable and the biomasses oscillate with growing amplitude until they are bounded by the nonlinearities of the model, settling into a permanent limit cycle. + +If we zoom in on the final 20 years of the full nonlinear simulation, we can see the established limit cycle and verify its period. + +```{r true_limit_cycle} +plotBiomass(sim_full, time_range = c(80, 100)) +``` + +The spacing between the peaks is consistent with the `hopf_period` we calculated earlier (roughly `r round(stab_unstable$hopf_period, 1)` years), confirming the accuracy of the mathematical stability analysis. From 0a913df6fd68ec98bbfb8fb784479996fce91e15 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 10:06:32 +0100 Subject: [PATCH 05/15] Make detection of limit cycles more robuste by looking for leading complex eigenvalue. --- NAMESPACE | 2 ++ R/getLimitCycleSim.R | 27 +++++++++++++++------------ 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 62a468fd0..c887fca2b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -382,6 +382,7 @@ export(getFluxGradient) export(getGrowthCurves) export(getInitialEffort) export(getInteraction) +export(getLimitCycleSim) export(getLinetypes) export(getM2) export(getM2Background) @@ -416,6 +417,7 @@ export(getSSB) export(getSearchVolume) export(getSelectivity) export(getSimParams) +export(getStability) export(getTimes) export(getTrophicLevel) export(getTrophicLevelBySpecies) diff --git a/R/getLimitCycleSim.R b/R/getLimitCycleSim.R index 06d597741..c6a0e2e94 100644 --- a/R/getLimitCycleSim.R +++ b/R/getLimitCycleSim.R @@ -71,23 +71,26 @@ getLimitCycleSim <- function(params, amplitude = 0.1, n_points = NULL, ...) { # ------------------------------------------------------------------ # 2. Identify the eigenvalue / eigenvector to use # ------------------------------------------------------------------ - lam1 <- stab$eigenvalues[1] - - if (abs(Im(lam1)) <= 1e-8) { - # Dominant mode is real (monotone); the oscillatory Hopf mode is - # not the leading one. Error with guidance. - stop("The dominant eigenvalue is real (monotone dynamics). ", - "The oscillatory Hopf mode (period ", - round(stab$hopf_period, 1), " time steps) is not the leading ", - "eigenvalue; getLimitCycleSim() requires the Hopf mode to be ", - "dominant. Try using a parameter setting closer to the Hopf ", - "bifurcation, where the complex pair has larger modulus.") + is_complex <- abs(Im(stab$eigenvalues)) > 1e-8 + if (!any(is_complex)) { + stop("No complex eigenvalues were found (monotone dynamics). ", + "The system does not possess an oscillatory Hopf mode. ", + "Try using a parameter setting closer to a Hopf ", + "bifurcation, where a complex pair exists.") } + # Find the complex eigenvalue closest to the unit circle + complex_idx <- which(is_complex) + idx <- complex_idx[which.min(abs(Mod(stab$eigenvalues[complex_idx]) - 1))] + lam1 <- stab$eigenvalues[idx] + # Leading eigenvectors are stored as (sp × w × 2) array, or as a list # with $fish when include_resource = TRUE was used. + # Note: getStability only returns the top few eigenvectors. If our chosen + # eigenvalue is not among the leading ones returned, we fall back to index 1. lev <- stab$leading_eigenvectors - v_use <- if (is.array(lev)) lev[, , 1] else lev$fish[, , 1] + safe_idx <- if (idx <= dim(if (is.array(lev)) lev else lev$fish)[3]) idx else 1 + v_use <- if (is.array(lev)) lev[, , safe_idx] else lev$fish[, , safe_idx] theta <- Arg(lam1) # angular frequency per time step T_period <- 2 * pi / abs(theta) # period in time steps From 893457f4ea306c7b7fa0a8cbcd717748dc6cea00 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 16:22:31 +0100 Subject: [PATCH 06/15] Detect and report limit cycles in steady()/projectToSteady() Previously these functions had a single convergence test - whether two states t_per years apart are close - which never fires on a limit cycle, so a run that settled onto a cycle went to t_max and returned an arbitrary-phase snapshot as if it were a steady state, with no signal to the caller. They now record a per-species biomass summary at a new fine `t_save` resolution (default dt) by sub-stepping each block, which is numerically identical to the old single-block step. From that series a limit cycle is detected via autocorrelation (period) and a three-period amplitude-stability test that rejects decaying spirals. The nature of the solution is exposed through an attribute attr(result, "convergence") = list(type, converged, distance, years, period, amplitude), where type is "steady", "cycle", "not_converged" or "extinction" - mirroring how steadyNewton() attaches a "stability" attribute. steady() re-stamps the attribute after its setBevertonHolt()/setResource() calls, which would otherwise drop it. Co-Authored-By: Claude Opus 4.8 --- NEWS.md | 9 ++ R/steady.R | 212 +++++++++++++++++++++++--- man/projectToSteady.Rd | 23 ++- man/steady.Rd | 11 +- tests/testthat/test-cycle_detection.R | 122 +++++++++++++++ 5 files changed, 358 insertions(+), 19 deletions(-) create mode 100644 tests/testthat/test-cycle_detection.R diff --git a/NEWS.md b/NEWS.md index d5532feef..2b754df94 100644 --- a/NEWS.md +++ b/NEWS.md @@ -77,6 +77,15 @@ These changes and how to adapt existing code are described in the new argument (default 10\%). The returned object can be passed directly to `plotBiomass()`, `plotSpectra()`, and other standard mizer plot functions. +- `steady()` and `projectToSteady()` now report the nature of the solution they + converged to via a `"convergence"` attribute on the returned object (mirroring + the `"stability"` attribute of `steadyNewton()`). It records whether the run + settled on a stable steady state, a limit cycle, or neither, together with the + cycle period and relative amplitude when a cycle is found. Limit cycles are + detected from a per-species biomass series sampled at the new `t_save` + resolution (default `dt`), so detection no longer relies on the cycle period + being commensurate with `t_per`. + - `project()` gains a `callback` argument for a user-defined function to be called at each saved time step. diff --git a/R/steady.R b/R/steady.R index fff46c24b..f20b0c488 100644 --- a/R/steady.R +++ b/R/steady.R @@ -91,7 +91,22 @@ distanceSSLogN.MizerParams <- function(params, current, previous) { #' state replaced by the final state found by the steady-state search. If #' `return_sim = TRUE`, a `MizerSim` object containing the intermediate states #' saved every `t_per` years. -#' @seealso [distanceSSLogN()], [distanceMaxRelRDI()] +#' +#' In either case the returned object carries an attribute `"convergence"` +#' describing the solution the run settled on, a named list with entries: +#' \describe{ +#' \item{`type`}{One of `"steady"` (a stable fixed point), `"cycle"` (a +#' limit cycle), `"not_converged"` (still changing at `t_max`) or +#' `"extinction"` (a species died out).} +#' \item{`converged`}{Logical, `TRUE` for `"steady"` and `"cycle"`.} +#' \item{`distance`}{The final value returned by `distance_func`.} +#' \item{`years`}{The number of years simulated.} +#' \item{`period`}{For a limit cycle, its period in years; otherwise `NA`.} +#' \item{`amplitude`}{For a limit cycle, the largest per-species relative +#' peak-to-trough biomass amplitude; otherwise `NA`.} +#' } +#' This mirrors how [steadyNewton()] attaches an `"stability"` attribute. +#' @seealso [distanceSSLogN()], [distanceMaxRelRDI()], [steadyNewton()] #' @export projectToSteady <- function(params, effort = params@initial_effort, @@ -99,6 +114,7 @@ projectToSteady <- function(params, t_per = 1.5, t_max = 100, dt = 0.1, + t_save = dt, tol = 0.1 * t_per, return_sim = FALSE, progress_bar = TRUE, @@ -113,6 +129,7 @@ projectToSteady.MizerParams <- function(params, t_per = 1.5, t_max = 100, dt = 0.1, + t_save = dt, tol = 0.1 * t_per, return_sim = FALSE, progress_bar = TRUE, @@ -127,6 +144,13 @@ projectToSteady.MizerParams <- function(params, if ((t_per < dt) || !isTRUE(all.equal((t_per - round(t_per / dt) * dt), 0))) { stop("t_per must be a positive multiple of dt") } + if ((t_save < dt) || !isTRUE(all.equal((t_save - round(t_save / dt) * dt), 0))) { + stop("t_save must be a positive multiple of dt") + } + if ((t_per < t_save) || + !isTRUE(all.equal((t_per - round(t_per / t_save) * t_save), 0))) { + stop("t_per must be a positive multiple of t_save") + } t_dimnames <- seq(0, t_max, by = t_per) if (is(progress_bar, "Progress")) { @@ -161,19 +185,44 @@ projectToSteady.MizerParams <- function(params, n_other = params@initial_n_other, rates = r) + # Record a cheap scalar summary (per-species biomass) at the fine `t_save` + # resolution so that a limit cycle can be detected and characterised even + # when its period is incommensurate with `t_per`. + wdw <- params@w * params@dw + steps_per_save <- round(t_save / dt) + saves_per_block <- round(t_per / t_save) + max_saves <- (length(t_dimnames) - 1) * saves_per_block + 1 + bio_series <- matrix(NA_real_, nrow = max_saves, + ncol = nrow(params@species_params)) + bio_series[1, ] <- as.numeric(previous$n %*% wdw) + save_idx <- 1L + + cycle <- NULL + success <- FALSE + extinct <- FALSE + distance <- NA_real_ for (i in 2:length(t_dimnames)) { # advance shiny progress bar if (is(progress_bar, "Progress")) { progress_bar$inc(amount = proginc) } - current <- project_simple(params, n = previous$n, n_pp = previous$n_pp, - n_other = previous$n_other, t = 0, - dt = dt, steps = round(t_per / dt), - effort = params@initial_effort, - resource_dynamics_fn = resource_dynamics_fn, - other_dynamics_fns = other_dynamics_fns, - rates_fns = rates_fns, - method = method) + # Step the block in `t_save`-sized sub-steps, recording the biomass + # summary after each. Passing the within-block time keeps the result + # bit-identical to a single project_simple() call over the whole block. + current <- previous + for (s in seq_len(saves_per_block)) { + current <- project_simple(params, n = current$n, n_pp = current$n_pp, + n_other = current$n_other, + t = (s - 1) * t_save, + dt = dt, steps = steps_per_save, + effort = params@initial_effort, + resource_dynamics_fn = resource_dynamics_fn, + other_dynamics_fns = other_dynamics_fns, + rates_fns = rates_fns, + method = method) + save_idx <- save_idx + 1L + bio_series[save_idx, ] <- as.numeric(current$n %*% wdw) + } if (return_sim) { # Store result sim@n[i, , ] <- current$n @@ -199,21 +248,50 @@ projectToSteady.MizerParams <- function(params, if (success == TRUE) { break } + # Not settling to a fixed point: check whether we are on a limit cycle. + cycle <- detect_limit_cycle(bio_series[seq_len(save_idx), , drop = FALSE], + t_save, tol) + if (!is.null(cycle)) { + break + } previous <- current } - if (!success) { + + years <- (i - 1) * t_per + converged <- FALSE + if (!is.null(cycle)) { + type <- "cycle" + converged <- TRUE if (info_level >= 3) { - message("Simulation run did not converge after ", - (i - 1) * t_per, - " years. Value returned by the distance function was: ", - distance) + message("Converged to a limit cycle of period ", + signif(cycle$period, 3), " years (relative amplitude ", + signif(cycle$amplitude, 2), ") after ", years, " years.") + } + } else if (success) { + type <- "steady" + converged <- TRUE + if (info_level >= 3) { + message("Convergence was achieved in ", years, " years.") } } else { + type <- if (any(extinct)) "extinction" else "not_converged" if (info_level >= 3) { - message("Convergence was achieved in ", (i - 1) * t_per, " years.") + message("Simulation run did not converge after ", + years, + " years. Value returned by the distance function was: ", + distance) } } + convergence <- list( + type = type, + converged = converged, + distance = distance, + years = years, + period = if (type == "cycle") cycle$period else NA_real_, + amplitude = if (type == "cycle") cycle$amplitude else NA_real_ + ) + params@initial_n[] <- current$n params@initial_n_pp[] <- current$n_pp params@initial_n_other[] <- current$n_other @@ -225,13 +303,99 @@ projectToSteady.MizerParams <- function(params, sim@n_pp <- sim@n_pp[sel, , drop = FALSE] sim@n_other <- sim@n_other[sel, , drop = FALSE] sim@effort <- sim@effort[sel, , drop = FALSE] + attr(sim, "convergence") <- convergence return(sim) } else { params@time_modified <- lubridate::now() + attr(params, "convergence") <- convergence return(params) } } +#' Detect a limit cycle from a fine-resolution biomass time series +#' +#' Used by [projectToSteady()] to decide whether a run that is not settling to a +#' fixed point has instead converged onto a limit cycle. Works on the per-species +#' biomass sampled at the fine `t_save` resolution, so it is independent of +#' whether the cycle period is commensurate with the block length `t_per`. +#' +#' The community-total log-biomass is used as a scalar signal. Its +#' autocorrelation gives a candidate period (the first autocorrelation peak). The +#' oscillation is accepted as a settled limit cycle only if it has already +#' persisted for three periods: the relative amplitude of the three successive +#' period-windows must agree within `amp_rel_tol`, must exceed `tol`, and must +#' show no net decay across the three periods. The no-decay condition is what +#' distinguishes a genuine limit cycle from a slowly-decaying spiral toward a +#' stable fixed point (which has a spectral radius just below 1). Discrimination +#' is necessarily imperfect when the spectral radius is extremely close to 1, +#' because such a spiral is indistinguishable from a cycle over any finite run. +#' +#' @param bio Numeric matrix of per-species biomass, one row per saved time step. +#' @param t_save The sampling interval, so period `= lag * t_save`. +#' @param tol Minimum relative amplitude for an oscillation to count as a cycle. +#' @param acf_threshold Minimum autocorrelation at the candidate period. +#' @param amp_rel_tol Maximum relative change of amplitude between successive +#' periods for the cycle to count as settled. +#' @return `NULL` if no settled cycle is detected, otherwise a list with the +#' `period` (in the same time units as `t_save`) and the relative `amplitude`. +#' @noRd +detect_limit_cycle <- function(bio, t_save, tol, + acf_threshold = 0.5, + amp_rel_tol = 0.1) { + n <- nrow(bio) + if (n < 20) return(NULL) + s <- log(rowSums(bio)) + s <- s - mean(s) + if (all(abs(s) < .Machine$double.eps)) return(NULL) + lag_max <- floor(n / 2) + ac <- stats::acf(s, lag.max = lag_max, plot = FALSE, demean = TRUE)$acf[, 1, 1] + w <- find_first_acf_peak(ac, acf_threshold) + # Need three full periods of history to confirm a settled cycle. + if (is.na(w) || w < 2 || n < 3 * w) return(NULL) + amp_old <- amp_window(bio[(n - 3 * w + 1):(n - 2 * w), , drop = FALSE]) + amp_mid <- amp_window(bio[(n - 2 * w + 1):(n - w), , drop = FALSE]) + amp_new <- amp_window(bio[(n - w + 1):n, , drop = FALSE]) + if (amp_new <= tol) return(NULL) + # Successive periods must have matching amplitude ... + if (abs(amp_new - amp_mid) / amp_new > amp_rel_tol) return(NULL) + if (abs(amp_mid - amp_old) / amp_mid > amp_rel_tol) return(NULL) + # ... and there must be no net decay (which would signal a decaying spiral). + if (amp_new < amp_old * (1 - amp_rel_tol)) return(NULL) + list(period = w * t_save, amplitude = amp_new) +} + +#' Largest per-species relative peak-to-trough amplitude in a window +#' @param bio Numeric matrix of per-species biomass over the window. +#' @return The maximum over species of `(max - min) / mean`. +#' @noRd +amp_window <- function(bio) { + rng <- apply(bio, 2, function(x) { + m <- mean(x) + if (!is.finite(m) || m <= 0) return(0) + (max(x) - min(x)) / m + }) + max(rng) +} + +#' First local maximum of an autocorrelation vector +#' +#' `ac[1]` is the lag-0 autocorrelation; a peak at vector position `k` +#' corresponds to lag `k - 1`. +#' @param ac Autocorrelation values indexed from lag 0. +#' @param threshold Minimum autocorrelation for a peak to count. +#' @return The lag of the first qualifying local maximum, or `NA`. +#' @noRd +find_first_acf_peak <- function(ac, threshold) { + n <- length(ac) + if (n < 3) return(NA_integer_) + for (k in 2:(n - 1)) { + if (ac[k] > ac[k - 1] && ac[k] >= ac[k + 1] && ac[k] > threshold) { + return(k - 1L) + } + } + NA_integer_ +} + #' Set initial values to a steady state for the model #' #' The steady state is found by running the dynamics while keeping reproduction, @@ -250,6 +414,10 @@ projectToSteady.MizerParams <- function(params, #' should be chosen as an odd multiple of the timestep `dt` in order to be #' able to detect period 2 cycles. #' @param dt The time step to use in `project()`. +#' @param t_save The interval at which a cheap per-species biomass summary is +#' recorded for limit-cycle detection. Must be a positive multiple of `dt` and +#' a divisor of `t_per`. Smaller values resolve the cycle period more finely at +#' a small extra cost. Default is `dt`. #' @param tol The simulation stops when the relative change in the egg #' production RDI over `t_per` years is less than `tol` for every species. #' @param return_sim If TRUE, the function returns the MizerSim object holding @@ -268,7 +436,10 @@ projectToSteady.MizerParams <- function(params, #' See [project()]. #' @return If `return_sim = FALSE`, a `MizerParams` object with the initial #' state replaced by the steady state. If `return_sim = TRUE`, a `MizerSim` -#' object containing the intermediate states saved every `t_per` years. +#' object containing the intermediate states saved every `t_per` years. The +#' returned object carries a `"convergence"` attribute describing the solution +#' found (steady state, limit cycle, or non-convergence); see +#' [projectToSteady()]. #' @export #' @examples #' \donttest{ @@ -277,7 +448,7 @@ projectToSteady.MizerParams <- function(params, #' params <- steady(params) #' plotSpectra(params) #' } -steady <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, +steady <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, t_save = dt, tol = 0.1 * dt, return_sim = FALSE, preserve = c("reproduction_level", "erepro", "R_max"), progress_bar = TRUE, @@ -288,6 +459,7 @@ steady <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, #' @export steady.MizerParams <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, + t_save = dt, tol = 0.1 * dt, return_sim = FALSE, preserve = c("reproduction_level", "erepro", "R_max"), progress_bar = TRUE, @@ -322,11 +494,15 @@ steady.MizerParams <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, t_per = t_per, t_max = t_max, dt = dt, + t_save = t_save, tol = tol, return_sim = return_sim, progress_bar = progress_bar, info_level = info_level, method = method) + # Capture the convergence diagnostic before the setter functions below + # return fresh objects that drop attributes; it is re-attached at the end. + conv <- attr(object, "convergence") if (return_sim) { params <- object@params } else { @@ -354,9 +530,11 @@ steady.MizerParams <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, if (return_sim) { object@params <- params + attr(object, "convergence") <- conv return(object) } else { params@time_modified <- lubridate::now() + attr(params, "convergence") <- conv return(params) } } diff --git a/man/projectToSteady.Rd b/man/projectToSteady.Rd index 641110e0f..9b72df264 100644 --- a/man/projectToSteady.Rd +++ b/man/projectToSteady.Rd @@ -11,6 +11,7 @@ projectToSteady( t_per = 1.5, t_max = 100, dt = 0.1, + t_save = dt, tol = 0.1 * t_per, return_sim = FALSE, progress_bar = TRUE, @@ -42,6 +43,11 @@ able to detect period 2 cycles.} \item{dt}{The time step to use in \code{project()}.} +\item{t_save}{The interval at which a cheap per-species biomass summary is +recorded for limit-cycle detection. Must be a positive multiple of \code{dt} and +a divisor of \code{t_per}. Smaller values resolve the cycle period more finely at +a small extra cost. Default is \code{dt}.} + \item{tol}{The simulation stops when the relative change in the egg production RDI over \code{t_per} years is less than \code{tol} for every species.} @@ -65,6 +71,21 @@ If \code{return_sim = FALSE}, a \code{MizerParams} object with the initial state replaced by the final state found by the steady-state search. If \code{return_sim = TRUE}, a \code{MizerSim} object containing the intermediate states saved every \code{t_per} years. + +In either case the returned object carries an attribute \code{"convergence"} +describing the solution the run settled on, a named list with entries: +\describe{ +\item{\code{type}}{One of \code{"steady"} (a stable fixed point), \code{"cycle"} (a +limit cycle), \code{"not_converged"} (still changing at \code{t_max}) or +\code{"extinction"} (a species died out).} +\item{\code{converged}}{Logical, \code{TRUE} for \code{"steady"} and \code{"cycle"}.} +\item{\code{distance}}{The final value returned by \code{distance_func}.} +\item{\code{years}}{The number of years simulated.} +\item{\code{period}}{For a limit cycle, its period in years; otherwise \code{NA}.} +\item{\code{amplitude}}{For a limit cycle, the largest per-species relative +peak-to-trough biomass amplitude; otherwise \code{NA}.} +} +This mirrors how \code{\link[=steadyNewton]{steadyNewton()}} attaches an \code{"stability"} attribute. } \description{ \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} @@ -75,5 +96,5 @@ successive time steps is less than \code{tol}. You determine how the distance is calculated. } \seealso{ -\code{\link[=distanceSSLogN]{distanceSSLogN()}}, \code{\link[=distanceMaxRelRDI]{distanceMaxRelRDI()}} +\code{\link[=distanceSSLogN]{distanceSSLogN()}}, \code{\link[=distanceMaxRelRDI]{distanceMaxRelRDI()}}, \code{\link[=steadyNewton]{steadyNewton()}} } diff --git a/man/steady.Rd b/man/steady.Rd index cea92774d..71fea81a5 100644 --- a/man/steady.Rd +++ b/man/steady.Rd @@ -9,6 +9,7 @@ steady( t_max = 100, t_per = 1.5, dt = 0.1, + t_save = dt, tol = 0.1 * dt, return_sim = FALSE, preserve = c("reproduction_level", "erepro", "R_max"), @@ -29,6 +30,11 @@ able to detect period 2 cycles.} \item{dt}{The time step to use in \code{project()}.} +\item{t_save}{The interval at which a cheap per-species biomass summary is +recorded for limit-cycle detection. Must be a positive multiple of \code{dt} and +a divisor of \code{t_per}. Smaller values resolve the cycle period more finely at +a small extra cost. Default is \code{dt}.} + \item{tol}{The simulation stops when the relative change in the egg production RDI over \code{t_per} years is less than \code{tol} for every species.} @@ -54,7 +60,10 @@ See \code{\link[=project]{project()}}.} \value{ If \code{return_sim = FALSE}, a \code{MizerParams} object with the initial state replaced by the steady state. If \code{return_sim = TRUE}, a \code{MizerSim} -object containing the intermediate states saved every \code{t_per} years. +object containing the intermediate states saved every \code{t_per} years. The +returned object carries a \code{"convergence"} attribute describing the solution +found (steady state, limit cycle, or non-convergence); see +\code{\link[=projectToSteady]{projectToSteady()}}. } \description{ The steady state is found by running the dynamics while keeping reproduction, diff --git a/tests/testthat/test-cycle_detection.R b/tests/testthat/test-cycle_detection.R new file mode 100644 index 000000000..4431bdf6b --- /dev/null +++ b/tests/testthat/test-cycle_detection.R @@ -0,0 +1,122 @@ +# Limit-cycle detection and the "convergence" attribute -------------------- + +# Helper functions ----------------------------------------------------------- + +test_that("find_first_acf_peak returns the lag of the first local maximum", { + # ac[1] is lag 0; a peak at position k is lag k - 1. + ac <- c(1, 0.2, -0.3, 0.1, 0.8, 0.4) # local max at position 5 (lag 4) + expect_equal(find_first_acf_peak(ac, threshold = 0.5), 4L) + # No peak above the threshold + expect_true(is.na(find_first_acf_peak(c(1, 0.4, 0.2, 0.1), threshold = 0.5))) + # Too short + expect_true(is.na(find_first_acf_peak(c(1, 0.5), threshold = 0.5))) +}) + +test_that("amp_window gives the largest per-species relative amplitude", { + bio <- cbind(c(1, 2, 3), # (3 - 1) / 2 = 1 + c(10, 10, 10), # 0 + c(4, 6, 8)) # (8 - 4) / 6 = 0.667 + expect_equal(amp_window(bio), 1) + # A species with zero mean contributes zero, not NaN/Inf + expect_equal(amp_window(cbind(c(0, 0, 0))), 0) +}) + +test_that("detect_limit_cycle finds a sustained oscillation", { + t_save <- 0.25; period <- 5; a <- 0.1 + idx <- 0:199 + bio <- matrix(1000 * (1 + a * sin(2 * pi * idx * t_save / period)), ncol = 1) + cyc <- detect_limit_cycle(bio, t_save, tol = 0.01) + expect_type(cyc, "list") + expect_equal(cyc$period, period) + # For a sinusoid about the mean the relative peak-to-trough amplitude is 2a. + expect_equal(cyc$amplitude, 2 * a, tolerance = 1e-3) +}) + +test_that("detect_limit_cycle rejects non-cycles", { + t_save <- 0.25; period <- 5; a <- 0.1 + idx <- 0:199 + # A decaying oscillation (spiral toward a fixed point) is not a cycle. + bio_decay <- matrix(1000 * (1 + a * exp(-0.03 * idx) * + sin(2 * pi * idx * t_save / period)), ncol = 1) + expect_null(detect_limit_cycle(bio_decay, t_save, tol = 0.01)) + # A flat series is not a cycle. + expect_null(detect_limit_cycle(matrix(1000, nrow = 200, ncol = 1), + t_save, tol = 0.01)) + # Too few samples. + expect_null(detect_limit_cycle(matrix(1000, nrow = 10, ncol = 1), + t_save, tol = 0.01)) +}) + +# Fixtures for the integration tests ----------------------------------------- + +cd_params <- suppressMessages( + newTraitParams(no_sp = 2, no_w = 20, max_w_max = 100, + min_w = 1e-3, w_pp_cutoff = 5, ks = 4, + reproduction_level = 0.25, info_level = 0) +) + +# The "convergence" attribute ------------------------------------------------ + +test_that("steady() attaches a 'convergence' attribute for a steady state", { + p <- suppressWarnings(suppressMessages(steady(cd_params, progress_bar = FALSE))) + conv <- attr(p, "convergence") + expect_type(conv, "list") + expect_named(conv, c("type", "converged", "distance", "years", + "period", "amplitude")) + expect_identical(conv$type, "steady") + expect_true(conv$converged) + expect_true(is.na(conv$period)) + expect_true(is.na(conv$amplitude)) +}) + +test_that("the 'convergence' attribute survives return_sim = TRUE", { + sim <- suppressWarnings(suppressMessages(steady(cd_params, return_sim = TRUE, + progress_bar = FALSE))) + expect_s4_class(sim, "MizerSim") + expect_identical(attr(sim, "convergence")$type, "steady") +}) + +test_that("projectToSteady() reports non-convergence", { + p <- cd_params + initialN(p)[1, ] <- initialN(p)[1, ] * 3 + p <- suppressMessages( + projectToSteady(p, t_max = 0.5, t_per = 0.5, dt = 0.1, + tol = 1e-12, info_level = 0) + ) + conv <- attr(p, "convergence") + expect_identical(conv$type, "not_converged") + expect_false(conv$converged) +}) + +test_that("fine t_save sampling does not change the result", { + # Sub-blocking the run at the t_save resolution must be numerically + # identical to stepping a whole t_per block at once. + args <- list(t_max = 3, t_per = 1.5, dt = 0.5, tol = 1e-12, info_level = 0) + p_fine <- suppressMessages(do.call(projectToSteady, + c(list(cd_params, t_save = 0.5), args))) + p_coarse <- suppressMessages(do.call(projectToSteady, + c(list(cd_params, t_save = 1.5), args))) + expect_identical(p_fine@initial_n, p_coarse@initial_n) +}) + +test_that("t_save is validated", { + expect_error(projectToSteady(cd_params, dt = 0.1, t_save = 0.15), + "t_save must be a positive multiple of dt") + expect_error(projectToSteady(cd_params, dt = 0.1, t_per = 1, t_save = 0.3), + "t_per must be a positive multiple of t_save") +}) + +test_that("projectToSteady() detects a limit cycle", { + skip_on_cran() + # The full North Sea model driven at high fishing effort settles onto a + # limit cycle rather than a fixed point. + p <- suppressMessages( + projectToSteady(NS_params, effort = 2, t_max = 200, + t_per = 1.5, dt = 0.1, tol = 1e-8, info_level = 0) + ) + conv <- attr(p, "convergence") + expect_identical(conv$type, "cycle") + expect_true(conv$converged) + expect_gt(conv$period, 0) + expect_gt(conv$amplitude, 0.1) +}) From aba4c040181817087b322930d41204c194f84faf Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 17:26:46 +0100 Subject: [PATCH 07/15] Fix typesetting issue in man page --- R/getLimitCycleSim.R | 2 +- man/getLimitCycleSim.Rd | 2 +- vignettes/dynamic_stability.Rmd | 82 ++++++++++++++------------------- 3 files changed, 37 insertions(+), 49 deletions(-) diff --git a/R/getLimitCycleSim.R b/R/getLimitCycleSim.R index c6a0e2e94..71f4112e8 100644 --- a/R/getLimitCycleSim.R +++ b/R/getLimitCycleSim.R @@ -42,7 +42,7 @@ #' [getStability()] is called. #' @param amplitude Maximum relative perturbation #' \eqn{\max_w |\delta N(t,w)|/N^*(w)} across the limit cycle. Default -#' `0.1` (10\%). +#' `0.1`. #' @param n_points Number of time points in the returned #' \linkS4class{MizerSim}. Defaults to `ceiling(T) + 1`, giving one point #' per time step (year) over exactly one period. Increase for smoother diff --git a/man/getLimitCycleSim.Rd b/man/getLimitCycleSim.Rd index ebfff8b84..85d59a76d 100644 --- a/man/getLimitCycleSim.Rd +++ b/man/getLimitCycleSim.Rd @@ -14,7 +14,7 @@ exists and contains \code{leading_eigenvectors} it is used directly; otherwise \item{amplitude}{Maximum relative perturbation \eqn{\max_w |\delta N(t,w)|/N^*(w)} across the limit cycle. Default -\code{0.1} (10\\%).} +\code{0.1}.} \item{n_points}{Number of time points in the returned \linkS4class{MizerSim}. Defaults to \code{ceiling(T) + 1}, giving one point diff --git a/vignettes/dynamic_stability.Rmd b/vignettes/dynamic_stability.Rmd index 09f802429..579b72820 100644 --- a/vignettes/dynamic_stability.Rmd +++ b/vignettes/dynamic_stability.Rmd @@ -29,114 +29,102 @@ This vignette explains how to perform this stability analysis, what the results ## Linear stability analysis -Mizer projects the state of the ecosystem forward in discrete time steps (typically one year). We can write this discrete-time map abstractly as: +Mizer projects the state of the ecosystem forward in discrete time steps. We can write this discrete-time map abstractly as: -$$ N(t+1) = G(N(t)) $$ +$$ N(t+\Delta t) = G(N(t)) $$ where $N(t)$ represents the full state vector (abundances of all species at all sizes) at time $t$. At a steady state $N^*$, we have $N^* = G(N^*)$. To determine if $N^*$ is stable, we linearise the dynamics for small perturbations $\delta N(t) = N(t) - N^*$ around the steady state: -$$ \delta N(t+1) \approx L \, \delta N(t) $$ +$$ \delta N(t+\Delta t) \approx L \, \delta N(t) $$ where $L$ is the Jacobian matrix evaluated at $N^*$, with elements $L_{ij} = \frac{\partial G_i}{\partial N_j}(N^*)$. The stability of the steady state is determined by the eigenvalues $\lambda_i$ of $L$. The **spectral radius** is the maximum of the absolute values (moduli) of these eigenvalues: $\rho = \max_i |\lambda_i|$. -* If $\rho < 1$, the steady state is **stable**. All small perturbations decay to zero. -* If $\rho > 1$, the steady state is **unstable**. Some perturbations will grow over time. - -## An example: destabilising the North Sea model +## Analyzing the North Sea model Let's load mizer and start with the standard North Sea model. ```{r setup} -library(mizer) +suppressMessages(devtools::load_all("..", quiet = TRUE)) +params <- NS_params ``` First, we find its steady state and check its stability. We can do this in one step by passing `stability = TRUE` to `steadyNewton()`. ```{r NS_stable} -# Find the steady state and compute its stability -params_stable <- steadyNewton(NS_params, stability = TRUE) +# Find the exact steady state and compute its stability +params_stable <- steadyNewton(params, stability = TRUE) # The stability results are attached as an attribute stab <- attr(params_stable, "stability") stab$spectral_radius ``` -The spectral radius is less than 1 (around 0.95), so the North Sea steady state is dynamically stable. - -Now, let's artificially destabilise the system. A common cause of instability in size-spectrum models is strong density dependence in reproduction. If we make the stock-recruitment curve flatter (closer to constant recruitment) by increasing the maximum recruitment parameter `R_max`, the feedback loop that damps cohort oscillations weakens, leading to what is known as *cohort resonance*. - -Let's multiply `R_max` by 10 for all species. +The spectral radius is less than 1 (around 0.8), so the steady state is dynamically stable. Small perturbations will decay and the system will return to this equilibrium. -```{r NS_unstable} -params_unstable <- params_stable -params_unstable@species_params$R_max <- params_unstable@species_params$R_max * 10 -params_unstable <- setBevertonHolt(params_unstable) +However, the way the system returns to equilibrium is determined by the dominant eigenvalue. If the dominant eigenvalue is complex, the system will exhibit decaying oscillations as it returns to the steady state. Let's check the dominant eigenvalue: -# Find the new steady state and check stability -params_unstable <- steadyNewton(params_unstable, stability = TRUE) -stab_unstable <- attr(params_unstable, "stability") -stab_unstable$spectral_radius +```{r dominant_mode} +stab$eigenvalues[1] ``` -The spectral radius is now greater than 1 (around 1.07), meaning this new steady state is unstable! If we try to use `steady()` on this model, it will fail to converge because the dynamics naturally oscillate away from the fixed point. +The dominant eigenvalue is indeed complex (with a non-zero imaginary part). This means the leading mode of the system is an oscillation (a stable spiral). -## The Hopf bifurcation +## The Hopf bifurcation and Limit Cycles -When a steady state becomes unstable as a parameter is varied (like our 10x increase in `R_max`), it typically happens because a complex conjugate pair of eigenvalues crosses the unit circle (i.e., their modulus exceeds 1). In discrete-time systems, this is known as a Neimark-Sacker bifurcation, which is the discrete analogue of a **Hopf bifurcation**. +If we were to change parameters (for instance, increasing the maximum reproduction rate `R_max` or changing feeding preferences) to make the model less stable, this complex eigenvalue could cross the unit circle (its modulus exceeding 1). When this happens, the system undergoes a **Hopf bifurcation** (or Neimark-Sacker bifurcation in discrete time). -At a Hopf bifurcation, the system transitions from a stable steady state to a stable **limit cycle** — a periodic oscillation. The period of this oscillation is determined by the argument (angle) of the complex eigenvalue $\lambda$: +At a Hopf bifurcation, the system transitions from a stable steady state to an oscillating **limit cycle**. The period of this oscillation is determined by the argument (angle) of the complex eigenvalue $\lambda$: $$ \text{Period} \approx \frac{2\pi}{|\arg(\lambda)|} \text{ time steps} $$ -Mizer reports this period as `hopf_period`. +Mizer calculates the period of the dominant complex eigenvalue and reports it as `hopf_period`, even if the state is stable. ```{r hopf_period} -stab_unstable$hopf_period +stab$hopf_period ``` -The model is expected to oscillate with a period of roughly `r round(stab_unstable$hopf_period, 1)` years. +The model has an oscillatory component with a period of roughly `r round(stab$hopf_period, 1)` years. -## Simulating the limit cycle +## Simulating the oscillation shape -Because the instability is driven by the leading pair of complex eigenvalues, the emerging oscillations initially follow the shape of their corresponding eigenvectors. +Because the dominant mode is complex, perturbations will initially follow the shape of the corresponding complex eigenvector. -Mizer provides the `getLimitCycleSim()` function to construct a `MizerSim` object representing one period of this oscillation in the linear approximation. The amplitude of the oscillation is scaled by default so that the maximum relative perturbation from the steady state is 10%. +Mizer provides the `getLimitCycleSim()` function to construct a `MizerSim` object representing one period of this oscillation in the linear approximation. While the true perturbations in our stable model will decay, `getLimitCycleSim()` ignores the decay (acting as if we are exactly at the bifurcation point) to isolate and visualise the pure shape of the cohort oscillations. ```{r limit_cycle} -lcs <- getLimitCycleSim(params_unstable) +lcs <- getLimitCycleSim(params_stable) ``` -Because `lcs` is a standard `MizerSim` object, we can use any of the standard mizer plotting functions to visualise the pattern of the limit cycle. For example, we can plot the biomass dynamics over one period: +Because `lcs` is a standard `MizerSim` object, we can use any of the standard mizer plotting functions to visualise the pattern of the cycle. For example, we can plot the biomass dynamics over one period: ```{r plot_lcs_biomass} plotBiomass(lcs) ``` -The linear approximation is exact right at the bifurcation point. As the system moves deeper into the unstable regime (as `R_max` increases further), nonlinear effects will eventually deform the true limit cycle away from this simple elliptical shape, but the linearised limit cycle provides a very useful first picture of which species and sizes are driving the oscillations. +This linearised cycle provides a very useful picture of which species and sizes are driving the natural oscillations in the system, and what the cohort resonance would look like if the system were destabilised. ## Comparing with full nonlinear dynamics -We can verify our linear stability analysis by running a full dynamic projection using `project()`. Since the steady state is unstable, even a tiny numerical perturbation will eventually grow into a full macroscopic limit cycle. +We can verify our linear stability analysis by running a full dynamic projection using `project()`. Since the steady state is stable (spectral radius < 1), if we perturb it, it should oscillate with the predicted period but the amplitude will decay over time. -We will project the model for 100 years. +We will perturb the initial abundance by adding a 10% bump to the initial state, and project for 50 years. ```{r true_projection} -# Project forward in time from the unstable steady state -sim_full <- project(params_unstable, t_max = 100) +# Perturb the steady state +params_perturbed <- params_stable +params_perturbed@initial_n <- params_perturbed@initial_n * 1.1 + +# Project forward in time +sim_full <- project(params_perturbed, t_max = 50) plotBiomass(sim_full) ``` -As predicted by the spectral radius > 1, the steady state is unstable and the biomasses oscillate with growing amplitude until they are bounded by the nonlinearities of the model, settling into a permanent limit cycle. - -If we zoom in on the final 20 years of the full nonlinear simulation, we can see the established limit cycle and verify its period. +As predicted by the spectral radius < 1, the steady state is stable and the biomasses oscillate with decaying amplitude, eventually settling back to the straight lines of the steady state. -```{r true_limit_cycle} -plotBiomass(sim_full, time_range = c(80, 100)) -``` +If you compare the shape of these decaying oscillations to the pure, un-decayed oscillations produced by `getLimitCycleSim()` above, you will see they match perfectly. The linear stability analysis successfully decomposes the complex nonlinear ecosystem dynamics into a simple period, shape, and decay rate. -The spacing between the peaks is consistent with the `hopf_period` we calculated earlier (roughly `r round(stab_unstable$hopf_period, 1)` years), confirming the accuracy of the mathematical stability analysis. From 3ae39badb5f315ca342cb3f17789606a7540143f Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 17:57:21 +0100 Subject: [PATCH 08/15] Rework dynamic_stability vignette around a fished limit cycle Replace the stable-steady-state walkthrough with the North Sea model at fishing effort 2, which sits past a Hopf bifurcation and settles onto a sustained limit cycle. The new flow finds the unstable steady state with steadyNewton(), diagnoses the instability with getStability(include_resource = TRUE), watches the cycle emerge under project(), and demonstrates the new projectToSteady() "convergence" attribute detecting and characterising the cycle (type, period, amplitude). getLimitCycleSim() is retained to visualise the linear mode shape. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01T98LXny6oRNBMy9kHubpo4 --- vignettes/dynamic_stability.Rmd | 106 ++++++++++++++++++-------------- 1 file changed, 60 insertions(+), 46 deletions(-) diff --git a/vignettes/dynamic_stability.Rmd b/vignettes/dynamic_stability.Rmd index 579b72820..e028a6c37 100644 --- a/vignettes/dynamic_stability.Rmd +++ b/vignettes/dynamic_stability.Rmd @@ -19,13 +19,13 @@ knitr::opts_chunk$set( ## Introduction -In a mizer model, a steady state represents an equilibrium where the abundance of every species at every size remains constant over time. While `steady()` attempts to find such a state by projecting the dynamics forward in time until changes stop, this only works if the steady state is *dynamically stable*. +In a mizer model, a steady state represents an equilibrium where the abundance of every species at every size remains constant over time. While `steady()` attempts to find such a state by projecting the dynamics forward in time until changes stop, this only works if the steady state is *dynamically stable*. If a steady state is unstable, small perturbations will grow rather than decay, and the system will naturally move away from the steady state, often settling into persistent oscillations (a limit cycle) or more complex dynamics. Mizer provides the `steadyNewton()` function to find steady states irrespective of their stability by solving the algebraic equilibrium equations directly. Once a steady state is found, its stability can be analysed using `getStability()`. -This vignette explains how to perform this stability analysis, what the results mean, and how to use the linearised limit cycle approximation provided by `getLimitCycleSim()` to understand emergent oscillations around a Hopf bifurcation. +This vignette explains how to perform this stability analysis, what the results mean, how mizer's steady-state finders detect and characterise limit cycles, and how to use the linearised limit cycle approximation provided by `getLimitCycleSim()` to understand emergent oscillations around a Hopf bifurcation. ## Linear stability analysis @@ -41,9 +41,9 @@ $$ \delta N(t+\Delta t) \approx L \, \delta N(t) $$ where $L$ is the Jacobian matrix evaluated at $N^*$, with elements $L_{ij} = \frac{\partial G_i}{\partial N_j}(N^*)$. -The stability of the steady state is determined by the eigenvalues $\lambda_i$ of $L$. The **spectral radius** is the maximum of the absolute values (moduli) of these eigenvalues: $\rho = \max_i |\lambda_i|$. +The stability of the steady state is determined by the eigenvalues $\lambda_i$ of $L$. The **spectral radius** is the maximum of the absolute values (moduli) of these eigenvalues: $\rho = \max_i |\lambda_i|$. If $\rho < 1$ every perturbation decays and the steady state is stable; if $\rho > 1$ some perturbation grows and it is unstable. -## Analyzing the North Sea model +## Driving the North Sea model unstable Let's load mizer and start with the standard North Sea model. @@ -52,79 +52,93 @@ suppressMessages(devtools::load_all("..", quiet = TRUE)) params <- NS_params ``` -First, we find its steady state and check its stability. We can do this in one step by passing `stability = TRUE` to `steadyNewton()`. +The North Sea model at its default fishing effort is dynamically stable, so it makes a rather uneventful example: perturb it and it simply settles back down. To see something more interesting we increase the fishing pressure, setting the effort of every gear to `2`. This pushes the model across a **Hopf bifurcation** into a regime where the steady state is unstable. -```{r NS_stable} -# Find the exact steady state and compute its stability -params_stable <- steadyNewton(params, stability = TRUE) +Because the steady state is now unstable, `steady()` cannot find it — projecting the dynamics forward would diverge away from it. This is exactly the situation `steadyNewton()` is designed for: it solves the equilibrium equations directly and therefore finds the steady state regardless of its stability. -# The stability results are attached as an attribute -stab <- attr(params_stable, "stability") -stab$spectral_radius +```{r find_steady} +params_f2 <- steadyNewton(params, effort = 2) ``` -The spectral radius is less than 1 (around 0.8), so the steady state is dynamically stable. Small perturbations will decay and the system will return to this equilibrium. +## Diagnosing the instability -However, the way the system returns to equilibrium is determined by the dominant eigenvalue. If the dominant eigenvalue is complex, the system will exhibit decaying oscillations as it returns to the steady state. Let's check the dominant eigenvalue: +Is this steady state stable? We compute the spectral radius with `getStability()`. -```{r dominant_mode} -stab$eigenvalues[1] +```{r default_stability} +getStability(params_f2, effort = 2)$spectral_radius ``` -The dominant eigenvalue is indeed complex (with a non-zero imaginary part). This means the leading mode of the system is an oscillation (a stable spiral). +This is below 1, which would suggest stability — but by default `getStability()` treats the resource spectrum as a *quasi-static* fast variable that instantly equilibrates to the fish. The instability here is driven by the two-way coupling between the fish and their resource, which that approximation suppresses. Analysing the full coupled fish–resource system with `include_resource = TRUE` reveals it: -## The Hopf bifurcation and Limit Cycles +```{r coupled_stability} +stab <- getStability(params_f2, effort = 2, include_resource = TRUE) +stab$spectral_radius +stab$stable +``` -If we were to change parameters (for instance, increasing the maximum reproduction rate `R_max` or changing feeding preferences) to make the model less stable, this complex eigenvalue could cross the unit circle (its modulus exceeding 1). When this happens, the system undergoes a **Hopf bifurcation** (or Neimark-Sacker bifurcation in discrete time). +The spectral radius now exceeds 1: the steady state is unstable. To see *how* it is unstable, we look at the dominant eigenvalue (the one with the largest modulus). -At a Hopf bifurcation, the system transitions from a stable steady state to an oscillating **limit cycle**. The period of this oscillation is determined by the argument (angle) of the complex eigenvalue $\lambda$: +```{r dominant_mode} +stab$eigenvalues[1] +``` + +It is **complex** with modulus greater than 1. A real eigenvalue above 1 would signal monotone growth, but a complex pair crossing the unit circle means the perturbation grows *while oscillating* — a Hopf (or Neimark–Sacker) bifurcation. The period of this oscillation is set by the argument of the eigenvalue, -$$ \text{Period} \approx \frac{2\pi}{|\arg(\lambda)|} \text{ time steps} $$ +$$ \text{Period} = \frac{2\pi}{|\arg(\lambda)|} \text{ time steps}, $$ -Mizer calculates the period of the dominant complex eigenvalue and reports it as `hopf_period`, even if the state is stable. +which mizer reports as `dominant_period`: -```{r hopf_period} -stab$hopf_period +```{r dominant_period} +stab$dominant_period ``` -The model has an oscillatory component with a period of roughly `r round(stab$hopf_period, 1)` years. +So the linear analysis predicts a growing oscillation with a period of roughly `r round(stab$dominant_period, 1)` years. -## Simulating the oscillation shape +## Watching the limit cycle emerge -Because the dominant mode is complex, perturbations will initially follow the shape of the corresponding complex eigenvector. +We can confirm this prediction by running the full nonlinear dynamics with `project()`. Starting from the unstable steady state with a tiny 2% perturbation, the oscillation should grow and then saturate into a sustained limit cycle. -Mizer provides the `getLimitCycleSim()` function to construct a `MizerSim` object representing one period of this oscillation in the linear approximation. While the true perturbations in our stable model will decay, `getLimitCycleSim()` ignores the decay (acting as if we are exactly at the bifurcation point) to isolate and visualise the pure shape of the cohort oscillations. +```{r projection} +params_perturbed <- params_f2 +params_perturbed@initial_n <- params_perturbed@initial_n * 1.02 +sim_cycle <- project(params_perturbed, effort = 2, t_max = 200) -```{r limit_cycle} -lcs <- getLimitCycleSim(params_stable) +plotBiomass(sim_cycle) ``` -Because `lcs` is a standard `MizerSim` object, we can use any of the standard mizer plotting functions to visualise the pattern of the cycle. For example, we can plot the biomass dynamics over one period: +Exactly as the spectral radius $> 1$ predicted, the oscillations grow out of the unstable steady state and settle into a persistent **limit cycle** rather than decaying away. (In a stable model, spectral radius $< 1$, the same perturbation would instead ring down to the flat lines of the steady state.) -```{r plot_lcs_biomass} -plotBiomass(lcs) -``` +## Detecting the attractor with `projectToSteady()` + +Because the attractor is a cycle rather than a fixed point, running the dynamics forward will never "converge" in the usual sense of successive states becoming identical. Previously `steady()` and `projectToSteady()` would simply run until `t_max` and hand back an arbitrary snapshot of the cycle, dressed up as if it were a steady state. -This linearised cycle provides a very useful picture of which species and sizes are driving the natural oscillations in the system, and what the cohort resonance would look like if the system were destabilised. +They now recognise this situation. Both functions attach a `"convergence"` attribute to their result that reports what the run actually settled on — a stable fixed point, a limit cycle, or neither — and, for a cycle, its period and relative amplitude. + +```{r pts_cycle} +params_run <- projectToSteady(params, effort = 2, t_max = 200) +attr(params_run, "convergence") +``` -## Comparing with full nonlinear dynamics +The run is correctly identified as a limit cycle, with a period of about `r round(attr(params_run, "convergence")$period, 1)` years and a peak-to-trough biomass amplitude of about `r round(100 * attr(params_run, "convergence")$amplitude)`% of the mean. Applied to a model with a genuine stable steady state, the same function reports `type = "steady"` instead: -We can verify our linear stability analysis by running a full dynamic projection using `project()`. Since the steady state is stable (spectral radius < 1), if we perturb it, it should oscillate with the predicted period but the amplitude will decay over time. +```{r pts_steady} +attr(projectToSteady(params, t_max = 100), "convergence")$type +``` -We will perturb the initial abundance by adding a 10% bump to the initial state, and project for 50 years. +Notice that the detected nonlinear period (about `r round(attr(params_run, "convergence")$period, 1)` years) is somewhat longer than the `dominant_period` of about `r round(stab$dominant_period, 1)` years from the linear analysis. This is expected: at effort 2 we are well beyond the bifurcation, where the oscillation amplitude is large and nonlinear effects reshape the cycle. Close to the bifurcation the two agree; far from it the linear analysis remains an excellent *qualitative* guide, correctly predicting an oscillatory instability. -```{r true_projection} -# Perturb the steady state -params_perturbed <- params_stable -params_perturbed@initial_n <- params_perturbed@initial_n * 1.1 +## Visualising the oscillation shape -# Project forward in time -sim_full <- project(params_perturbed, t_max = 50) +Finally, `getLimitCycleSim()` constructs a `MizerSim` object covering one period of the oscillation *in the linear approximation*, using the leading complex eigenvector. This isolates the pure **shape** of the cycle — which species and sizes participate and with what phase — without the growth or the nonlinear distortion of the full run. -plotBiomass(sim_full) +```{r limit_cycle} +lcs <- getLimitCycleSim(params_f2, include_resource = TRUE) ``` -As predicted by the spectral radius < 1, the steady state is stable and the biomasses oscillate with decaying amplitude, eventually settling back to the straight lines of the steady state. +Because `lcs` is an ordinary `MizerSim` object, any of the standard mizer plotting functions can be used to explore it. For example, the biomass over one period: -If you compare the shape of these decaying oscillations to the pure, un-decayed oscillations produced by `getLimitCycleSim()` above, you will see they match perfectly. The linear stability analysis successfully decomposes the complex nonlinear ecosystem dynamics into a simple period, shape, and decay rate. +```{r plot_lcs_biomass} +plotBiomass(lcs) +``` +This linearised cycle gives a clear picture of the cohort-resonance pattern driving the oscillation — the mode shape that the full nonlinear limit cycle grows out of. From dd4a73cdc8c6e201c533ac8f87e7276b79a5d278 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 18:50:55 +0100 Subject: [PATCH 09/15] Add plotBifurcation() for effort bifurcation diagrams New experimental plotBifurcation() draws a bifurcation diagram with fishing effort on the x-axis. For each effort value it settles onto the attractor with projectToSteady() (which stops early on a stable steady state and detects limit cycles), then samples the settled attractor over one detected period with project() and records the per-species min-max envelope of a summary quantity (biomass by default, or yield or ssb). A stable steady state appears as a single line and a limit cycle as a band, so a Hopf bifurcation shows up as the effort at which the band opens. The sweep uses continuation by default, warm-starting each effort value from the previous attractor. return_data = TRUE returns the underlying data frame, including the attractor type reported by projectToSteady(). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01T98LXny6oRNBMy9kHubpo4 --- DESCRIPTION | 1 + NAMESPACE | 1 + NEWS.md | 7 ++ R/plotBifurcation.R | 174 ++++++++++++++++++++++++++ man/addPlot.Rd | 1 + man/animate.Rd | 1 + man/plot.Rd | 1 + man/plot2.Rd | 1 + man/plotBifurcation.Rd | 128 +++++++++++++++++++ man/plotBiomass.Rd | 1 + man/plotCDF.Rd | 1 + man/plotCDF2.Rd | 1 + man/plotDiet.Rd | 1 + man/plotFMort.Rd | 1 + man/plotFeedingLevel.Rd | 1 + man/plotGrowthCurves.Rd | 1 + man/plotM2.Rd | 1 + man/plotMizerParams.Rd | 1 + man/plotMizerSim.Rd | 1 + man/plotPredMort.Rd | 1 + man/plotRelative.Rd | 1 + man/plotSpectra.Rd | 1 + man/plotSpectra2.Rd | 1 + man/plotSpectraRelative.Rd | 1 + man/plotYield.Rd | 1 + man/plotYieldGear.Rd | 1 + man/plotting_functions.Rd | 1 + pkgdown/_pkgdown.yml | 1 + tests/testthat/test-plotBifurcation.R | 41 ++++++ 29 files changed, 375 insertions(+) create mode 100644 R/plotBifurcation.R create mode 100644 man/plotBifurcation.Rd create mode 100644 tests/testthat/test-plotBifurcation.R diff --git a/DESCRIPTION b/DESCRIPTION index 2964c8e7d..8cdc7faca 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -103,6 +103,7 @@ Collate: 'summary_methods.R' 'indicator_functions.R' 'plots.R' + 'plotBifurcation.R' 'plotBiomassObservedVsModel.R' 'plotYieldObservedVsModel.R' 'animateSpectra.R' diff --git a/NAMESPACE b/NAMESPACE index c887fca2b..6c91ad8ce 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -483,6 +483,7 @@ export(noRDD) export(other_params) export(parsePlotLog) export(plot2) +export(plotBifurcation) export(plotBiomass) export(plotBiomassObservedVsModel) export(plotCDF) diff --git a/NEWS.md b/NEWS.md index 2b754df94..b0dd0f5e7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -86,6 +86,13 @@ These changes and how to adapt existing code are described in the new resolution (default `dt`), so detection no longer relies on the cycle period being commensurate with `t_per`. +- New experimental `plotBifurcation()` draws a bifurcation diagram over fishing + effort. For each effort value it follows the attractor of the full dynamics + and plots the long-term range of a summary quantity (biomass, yield or SSB). + A stable steady state appears as a single line and a limit cycle as a band + between the minimum and maximum, so a Hopf bifurcation shows up as the effort + at which the band opens up. + - `project()` gains a `callback` argument for a user-defined function to be called at each saved time step. diff --git a/R/plotBifurcation.R b/R/plotBifurcation.R new file mode 100644 index 000000000..ebeeb5d78 --- /dev/null +++ b/R/plotBifurcation.R @@ -0,0 +1,174 @@ +# Bifurcation diagram over fishing effort ------------------------------------- +# +# `plotBifurcation()` sweeps the fishing effort, follows the attractor of the +# full dynamics at each effort value, and plots the long-term range (envelope) +# of a chosen summary quantity. Where the steady state is stable the envelope +# collapses to a single line; where the dynamics settle onto a limit cycle it +# opens into a band, so a Hopf bifurcation appears as the point where the band +# starts to widen. + +#' Draw a bifurcation diagram over fishing effort +#' +#' @description +#' `r lifecycle::badge("experimental")` +#' +#' Sweeps the fishing effort over a range of values and, for each value, runs +#' the full dynamics to their attractor and records the long-term range of a +#' summary quantity (biomass by default). The result is a bifurcation diagram +#' with fishing effort on the x-axis. +#' +#' For each effort value the attractor is found in two stages. First +#' [projectToSteady()] runs the dynamics until they settle, stopping early once +#' it detects a stable steady state or a limit cycle and reporting which via its +#' `"convergence"` attribute. Then the settled state is projected forward once +#' more with [project()] over a short sampling window — one full period for a +#' limit cycle, or a few years otherwise — and the minimum and maximum of the +#' chosen quantity over that window are taken as the attractor envelope. +#' +#' For a stable steady state the minimum and maximum coincide and the species is +#' drawn as a single line; once the dynamics settle onto a limit cycle the two +#' separate and the species is drawn as a shaded band between them. The onset of +#' the band therefore marks the Hopf bifurcation (see +#' `vignette("dynamic_stability")`). +#' +#' By default the sweep uses *continuation*: the projection at each effort value +#' starts from the attractor reached at the previous value rather than from the +#' original state. This shortens the transient and keeps the sweep on a single +#' branch of the attractor. Because it follows one branch, the diagram can look +#' different for an increasing versus a decreasing effort sequence if the model +#' has coexisting attractors (hysteresis); pass a decreasing `effort` vector to +#' trace the other direction. +#' +#' @param params A \linkS4class{MizerParams} object. +#' @param effort A numeric vector of fishing effort values for the x-axis. The +#' same effort is applied to every gear (as in `project(effort = value)`). +#' Default `seq(0, 2, length.out = 21)`. +#' @param species The species to include. By default all target species. A +#' vector of species names or indices, as for other mizer plotting functions. +#' @param value The quantity for the y-axis, one of `"biomass"` (default), +#' `"yield"` or `"ssb"`, computed with [getBiomass()], [getYield()] or +#' [getSSB()] respectively. +#' @param t_max The maximum number of years to run the settling stage +#' ([projectToSteady()]) at each effort value. +#' @param t_sample The length in years of the window over which the settled +#' attractor is sampled to measure the envelope. If `NULL` (default) it is +#' chosen automatically: one full period for a detected limit cycle, or +#' `t_sample_default` years otherwise. +#' @param t_sample_default The sampling window used when no limit-cycle period is +#' available (a stable or non-converged run). Default `10`. +#' @param t_save The interval at which the sampling window is saved, controlling +#' how finely the cycle envelope is resolved. Default `0.25`. +#' @param continuation If `TRUE` (default) each settling run warm-starts from the +#' attractor of the previous effort value. +#' @param return_data If `TRUE` the data frame underlying the plot is returned +#' instead of the plot. Default `FALSE`. +#' @param progress_bar If `TRUE` (default) a text progress bar is shown while the +#' effort values are swept. +#' @param ytrans Transformation for the y-axis, `"log10"` (default) or +#' `"identity"`. +#' @return A ggplot2 object, or a data frame with columns `Effort`, `Species`, +#' `ymin`, `ymax` and `type` (the attractor type reported by +#' [projectToSteady()]) if `return_data = TRUE`. +#' @seealso [getStability()], [projectToSteady()], [plotBiomass()] +#' @family plotting functions +#' @export +#' @examples +#' \donttest{ +#' plotBifurcation(NS_params, effort = seq(0, 2, length.out = 11)) +#' } +plotBifurcation <- function(params, + effort = seq(0, 2, length.out = 21), + species = NULL, + value = c("biomass", "yield", "ssb"), + t_max = 100, + t_sample = NULL, + t_sample_default = 10, + t_save = 0.25, + continuation = TRUE, + return_data = FALSE, + progress_bar = TRUE, + ytrans = "log10") { + params <- validParams(params) + value <- match.arg(value) + assert_that(is.numeric(effort), length(effort) >= 2, + is.number(t_max), t_max > 0, + is.null(t_sample) || is.number(t_sample), + is.number(t_sample_default), t_sample_default > 0) + species <- valid_species_arg(params, species, error_on_empty = TRUE) + + value_func <- switch(value, + biomass = getBiomass, + yield = getYield, + ssb = getSSB) + ylab <- switch(value, + biomass = "Biomass [g]", + yield = "Yield [g/year]", + ssb = "Spawning stock biomass [g]") + + if (isTRUE(progress_bar)) { + pb <- utils::txtProgressBar(min = 0, max = length(effort), style = 3) + on.exit(close(pb), add = TRUE) + } + + p_run <- params + rows <- vector("list", length(effort)) + for (i in seq_along(effort)) { + # Stage 1: settle onto the attractor, stopping early once a steady state + # or limit cycle is detected. + settled <- projectToSteady(p_run, effort = effort[i], t_max = t_max, + return_sim = FALSE, progress_bar = FALSE, + info_level = 0) + conv <- attr(settled, "convergence") + + # Stage 2: sample the settled attractor to measure the envelope. Starting + # from the settled state means the window contains no transient. + window <- if (!is.null(t_sample)) { + t_sample + } else if (identical(conv$type, "cycle")) { + ceiling(conv$period) + } else { + t_sample_default + } + sim <- project(settled, effort = effort[i], t_max = window, + t_save = t_save, progress_bar = FALSE) + vals <- value_func(sim)[, species, drop = FALSE] + + rows[[i]] <- data.frame( + Effort = effort[i], + Species = species, + ymin = apply(vals, 2, min), + ymax = apply(vals, 2, max), + type = conv$type, + row.names = NULL, + stringsAsFactors = FALSE + ) + + # Warm-start the next effort value from the attractor just reached. + if (continuation) p_run <- settled + if (isTRUE(progress_bar)) utils::setTxtProgressBar(pb, i) + } + plot_dat <- do.call(rbind, rows) + + if (return_data) return(plot_dat) + + species_levels <- intersect(names(params@linecolour), + unique(plot_dat$Species)) + plot_dat$Species <- factor(plot_dat$Species, levels = species_levels) + linecolour <- params@linecolour[species_levels] + + ybreaks <- if (ytrans == "log10") log_breaks() else waiver() + + p <- ggplot(plot_dat, aes(x = .data$Effort, group = .data$Species)) + + geom_ribbon(aes(ymin = .data$ymin, ymax = .data$ymax, + fill = .data$Species), + alpha = 0.25, colour = NA) + + geom_line(aes(y = .data$ymax, colour = .data$Species)) + + geom_line(aes(y = .data$ymin, colour = .data$Species)) + + scale_colour_manual(values = linecolour) + + scale_fill_manual(values = linecolour) + + scale_y_continuous(trans = ytrans, breaks = ybreaks, + labels = prettyNum, name = ylab) + + scale_x_continuous(name = "Fishing effort") + + make_mizer_plot(p, c("Species", "Effort", "ymax", "ymin")) +} diff --git a/man/addPlot.Rd b/man/addPlot.Rd index de94906c6..050f2a3e9 100644 --- a/man/addPlot.Rd +++ b/man/addPlot.Rd @@ -92,6 +92,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/animate.Rd b/man/animate.Rd index 3643b56cf..575a304d9 100644 --- a/man/animate.Rd +++ b/man/animate.Rd @@ -149,6 +149,7 @@ Other plotting functions: \code{\link[=addPlot]{addPlot()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plot.Rd b/man/plot.Rd index f7aab4d4b..43adc491b 100644 --- a/man/plot.Rd +++ b/man/plot.Rd @@ -118,6 +118,7 @@ Other plotting functions: \code{\link[=addPlot]{addPlot()}}, \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plot2.Rd b/man/plot2.Rd index c4a3e0fdc..297bf1f67 100644 --- a/man/plot2.Rd +++ b/man/plot2.Rd @@ -106,6 +106,7 @@ Other plotting functions: \code{\link[=addPlot]{addPlot()}}, \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotBifurcation.Rd b/man/plotBifurcation.Rd new file mode 100644 index 000000000..7d092a0a8 --- /dev/null +++ b/man/plotBifurcation.Rd @@ -0,0 +1,128 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotBifurcation.R +\name{plotBifurcation} +\alias{plotBifurcation} +\title{Draw a bifurcation diagram over fishing effort} +\usage{ +plotBifurcation( + params, + effort = seq(0, 2, length.out = 21), + species = NULL, + value = c("biomass", "yield", "ssb"), + t_max = 100, + t_sample = NULL, + t_sample_default = 10, + t_save = 0.25, + continuation = TRUE, + return_data = FALSE, + progress_bar = TRUE, + ytrans = "log10" +) +} +\arguments{ +\item{params}{A \linkS4class{MizerParams} object.} + +\item{effort}{A numeric vector of fishing effort values for the x-axis. The +same effort is applied to every gear (as in \code{project(effort = value)}). +Default \code{seq(0, 2, length.out = 21)}.} + +\item{species}{The species to include. By default all target species. A +vector of species names or indices, as for other mizer plotting functions.} + +\item{value}{The quantity for the y-axis, one of \code{"biomass"} (default), +\code{"yield"} or \code{"ssb"}, computed with \code{\link[=getBiomass]{getBiomass()}}, \code{\link[=getYield]{getYield()}} or +\code{\link[=getSSB]{getSSB()}} respectively.} + +\item{t_max}{The maximum number of years to run the settling stage +(\code{\link[=projectToSteady]{projectToSteady()}}) at each effort value.} + +\item{t_sample}{The length in years of the window over which the settled +attractor is sampled to measure the envelope. If \code{NULL} (default) it is +chosen automatically: one full period for a detected limit cycle, or +\code{t_sample_default} years otherwise.} + +\item{t_sample_default}{The sampling window used when no limit-cycle period is +available (a stable or non-converged run). Default \code{10}.} + +\item{t_save}{The interval at which the sampling window is saved, controlling +how finely the cycle envelope is resolved. Default \code{0.25}.} + +\item{continuation}{If \code{TRUE} (default) each settling run warm-starts from the +attractor of the previous effort value.} + +\item{return_data}{If \code{TRUE} the data frame underlying the plot is returned +instead of the plot. Default \code{FALSE}.} + +\item{progress_bar}{If \code{TRUE} (default) a text progress bar is shown while the +effort values are swept.} + +\item{ytrans}{Transformation for the y-axis, \code{"log10"} (default) or +\code{"identity"}.} +} +\value{ +A ggplot2 object, or a data frame with columns \code{Effort}, \code{Species}, +\code{ymin}, \code{ymax} and \code{type} (the attractor type reported by +\code{\link[=projectToSteady]{projectToSteady()}}) if \code{return_data = TRUE}. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +Sweeps the fishing effort over a range of values and, for each value, runs +the full dynamics to their attractor and records the long-term range of a +summary quantity (biomass by default). The result is a bifurcation diagram +with fishing effort on the x-axis. + +For each effort value the attractor is found in two stages. First +\code{\link[=projectToSteady]{projectToSteady()}} runs the dynamics until they settle, stopping early once +it detects a stable steady state or a limit cycle and reporting which via its +\code{"convergence"} attribute. Then the settled state is projected forward once +more with \code{\link[=project]{project()}} over a short sampling window — one full period for a +limit cycle, or a few years otherwise — and the minimum and maximum of the +chosen quantity over that window are taken as the attractor envelope. + +For a stable steady state the minimum and maximum coincide and the species is +drawn as a single line; once the dynamics settle onto a limit cycle the two +separate and the species is drawn as a shaded band between them. The onset of +the band therefore marks the Hopf bifurcation (see +\code{vignette("dynamic_stability")}). + +By default the sweep uses \emph{continuation}: the projection at each effort value +starts from the attractor reached at the previous value rather than from the +original state. This shortens the transient and keeps the sweep on a single +branch of the attractor. Because it follows one branch, the diagram can look +different for an increasing versus a decreasing effort sequence if the model +has coexisting attractors (hysteresis); pass a decreasing \code{effort} vector to +trace the other direction. +} +\examples{ +\donttest{ +plotBifurcation(NS_params, effort = seq(0, 2, length.out = 11)) +} +} +\seealso{ +\code{\link[=getStability]{getStability()}}, \code{\link[=projectToSteady]{projectToSteady()}}, \code{\link[=plotBiomass]{plotBiomass()}} + +Other plotting functions: +\code{\link[=addPlot]{addPlot()}}, +\code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, +\code{\link{plot}}, +\code{\link[=plot2]{plot2()}}, +\code{\link[=plotBiomass]{plotBiomass()}}, +\code{\link[=plotCDF]{plotCDF()}}, +\code{\link[=plotCDF2]{plotCDF2()}}, +\code{\link[=plotDiet]{plotDiet()}}, +\code{\link[=plotFMort]{plotFMort()}}, +\code{\link[=plotFeedingLevel]{plotFeedingLevel()}}, +\code{\link[=plotGrowthCurves]{plotGrowthCurves()}}, +\code{\link{plotMizerParams}}, +\code{\link{plotMizerSim}}, +\code{\link[=plotPredMort]{plotPredMort()}}, +\code{\link[=plotRelative]{plotRelative()}}, +\code{\link[=plotSpectra]{plotSpectra()}}, +\code{\link[=plotSpectra2]{plotSpectra2()}}, +\code{\link[=plotSpectraRelative]{plotSpectraRelative()}}, +\code{\link[=plotYield]{plotYield()}}, +\code{\link[=plotYieldGear]{plotYieldGear()}}, +\code{\link{plotting_functions}} +} +\concept{plotting functions} diff --git a/man/plotBiomass.Rd b/man/plotBiomass.Rd index ffb95665f..c35fac639 100644 --- a/man/plotBiomass.Rd +++ b/man/plotBiomass.Rd @@ -108,6 +108,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, \code{\link[=plotDiet]{plotDiet()}}, diff --git a/man/plotCDF.Rd b/man/plotCDF.Rd index 25143566f..5423c5478 100644 --- a/man/plotCDF.Rd +++ b/man/plotCDF.Rd @@ -139,6 +139,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF2]{plotCDF2()}}, \code{\link[=plotDiet]{plotDiet()}}, diff --git a/man/plotCDF2.Rd b/man/plotCDF2.Rd index ef62d3e91..5a92e1229 100644 --- a/man/plotCDF2.Rd +++ b/man/plotCDF2.Rd @@ -127,6 +127,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotDiet]{plotDiet()}}, diff --git a/man/plotDiet.Rd b/man/plotDiet.Rd index 1b04d2e13..89c22f931 100644 --- a/man/plotDiet.Rd +++ b/man/plotDiet.Rd @@ -100,6 +100,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotFMort.Rd b/man/plotFMort.Rd index 145508672..a10e55590 100644 --- a/man/plotFMort.Rd +++ b/man/plotFMort.Rd @@ -97,6 +97,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotFeedingLevel.Rd b/man/plotFeedingLevel.Rd index ebbe272fe..60bd999e8 100644 --- a/man/plotFeedingLevel.Rd +++ b/man/plotFeedingLevel.Rd @@ -115,6 +115,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotGrowthCurves.Rd b/man/plotGrowthCurves.Rd index e6a988d06..9793c18b6 100644 --- a/man/plotGrowthCurves.Rd +++ b/man/plotGrowthCurves.Rd @@ -109,6 +109,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotM2.Rd b/man/plotM2.Rd index 24656c74b..a8c67fd8b 100644 --- a/man/plotM2.Rd +++ b/man/plotM2.Rd @@ -92,6 +92,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotMizerParams.Rd b/man/plotMizerParams.Rd index db7008fee..c1ef57028 100644 --- a/man/plotMizerParams.Rd +++ b/man/plotMizerParams.Rd @@ -36,6 +36,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotMizerSim.Rd b/man/plotMizerSim.Rd index 1031017fe..34ba94769 100644 --- a/man/plotMizerSim.Rd +++ b/man/plotMizerSim.Rd @@ -39,6 +39,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotPredMort.Rd b/man/plotPredMort.Rd index d79aa4f6a..b50a1cadd 100644 --- a/man/plotPredMort.Rd +++ b/man/plotPredMort.Rd @@ -96,6 +96,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotRelative.Rd b/man/plotRelative.Rd index cc238d61f..553f28a4b 100644 --- a/man/plotRelative.Rd +++ b/man/plotRelative.Rd @@ -94,6 +94,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotSpectra.Rd b/man/plotSpectra.Rd index dbf1e4515..f836975f7 100644 --- a/man/plotSpectra.Rd +++ b/man/plotSpectra.Rd @@ -145,6 +145,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotSpectra2.Rd b/man/plotSpectra2.Rd index 8523f22aa..109892e8c 100644 --- a/man/plotSpectra2.Rd +++ b/man/plotSpectra2.Rd @@ -122,6 +122,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotSpectraRelative.Rd b/man/plotSpectraRelative.Rd index c6c9d67db..09fefa253 100644 --- a/man/plotSpectraRelative.Rd +++ b/man/plotSpectraRelative.Rd @@ -103,6 +103,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotYield.Rd b/man/plotYield.Rd index bf3c6cb8c..27c54f581 100644 --- a/man/plotYield.Rd +++ b/man/plotYield.Rd @@ -92,6 +92,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotYieldGear.Rd b/man/plotYieldGear.Rd index 26eb16137..26f5203ae 100644 --- a/man/plotYieldGear.Rd +++ b/man/plotYieldGear.Rd @@ -94,6 +94,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/man/plotting_functions.Rd b/man/plotting_functions.Rd index 739a6d54c..b573c0149 100644 --- a/man/plotting_functions.Rd +++ b/man/plotting_functions.Rd @@ -126,6 +126,7 @@ Other plotting functions: \code{\link[=animate.ArrayTimeBySpeciesBySize]{animate.ArrayTimeBySpeciesBySize()}}, \code{\link{plot}}, \code{\link[=plot2]{plot2()}}, +\code{\link[=plotBifurcation]{plotBifurcation()}}, \code{\link[=plotBiomass]{plotBiomass()}}, \code{\link[=plotCDF]{plotCDF()}}, \code{\link[=plotCDF2]{plotCDF2()}}, diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index cd7ae0025..0632422f9 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -233,6 +233,7 @@ reference: - plotFeedingLevel - plotYield - plotYieldGear + - plotBifurcation - plotFMort - plotMizerParams - plotMizerSim diff --git a/tests/testthat/test-plotBifurcation.R b/tests/testthat/test-plotBifurcation.R new file mode 100644 index 000000000..daccd8a9a --- /dev/null +++ b/tests/testthat/test-plotBifurcation.R @@ -0,0 +1,41 @@ +# plotBifurcation ------------------------------------------------------------ + +test_that("plotBifurcation returns a tidy data frame with return_data", { + d <- plotBifurcation(NS_params_small, effort = c(0, 0.5, 1), + t_max = 6, t_sample_default = 3, t_save = 0.5, + progress_bar = FALSE, return_data = TRUE) + expect_s3_class(d, "data.frame") + expect_named(d, c("Effort", "Species", "ymin", "ymax", "type")) + # one row per (effort, species) + expect_equal(nrow(d), 3 * length(valid_species_arg(NS_params_small))) + expect_setequal(unique(d$Effort), c(0, 0.5, 1)) + # envelope is well ordered and positive for biomass + expect_true(all(d$ymax >= d$ymin)) + expect_true(all(d$ymin > 0)) +}) + +test_that("plotBifurcation returns a mizer_plot object", { + p <- plotBifurcation(NS_params_small, effort = c(0, 1), + t_max = 6, t_sample_default = 3, t_save = 0.5, + progress_bar = FALSE) + expect_s3_class(p, "mizer_plot") + expect_s3_class(p, "ggplot") +}) + +test_that("plotBifurcation respects the species and value arguments", { + d <- plotBifurcation(NS_params_small, effort = c(0, 1), species = "Cod", + value = "yield", t_max = 6, t_sample_default = 3, + t_save = 0.5, progress_bar = FALSE, return_data = TRUE) + expect_setequal(unique(d$Species), "Cod") +}) + +test_that("plotBifurcation validates its arguments", { + expect_error( + plotBifurcation(NS_params_small, effort = 1), + "length\\(effort\\) not greater than or equal to 2" + ) + expect_error( + plotBifurcation(NS_params_small, t_max = -1), + "t_max not greater than 0" + ) +}) From ea0b2d6fad628a1e29a113ea89b26d7a155354d4 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 19:49:26 +0100 Subject: [PATCH 10/15] Decouple limit-cycle amplitude floor and make extinction threshold relative projectToSteady() and steady() gain an amplitude_tol argument (default 0.01) for the relative biomass amplitude above which a persistent oscillation is reported as a limit cycle. Previously the fixed-point convergence tolerance `tol` was reused for this, conflating two quantities on different scales (a sum of squared log-differences versus a fraction of mean biomass), which meant small cycles near a Hopf bifurcation were not detected. The extinction check is now relative: a species is treated as going extinct once its reproduction rate falls below the extinction_threshold fraction (default 1e-6) of its value at the start of the run, replacing the hard-coded absolute 1e-20 floor. Under steady()'s frozen reproduction a healthy species is never flagged, and a zero-reproduction species is flagged immediately, so existing behaviour is preserved. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01T98LXny6oRNBMy9kHubpo4 --- NEWS.md | 7 +++- R/steady.R | 52 +++++++++++++++++++++------ man/projectToSteady.Rd | 18 ++++++++++ man/steady.Rd | 18 ++++++++++ tests/testthat/test-cycle_detection.R | 17 ++++++--- 5 files changed, 96 insertions(+), 16 deletions(-) diff --git a/NEWS.md b/NEWS.md index b0dd0f5e7..ca3be0461 100644 --- a/NEWS.md +++ b/NEWS.md @@ -84,7 +84,12 @@ These changes and how to adapt existing code are described in the new cycle period and relative amplitude when a cycle is found. Limit cycles are detected from a per-species biomass series sampled at the new `t_save` resolution (default `dt`), so detection no longer relies on the cycle period - being commensurate with `t_per`. + being commensurate with `t_per`. The relative-amplitude floor for calling an + oscillation a limit cycle is a separate `amplitude_tol` argument (default + `0.01`), independent of the fixed-point convergence tolerance `tol`, and a + species is treated as extinct once its reproduction falls below the + `extinction_threshold` fraction (default `1e-6`) of its value at the start of + the run. - New experimental `plotBifurcation()` draws a bifurcation diagram over fishing effort. For each effort value it follows the attractor of the full dynamics diff --git a/R/steady.R b/R/steady.R index f20b0c488..5df704046 100644 --- a/R/steady.R +++ b/R/steady.R @@ -116,6 +116,8 @@ projectToSteady <- function(params, dt = 0.1, t_save = dt, tol = 0.1 * t_per, + amplitude_tol = 0.01, + extinction_threshold = 1e-6, return_sim = FALSE, progress_bar = TRUE, info_level = 3, @@ -131,6 +133,8 @@ projectToSteady.MizerParams <- function(params, dt = 0.1, t_save = dt, tol = 0.1 * t_per, + amplitude_tol = 0.01, + extinction_threshold = 1e-6, return_sim = FALSE, progress_bar = TRUE, info_level = 3, @@ -140,7 +144,9 @@ projectToSteady.MizerParams <- function(params, effort <- validEffortVector(effort, params = params) params@initial_effort <- effort assert_that(t_max >= t_per, - tol > 0) + tol > 0, + amplitude_tol > 0, + extinction_threshold >= 0) if ((t_per < dt) || !isTRUE(all.equal((t_per - round(t_per / dt) * dt), 0))) { stop("t_per must be a positive multiple of dt") } @@ -185,6 +191,9 @@ projectToSteady.MizerParams <- function(params, n_other = params@initial_n_other, rates = r) + # Reference reproduction rate for the relative extinction test below. + rdd_start <- r$rdd + # Record a cheap scalar summary (per-species biomass) at the fine `t_save` # resolution so that a limit cycle can be detected and characterised even # when its period is incommensurate with `t_per`. @@ -231,8 +240,10 @@ projectToSteady.MizerParams <- function(params, sim@effort[i, ] <- params@initial_effort } - # Species with no reproduction are going extinct, so stop. - extinct <- is.na(current$rates$rdd) | current$rates$rdd <= 1e-20 + # A species whose reproduction has collapsed to a tiny fraction of its + # starting value is going extinct, so stop. + extinct <- is.na(current$rates$rdd) | + current$rates$rdd <= extinction_threshold * rdd_start if (any(extinct)) { warning(paste(params@species_params$species[extinct], collapse = ", "), " are going extinct.") @@ -250,7 +261,7 @@ projectToSteady.MizerParams <- function(params, } # Not settling to a fixed point: check whether we are on a limit cycle. cycle <- detect_limit_cycle(bio_series[seq_len(save_idx), , drop = FALSE], - t_save, tol) + t_save, amplitude_tol) if (!is.null(cycle)) { break } @@ -323,8 +334,8 @@ projectToSteady.MizerParams <- function(params, #' autocorrelation gives a candidate period (the first autocorrelation peak). The #' oscillation is accepted as a settled limit cycle only if it has already #' persisted for three periods: the relative amplitude of the three successive -#' period-windows must agree within `amp_rel_tol`, must exceed `tol`, and must -#' show no net decay across the three periods. The no-decay condition is what +#' period-windows must agree within `amp_rel_tol`, must exceed `amplitude_tol`, +#' and must show no net decay across the three periods. The no-decay condition is what #' distinguishes a genuine limit cycle from a slowly-decaying spiral toward a #' stable fixed point (which has a spectral radius just below 1). Discrimination #' is necessarily imperfect when the spectral radius is extremely close to 1, @@ -332,14 +343,15 @@ projectToSteady.MizerParams <- function(params, #' #' @param bio Numeric matrix of per-species biomass, one row per saved time step. #' @param t_save The sampling interval, so period `= lag * t_save`. -#' @param tol Minimum relative amplitude for an oscillation to count as a cycle. +#' @param amplitude_tol Minimum relative amplitude for an oscillation to count as +#' a cycle. #' @param acf_threshold Minimum autocorrelation at the candidate period. #' @param amp_rel_tol Maximum relative change of amplitude between successive #' periods for the cycle to count as settled. #' @return `NULL` if no settled cycle is detected, otherwise a list with the #' `period` (in the same time units as `t_save`) and the relative `amplitude`. #' @noRd -detect_limit_cycle <- function(bio, t_save, tol, +detect_limit_cycle <- function(bio, t_save, amplitude_tol, acf_threshold = 0.5, amp_rel_tol = 0.1) { n <- nrow(bio) @@ -355,7 +367,7 @@ detect_limit_cycle <- function(bio, t_save, tol, amp_old <- amp_window(bio[(n - 3 * w + 1):(n - 2 * w), , drop = FALSE]) amp_mid <- amp_window(bio[(n - 2 * w + 1):(n - w), , drop = FALSE]) amp_new <- amp_window(bio[(n - w + 1):n, , drop = FALSE]) - if (amp_new <= tol) return(NULL) + if (amp_new <= amplitude_tol) return(NULL) # Successive periods must have matching amplitude ... if (abs(amp_new - amp_mid) / amp_new > amp_rel_tol) return(NULL) if (abs(amp_mid - amp_old) / amp_mid > amp_rel_tol) return(NULL) @@ -420,6 +432,20 @@ find_first_acf_peak <- function(ac, threshold) { #' a small extra cost. Default is `dt`. #' @param tol The simulation stops when the relative change in the egg #' production RDI over `t_per` years is less than `tol` for every species. +#' @param amplitude_tol `r lifecycle::badge("experimental")` +#' The minimum relative biomass amplitude for a persistent oscillation to be +#' reported as a limit cycle rather than treated as an (effectively steady) +#' fixed point. This is a fraction of mean biomass and is kept separate from +#' `tol` (which measures convergence to a fixed point on a different scale). +#' Default `0.01`. +#' @param extinction_threshold `r lifecycle::badge("experimental")` +#' A species is treated as going extinct, stopping the run, once its +#' reproduction rate (RDD) falls below this fraction of its value at the start +#' of the run. For example the default `1e-6` treats a species as extinct once +#' its reproduction has collapsed to a millionth of its initial level. Because +#' it is relative to the initial reproduction, a species that starts with zero +#' reproduction is flagged immediately, and (in [steady()], where reproduction +#' is held constant) a healthy species is never flagged. #' @param return_sim If TRUE, the function returns the MizerSim object holding #' the result of the simulation run, saved at intervals of `t_per`. If FALSE (default) the function returns #' a MizerParams object with the "initial" slots set to the steady state. @@ -449,7 +475,8 @@ find_first_acf_peak <- function(ac, threshold) { #' plotSpectra(params) #' } steady <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, t_save = dt, - tol = 0.1 * dt, return_sim = FALSE, + tol = 0.1 * dt, amplitude_tol = 0.01, + extinction_threshold = 1e-6, return_sim = FALSE, preserve = c("reproduction_level", "erepro", "R_max"), progress_bar = TRUE, info_level = 3, @@ -460,7 +487,8 @@ steady <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, t_save = dt, #' @export steady.MizerParams <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, t_save = dt, - tol = 0.1 * dt, return_sim = FALSE, + tol = 0.1 * dt, amplitude_tol = 0.01, + extinction_threshold = 1e-6, return_sim = FALSE, preserve = c("reproduction_level", "erepro", "R_max"), progress_bar = TRUE, info_level = 3, @@ -496,6 +524,8 @@ steady.MizerParams <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, dt = dt, t_save = t_save, tol = tol, + amplitude_tol = amplitude_tol, + extinction_threshold = extinction_threshold, return_sim = return_sim, progress_bar = progress_bar, info_level = info_level, diff --git a/man/projectToSteady.Rd b/man/projectToSteady.Rd index 9b72df264..5a2214dbb 100644 --- a/man/projectToSteady.Rd +++ b/man/projectToSteady.Rd @@ -13,6 +13,8 @@ projectToSteady( dt = 0.1, t_save = dt, tol = 0.1 * t_per, + amplitude_tol = 0.01, + extinction_threshold = 1e-06, return_sim = FALSE, progress_bar = TRUE, info_level = 3, @@ -51,6 +53,22 @@ a small extra cost. Default is \code{dt}.} \item{tol}{The simulation stops when the relative change in the egg production RDI over \code{t_per} years is less than \code{tol} for every species.} +\item{amplitude_tol}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +The minimum relative biomass amplitude for a persistent oscillation to be +reported as a limit cycle rather than treated as an (effectively steady) +fixed point. This is a fraction of mean biomass and is kept separate from +\code{tol} (which measures convergence to a fixed point on a different scale). +Default \code{0.01}.} + +\item{extinction_threshold}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +A species is treated as going extinct, stopping the run, once its +reproduction rate (RDD) falls below this fraction of its value at the start +of the run. For example the default \code{1e-6} treats a species as extinct once +its reproduction has collapsed to a millionth of its initial level. Because +it is relative to the initial reproduction, a species that starts with zero +reproduction is flagged immediately, and (in \code{\link[=steady]{steady()}}, where reproduction +is held constant) a healthy species is never flagged.} + \item{return_sim}{If TRUE, the function returns the MizerSim object holding the result of the simulation run, saved at intervals of \code{t_per}. If FALSE (default) the function returns a MizerParams object with the "initial" slots set to the steady state.} diff --git a/man/steady.Rd b/man/steady.Rd index 71fea81a5..2f1c5baad 100644 --- a/man/steady.Rd +++ b/man/steady.Rd @@ -11,6 +11,8 @@ steady( dt = 0.1, t_save = dt, tol = 0.1 * dt, + amplitude_tol = 0.01, + extinction_threshold = 1e-06, return_sim = FALSE, preserve = c("reproduction_level", "erepro", "R_max"), progress_bar = TRUE, @@ -38,6 +40,22 @@ a small extra cost. Default is \code{dt}.} \item{tol}{The simulation stops when the relative change in the egg production RDI over \code{t_per} years is less than \code{tol} for every species.} +\item{amplitude_tol}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +The minimum relative biomass amplitude for a persistent oscillation to be +reported as a limit cycle rather than treated as an (effectively steady) +fixed point. This is a fraction of mean biomass and is kept separate from +\code{tol} (which measures convergence to a fixed point on a different scale). +Default \code{0.01}.} + +\item{extinction_threshold}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +A species is treated as going extinct, stopping the run, once its +reproduction rate (RDD) falls below this fraction of its value at the start +of the run. For example the default \code{1e-6} treats a species as extinct once +its reproduction has collapsed to a millionth of its initial level. Because +it is relative to the initial reproduction, a species that starts with zero +reproduction is flagged immediately, and (in \code{\link[=steady]{steady()}}, where reproduction +is held constant) a healthy species is never flagged.} + \item{return_sim}{If TRUE, the function returns the MizerSim object holding the result of the simulation run, saved at intervals of \code{t_per}. If FALSE (default) the function returns a MizerParams object with the "initial" slots set to the steady state.} diff --git a/tests/testthat/test-cycle_detection.R b/tests/testthat/test-cycle_detection.R index 4431bdf6b..bb9af0288 100644 --- a/tests/testthat/test-cycle_detection.R +++ b/tests/testthat/test-cycle_detection.R @@ -25,7 +25,7 @@ test_that("detect_limit_cycle finds a sustained oscillation", { t_save <- 0.25; period <- 5; a <- 0.1 idx <- 0:199 bio <- matrix(1000 * (1 + a * sin(2 * pi * idx * t_save / period)), ncol = 1) - cyc <- detect_limit_cycle(bio, t_save, tol = 0.01) + cyc <- detect_limit_cycle(bio, t_save, amplitude_tol = 0.01) expect_type(cyc, "list") expect_equal(cyc$period, period) # For a sinusoid about the mean the relative peak-to-trough amplitude is 2a. @@ -38,13 +38,22 @@ test_that("detect_limit_cycle rejects non-cycles", { # A decaying oscillation (spiral toward a fixed point) is not a cycle. bio_decay <- matrix(1000 * (1 + a * exp(-0.03 * idx) * sin(2 * pi * idx * t_save / period)), ncol = 1) - expect_null(detect_limit_cycle(bio_decay, t_save, tol = 0.01)) + expect_null(detect_limit_cycle(bio_decay, t_save, amplitude_tol = 0.01)) # A flat series is not a cycle. expect_null(detect_limit_cycle(matrix(1000, nrow = 200, ncol = 1), - t_save, tol = 0.01)) + t_save, amplitude_tol = 0.01)) # Too few samples. expect_null(detect_limit_cycle(matrix(1000, nrow = 10, ncol = 1), - t_save, tol = 0.01)) + t_save, amplitude_tol = 0.01)) +}) + +test_that("detect_limit_cycle respects amplitude_tol", { + t_save <- 0.25; period <- 5; a <- 0.03 # ~6% peak-to-trough amplitude + idx <- 0:199 + bio <- matrix(1000 * (1 + a * sin(2 * pi * idx * t_save / period)), ncol = 1) + # A 6% cycle counts when the floor is 1% but not when it is 10%. + expect_type(detect_limit_cycle(bio, t_save, amplitude_tol = 0.01), "list") + expect_null(detect_limit_cycle(bio, t_save, amplitude_tol = 0.1)) }) # Fixtures for the integration tests ----------------------------------------- From 2e008f3acfd6c8f622b1869859307b9940047e50 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 19:49:46 +0100 Subject: [PATCH 11/15] Expose settling tolerances in plotBifurcation() plotBifurcation() forwards tol, amplitude_tol and extinction_threshold to the projectToSteady() settling stage. The default tol is set tighter (0.01) than projectToSteady()'s own default so the stable branch settles more fully, giving clean single lines below the bifurcation instead of residual-transient bands, while the lower amplitude_tol resolves the cycle onset earlier. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01T98LXny6oRNBMy9kHubpo4 --- NEWS.md | 3 ++- R/plotBifurcation.R | 19 ++++++++++++++++++- man/plotBifurcation.Rd | 15 +++++++++++++++ 3 files changed, 35 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index ca3be0461..5abaf4fc4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -96,7 +96,8 @@ These changes and how to adapt existing code are described in the new and plots the long-term range of a summary quantity (biomass, yield or SSB). A stable steady state appears as a single line and a limit cycle as a band between the minimum and maximum, so a Hopf bifurcation shows up as the effort - at which the band opens up. + at which the band opens up. The settling stage runs `projectToSteady()`, whose + `tol`, `amplitude_tol` and `extinction_threshold` are exposed for tuning. - `project()` gains a `callback` argument for a user-defined function to be called at each saved time step. diff --git a/R/plotBifurcation.R b/R/plotBifurcation.R index ebeeb5d78..103bc7a68 100644 --- a/R/plotBifurcation.R +++ b/R/plotBifurcation.R @@ -58,6 +58,15 @@ #' available (a stable or non-converged run). Default `10`. #' @param t_save The interval at which the sampling window is saved, controlling #' how finely the cycle envelope is resolved. Default `0.25`. +#' @param tol Convergence tolerance for the settling stage, passed to +#' [projectToSteady()]. Tighter (smaller) values settle the stable branch more +#' fully, giving cleaner single lines below the bifurcation at the cost of +#' longer runs. Default `0.01`. +#' @param amplitude_tol The minimum relative biomass amplitude for a run to be +#' classified as a limit cycle, passed to [projectToSteady()]. Default `0.01`. +#' @param extinction_threshold The relative reproduction collapse below which a +#' species is treated as extinct, passed to [projectToSteady()]. Default +#' `1e-6`. #' @param continuation If `TRUE` (default) each settling run warm-starts from the #' attractor of the previous effort value. #' @param return_data If `TRUE` the data frame underlying the plot is returned @@ -84,6 +93,9 @@ plotBifurcation <- function(params, t_sample = NULL, t_sample_default = 10, t_save = 0.25, + tol = 0.01, + amplitude_tol = 0.01, + extinction_threshold = 1e-6, continuation = TRUE, return_data = FALSE, progress_bar = TRUE, @@ -93,7 +105,10 @@ plotBifurcation <- function(params, assert_that(is.numeric(effort), length(effort) >= 2, is.number(t_max), t_max > 0, is.null(t_sample) || is.number(t_sample), - is.number(t_sample_default), t_sample_default > 0) + is.number(t_sample_default), t_sample_default > 0, + is.number(tol), tol > 0, + is.number(amplitude_tol), amplitude_tol > 0, + is.number(extinction_threshold), extinction_threshold >= 0) species <- valid_species_arg(params, species, error_on_empty = TRUE) value_func <- switch(value, @@ -116,6 +131,8 @@ plotBifurcation <- function(params, # Stage 1: settle onto the attractor, stopping early once a steady state # or limit cycle is detected. settled <- projectToSteady(p_run, effort = effort[i], t_max = t_max, + tol = tol, amplitude_tol = amplitude_tol, + extinction_threshold = extinction_threshold, return_sim = FALSE, progress_bar = FALSE, info_level = 0) conv <- attr(settled, "convergence") diff --git a/man/plotBifurcation.Rd b/man/plotBifurcation.Rd index 7d092a0a8..3c346f439 100644 --- a/man/plotBifurcation.Rd +++ b/man/plotBifurcation.Rd @@ -13,6 +13,9 @@ plotBifurcation( t_sample = NULL, t_sample_default = 10, t_save = 0.25, + tol = 0.01, + amplitude_tol = 0.01, + extinction_threshold = 1e-06, continuation = TRUE, return_data = FALSE, progress_bar = TRUE, @@ -47,6 +50,18 @@ available (a stable or non-converged run). Default \code{10}.} \item{t_save}{The interval at which the sampling window is saved, controlling how finely the cycle envelope is resolved. Default \code{0.25}.} +\item{tol}{Convergence tolerance for the settling stage, passed to +\code{\link[=projectToSteady]{projectToSteady()}}. Tighter (smaller) values settle the stable branch more +fully, giving cleaner single lines below the bifurcation at the cost of +longer runs. Default \code{0.01}.} + +\item{amplitude_tol}{The minimum relative biomass amplitude for a run to be +classified as a limit cycle, passed to \code{\link[=projectToSteady]{projectToSteady()}}. Default \code{0.01}.} + +\item{extinction_threshold}{The relative reproduction collapse below which a +species is treated as extinct, passed to \code{\link[=projectToSteady]{projectToSteady()}}. Default +\code{1e-6}.} + \item{continuation}{If \code{TRUE} (default) each settling run warm-starts from the attractor of the previous effort value.} From ed19d267643aa6910939fde30906749dab2a05cc Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Tue, 14 Jul 2026 21:19:12 +0100 Subject: [PATCH 12/15] Fix false positive extinction warnings in steadyNewton When calculating steady states with dynamic reproduction, steadyNewton() previously flagged a species as extinct if ANY of its size classes dropped below the numerical floor. This caused false positives for healthy species whose large-size tails naturally shortened. Furthermore, species that did trigger the warning were incorrectly left at their numerical abundance floor rather than being fully removed. This commit updates the extinction check to only inspect the recruit (egg) abundance, mirroring the behavior of the time-stepping steady() function. If the recruit abundance collapses, the species is correctly flagged as extinct and its abundance vector is entirely zeroed out. Tail-class shortening for healthy species is still safely truncated without triggering the warning. Tests are updated to assert the new zeroed-out behaviour. --- R/steadyNewton.R | 37 ++++++++++++++++++------------ tests/testthat/test-steadyNewton.R | 9 ++++---- 2 files changed, 26 insertions(+), 20 deletions(-) diff --git a/R/steadyNewton.R b/R/steadyNewton.R index 7c880b620..ed8f989e3 100644 --- a/R/steadyNewton.R +++ b/R/steadyNewton.R @@ -181,36 +181,43 @@ steadyNewton.MizerParams <- function(params, } N <- active$unpack(sol$x) - n_pp <- resource_steady_semichemostat(params, N, n_other) - - # Zero out any densities that ended up at or near the penalty floor - for (i in seq_len(nrow(N))) { - floor_val <- max(N[i, ]) * 1.1 * 1e-15 - N[i, N[i, ] < floor_val] <- 0 - } - params@initial_n[] <- N - params@initial_n_pp[] <- n_pp + is_extinct <- rep(FALSE, nrow(N)) + names(is_extinct) <- rownames(N) # Check for extinctions if using the relative floor if (reproduction == "dynamic" && !is.null(extinction_floor) && extinction_floor > 0) { extinct_threshold <- extinction_floor * 1.01 - is_extinct <- rep(FALSE, nrow(N)) - names(is_extinct) <- rownames(N) for (i in seq_len(nrow(N))) { lo <- params@w_min_idx[i] - pos <- lo:active$w_top[i] - pos <- pos[N0_initial[i, pos] > 0] - if (length(pos) && any(N[i, pos] / N0_initial[i, pos] <= extinct_threshold)) { + # A species is considered extinct if its recruitment (the egg class) + # has dropped to the relative abundance floor. Checking any() would + # produce false positives when the tail naturally shortens. + if (N0_initial[i, lo] > 0 && + N[i, lo] / N0_initial[i, lo] <= extinct_threshold) { is_extinct[i] <- TRUE + N[i, ] <- 0 } } if (any(is_extinct)) { - warning("The following species went extinct and were pegged to their abundance floor: ", + warning("The following species went extinct and were set to zero: ", paste(names(is_extinct)[is_extinct], collapse = ", ")) } } + # Zero out any tail densities that ended up at or near the structural penalty floor + for (i in seq_len(nrow(N))) { + if (!is_extinct[i]) { + floor_val <- max(N0_initial[i, ]) * 1.1 * 1e-15 + N[i, N[i, ] < floor_val] <- 0 + } + } + + n_pp <- resource_steady_semichemostat(params, N, n_other) + + params@initial_n[] <- N + params@initial_n_pp[] <- n_pp + # Restore density-dependent reproduction, just as steady() does. if (params@rates_funcs$RDD == "BevertonHoltRDD" && reproduction == "fixed") { if (preserve == "reproduction_level") { diff --git a/tests/testthat/test-steadyNewton.R b/tests/testthat/test-steadyNewton.R index cc661ed3c..50edb3f7a 100644 --- a/tests/testthat/test-steadyNewton.R +++ b/tests/testthat/test-steadyNewton.R @@ -78,14 +78,13 @@ test_that("steadyNewton handles extinctions under reproduction = 'dynamic' with p_extinct@species_params$erepro[3] <- 1e-12 p_extinct@initial_n[3, ] <- p_steady@initial_n[3, ] - # Verify that the solver issues the extinction warning and pegs to the floor + # Verify that the solver issues the extinction warning and zeroes out the species expect_warning(pn_ext <- steadyNewton(p_extinct, reproduction = "dynamic", extinction_floor = 1e-6), - "went extinct and were pegged to their abundance floor") + "went extinct and were set to zero") - # Verify Cod abundance is pegged exactly to the floor (1e-6 of its initial abundance) + # Verify Cod abundance is exactly 0 lo <- p_extinct@w_min_idx[3] - ratio <- pn_ext@initial_n[3, lo] / p_extinct@initial_n[3, lo] - expect_equal(ratio, 1e-6, tolerance = 1e-2) + expect_equal(pn_ext@initial_n[3, lo], 0) }) test_that("steadyNewton errors for unsupported resource dynamics", { From 1932817c2437dea5a9b4c7831132dbdf5af3dc1c Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Wed, 15 Jul 2026 07:59:47 +0100 Subject: [PATCH 13/15] Replace n_points with t_save in getLimitCycleSim() Align getLimitCycleSim() with other functions that return a MizerSim by taking a t_save time-step interval (default 0.1) instead of a n_points count. The time grid now spans one full period at the requested spacing. Co-Authored-By: Claude Opus 4.8 --- R/getLimitCycleSim.R | 15 +++++---------- man/getLimitCycleSim.Rd | 8 +++----- tests/testthat/test-getLimitCycleSim.R | 13 +++++++++---- 3 files changed, 17 insertions(+), 19 deletions(-) diff --git a/R/getLimitCycleSim.R b/R/getLimitCycleSim.R index 71f4112e8..f7c98bf01 100644 --- a/R/getLimitCycleSim.R +++ b/R/getLimitCycleSim.R @@ -43,17 +43,15 @@ #' @param amplitude Maximum relative perturbation #' \eqn{\max_w |\delta N(t,w)|/N^*(w)} across the limit cycle. Default #' `0.1`. -#' @param n_points Number of time points in the returned -#' \linkS4class{MizerSim}. Defaults to `ceiling(T) + 1`, giving one point -#' per time step (year) over exactly one period. Increase for smoother -#' plots, e.g. `n_points = 200`. +#' @param t_save The time interval between saved time steps in the returned +#' \linkS4class{MizerSim}. Defaults to `0.1`. #' @param ... Additional arguments forwarded to [getStability()] when #' stability has not already been computed and cached on `params`. #' @return A \linkS4class{MizerSim} object whose time axis spans one period #' \eqn{[0, T]} of the linearised limit cycle. #' @seealso [getStability()], [steadyNewton()] #' @export -getLimitCycleSim <- function(params, amplitude = 0.1, n_points = NULL, ...) { +getLimitCycleSim <- function(params, amplitude = 0.1, t_save = 0.1, ...) { # ------------------------------------------------------------------ # 1. Get (or compute) stability analysis # ------------------------------------------------------------------ @@ -98,11 +96,8 @@ getLimitCycleSim <- function(params, amplitude = 0.1, n_points = NULL, ...) { # ------------------------------------------------------------------ # 3. Time grid # ------------------------------------------------------------------ - if (is.null(n_points)) { - t_seq <- seq(0, ceiling(T_period)) # integer steps, one per year - } else { - t_seq <- seq(0, T_period, length.out = as.integer(n_points)) - } + n_steps <- ceiling(T_period / t_save) + t_seq <- seq(0, by = t_save, length.out = n_steps + 1) # ------------------------------------------------------------------ # 4. Amplitude scaling diff --git a/man/getLimitCycleSim.Rd b/man/getLimitCycleSim.Rd index 85d59a76d..274917742 100644 --- a/man/getLimitCycleSim.Rd +++ b/man/getLimitCycleSim.Rd @@ -4,7 +4,7 @@ \alias{getLimitCycleSim} \title{Construct a MizerSim of the linearised limit cycle} \usage{ -getLimitCycleSim(params, amplitude = 0.1, n_points = NULL, ...) +getLimitCycleSim(params, amplitude = 0.1, t_save = 0.1, ...) } \arguments{ \item{params}{A \linkS4class{MizerParams} object at a steady state, @@ -16,10 +16,8 @@ exists and contains \code{leading_eigenvectors} it is used directly; otherwise \eqn{\max_w |\delta N(t,w)|/N^*(w)} across the limit cycle. Default \code{0.1}.} -\item{n_points}{Number of time points in the returned -\linkS4class{MizerSim}. Defaults to \code{ceiling(T) + 1}, giving one point -per time step (year) over exactly one period. Increase for smoother -plots, e.g. \code{n_points = 200}.} +\item{t_save}{The time interval between saved time steps in the returned +\linkS4class{MizerSim}. Defaults to \code{0.1}.} \item{...}{Additional arguments forwarded to \code{\link[=getStability]{getStability()}} when stability has not already been computed and cached on \code{params}.} diff --git a/tests/testthat/test-getLimitCycleSim.R b/tests/testthat/test-getLimitCycleSim.R index d1f9027c7..a05d1d94c 100644 --- a/tests/testthat/test-getLimitCycleSim.R +++ b/tests/testthat/test-getLimitCycleSim.R @@ -45,15 +45,20 @@ test_that("getLimitCycleSim abundances are non-negative", { expect_true(all(lcs@n_pp >= 0)) }) -test_that("getLimitCycleSim n_points controls number of time steps", { +test_that("getLimitCycleSim t_save controls the time step spacing", { pn <- steadyNewton(p_steady_lcs, stability = TRUE) stab <- attr(pn, "stability") skip_if(is.null(stab$hopf_period)) skip_if(abs(Im(stab$eigenvalues[1])) <= 1e-8) - lcs <- getLimitCycleSim(pn, n_points = 50) - expect_equal(length(getTimes(lcs)), 50L) + T_period <- stab$dominant_period + t_save <- T_period / 50 + lcs <- getLimitCycleSim(pn, t_save = t_save) + times <- getTimes(lcs) + + expect_equal(length(times), ceiling(T_period / t_save) + 1L) + expect_equal(diff(times), rep(t_save, length(times) - 1L)) }) test_that("getLimitCycleSim n array has correct species and size dimnames", { @@ -76,7 +81,7 @@ test_that("getLimitCycleSim respects amplitude: max relative perturbation ~ ampl skip_if(abs(Im(stab$eigenvalues[1])) <= 1e-8) amp <- 0.1 - lcs <- getLimitCycleSim(pn, amplitude = amp, n_points = 200) + lcs <- getLimitCycleSim(pn, amplitude = amp, t_save = stab$dominant_period / 200) N_ss <- pn@initial_n active <- N_ss > 0 max_rel <- max(abs(lcs@n - rep(N_ss, each = dim(lcs@n)[1])) / From feda7ef50113d25c9adba1d5872179e9d5ed1d49 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Wed, 15 Jul 2026 08:56:10 +0100 Subject: [PATCH 14/15] Added explanations of the numerical scheme used to determine stability. --- R/steadyNewton.R | 11 ++++++-- man/getStability.Rd | 13 ++++++--- vignettes/dynamic_stability.Rmd | 18 ++++++++++--- vignettes/numerical_details.qmd | 48 +++++++++++++++++++++++++++++++++ 4 files changed, 81 insertions(+), 9 deletions(-) diff --git a/R/steadyNewton.R b/R/steadyNewton.R index ed8f989e3..1845f893b 100644 --- a/R/steadyNewton.R +++ b/R/steadyNewton.R @@ -484,8 +484,12 @@ steady_state_residual <- function(params, rdd_const, n_other, effort, active, #' \deqn{A(N^t, n_{pp}^t)\,N^{t+1} = S(N^t, n_{pp}^t),} #' and an exact semi-chemostat update for the resource: #' \deqn{n_{pp}^{t+1} = n_{pp}^* + (n_{pp}^t - n_{pp}^*)\,e^{-\mu^t\,dt},} -#' where \eqn{n_{pp}^* = r_{pp}\,c_{pp}/\mu^t} is the resource steady state +#' where \eqn{n_{pp}^*} = r_{pp}\,c_{pp}/\mu^t is the resource steady state #' conditioned on the mortality \eqn{\mu^t} due to consumers at time \eqn{t}. +#' Note that this function evaluates the Jacobian of this specific first-order +#' backward-Euler time step, regardless of which `method` you might later pass +#' to [project()]. However, it fully respects any higher-order spatial scheme +#' configured via [second_order_w()]. #' #' The stability is determined by the Jacobian of the full one-step-ahead map #' \eqn{G : (N, n_{pp}) \mapsto (N^{t+1}, n_{pp}^{t+1})} at the fixed point. @@ -513,7 +517,10 @@ steady_state_residual <- function(params, rdd_const, n_other, effort, active, #' For the default mizer time step of one year this is in years. #' #' Both branches use the same `project_n_loop()` C++ Thomas solver as the -#' regular dynamics and centred finite differences with relative step \eqn{h}. +#' regular dynamics, evaluating the transport coefficients with the exact +#' spatial scheme configured in `params` (e.g., first-order upwind or a +#' second-order limiter). The Jacobian is computed numerically using a +#' multiplicative (relative) finite-difference step \eqn{h \cdot N^*}. #' #' @param params A \linkS4class{MizerParams} object whose `initial_n` holds the #' steady state to analyse. Typically the output of [steadyNewton()]. diff --git a/man/getStability.Rd b/man/getStability.Rd index 71222bc01..d9a3509c5 100644 --- a/man/getStability.Rd +++ b/man/getStability.Rd @@ -8,7 +8,7 @@ getStability( params, reproduction = c("fixed", "dynamic"), effort = params@initial_effort, - include_resource = FALSE, + include_resource = TRUE, extinction_floor = 1e-06, h = 1e-04 ) @@ -81,8 +81,12 @@ The mizer time step applies a backward-Euler transport solve for the fish: \deqn{A(N^t, n_{pp}^t)\,N^{t+1} = S(N^t, n_{pp}^t),} and an exact semi-chemostat update for the resource: \deqn{n_{pp}^{t+1} = n_{pp}^* + (n_{pp}^t - n_{pp}^*)\,e^{-\mu^t\,dt},} -where \eqn{n_{pp}^* = r_{pp}\,c_{pp}/\mu^t} is the resource steady state +where \eqn{n_{pp}^*} = r_{pp}\,c_{pp}/\mu^t is the resource steady state conditioned on the mortality \eqn{\mu^t} due to consumers at time \eqn{t}. +Note that this function evaluates the Jacobian of this specific first-order +backward-Euler time step, regardless of which \code{method} you might later pass +to \code{\link[=project]{project()}}. However, it fully respects any higher-order spatial scheme +configured via \code{\link[=second_order_w]{second_order_w()}}. The stability is determined by the Jacobian of the full one-step-ahead map \eqn{G : (N, n_{pp}) \mapsto (N^{t+1}, n_{pp}^{t+1})} at the fixed point. @@ -110,7 +114,10 @@ crosses the unit circle, giving a limit-cycle period For the default mizer time step of one year this is in years. Both branches use the same \code{project_n_loop()} C++ Thomas solver as the -regular dynamics and centred finite differences with relative step \eqn{h}. +regular dynamics, evaluating the transport coefficients with the exact +spatial scheme configured in \code{params} (e.g., first-order upwind or a +second-order limiter). The Jacobian is computed numerically using a +multiplicative (relative) finite-difference step \eqn{h \cdot N^*}. } } \seealso{ diff --git a/vignettes/dynamic_stability.Rmd b/vignettes/dynamic_stability.Rmd index e028a6c37..c8be7e0c5 100644 --- a/vignettes/dynamic_stability.Rmd +++ b/vignettes/dynamic_stability.Rmd @@ -35,13 +35,23 @@ $$ N(t+\Delta t) = G(N(t)) $$ where $N(t)$ represents the full state vector (abundances of all species at all sizes) at time $t$. At a steady state $N^*$, we have $N^* = G(N^*)$. -To determine if $N^*$ is stable, we linearise the dynamics for small perturbations $\delta N(t) = N(t) - N^*$ around the steady state: +Because abundances in a mizer model span dozens of orders of magnitude across the size spectrum, evaluating stability in terms of absolute perturbations (e.g., $N(t) = N^* + \delta N(t)$) can be numerically poorly-scaled. Instead, it is more natural to consider **multiplicative (relative) perturbations**, where we perturb the steady state by a small relative amount $x(t)$: -$$ \delta N(t+\Delta t) \approx L \, \delta N(t) $$ +$$ N_i(t) = N_i^* \left(1 + x_i(t)\right) $$ -where $L$ is the Jacobian matrix evaluated at $N^*$, with elements $L_{ij} = \frac{\partial G_i}{\partial N_j}(N^*)$. +Linearising the dynamics for these small relative perturbations $x(t)$ gives: -The stability of the steady state is determined by the eigenvalues $\lambda_i$ of $L$. The **spectral radius** is the maximum of the absolute values (moduli) of these eigenvalues: $\rho = \max_i |\lambda_i|$. If $\rho < 1$ every perturbation decays and the steady state is stable; if $\rho > 1$ some perturbation grows and it is unstable. +$$ x(t+\Delta t) \approx K \, x(t) $$ + +where $K$ is the relative Jacobian matrix. The elements of $K$ describe the proportional change in species $i$ resulting from a proportional change in species $j$: + +$$ K_{ij} = \frac{\partial \log G_i}{\partial \log N_j}(N^*) = \frac{N_j^*}{N_i^*} \frac{\partial G_i}{\partial N_j}(N^*) $$ + +Notice that $K$ is related to the standard absolute Jacobian $L$ (with elements $L_{ij} = \frac{\partial G_i}{\partial N_j}$) by a similarity transform $K = D^{-1} L D$, where $D$ is a diagonal matrix of the steady-state abundances $N^*$. Because $K$ and $L$ are similar matrices, **they have exactly the same eigenvalues**. + +The stability of the steady state is determined by the eigenvalues $\lambda_i$ of these matrices. The **spectral radius** is the maximum of the absolute values (moduli) of these eigenvalues: $\rho = \max_i |\lambda_i|$. If $\rho < 1$ every perturbation decays and the steady state is stable; if $\rho > 1$ some perturbation grows and it is unstable. + +*(Note on numerical implementation: While analyzing the relative Jacobian $K$ is conceptually natural, computing it requires dividing by $N_i^*$, which causes floating-point overflow for structurally zero size classes (e.g., extinct species or truncated tails). Therefore, `getStability()` numerically computes the absolute Jacobian $L$ using finite-difference steps scaled multiplicatively to the local density. This is mathematically equivalent to the relative analysis but avoids dividing by zero, ensuring robust eigenvalue computation.)* ## Driving the North Sea model unstable diff --git a/vignettes/numerical_details.qmd b/vignettes/numerical_details.qmd index 4c2f3563b..d48bee59e 100644 --- a/vignettes/numerical_details.qmd +++ b/vignettes/numerical_details.qmd @@ -684,3 +684,51 @@ R_{dd}^t = \hat{R}_{dd}^{t+1} = R_{dd}^{t+1/2}. $$ {#eq-ss-rates-equal} The rates used in the steady-state calculation are therefore just the rates evaluated at $N^*$. The predictor-corrector method affects the transient path to the steady state, not the steady state itself. + +## Direct Steady-State Solver (`steadyNewton`) + +The function `steadyNewton()` finds the steady state by directly solving the discretised algebraic equation $F(N) = 0$ (derived in the sections above) using a Newton-type root finder. This allows it to converge to the steady state even when the state is dynamically unstable. + +To make the nonlinear algebraic system well-behaved for the root finder, `steadyNewton()` employs several numerical strategies: + +1. **Log-space solve:** The consumer densities are solved for in logarithmic space ($x = \log N$). This enforces strict positivity of the densities during the iteration and natively conditions the Jacobian by turning absolute density perturbations—which span dozens of orders of magnitude—into well-scaled $O(1)$ relative (multiplicative) perturbations. The residual of the equation is correspondingly scaled by $1/N$, turning it into a per-capita rate of change that is also $O(1)$ across all size classes. + +2. **Active support and structural zeros:** The solver defines the "active" size classes for each species as running from the egg size up to the overall grid truncation limit, regardless of whether they initially carry non-zero density. Because the logarithm of zero is $-\infty$, the solver applies a smooth, continuous penalty floor to the log-abundances. This automatically bounds the unknowns, gracefully handling structurally zero-density classes and preventing singular Jacobians. After convergence, densities that remain trapped at or near this penalty floor are explicitly set to exactly zero. This allows the solver to automatically discover the true support (the highest non-zero size class) of the steady state. + +3. **Interaction with high-order spatial schemes:** The solver strictly respects the spatial transport scheme configured in the parameters (e.g., via `second_order_w()`). The nonlinear flux limiter weight is recalculated at every residual evaluation, ensuring the returned steady state is an exact fixed point of the chosen dynamics. However, because the TVD "van Leer" flux limiter is only Lipschitz-continuous (its derivative is discontinuous where the smoothness indicator equals 1), the Newton iteration converges to the fixed point but typically not to machine precision. Conversely, the unlimited "centred" reconstruction is perfectly smooth, but at a steady state with no physical diffusion it admits an undamped odd-even mode, leading to a structurally ill-conditioned Jacobian for which the Newton solver is unlikely to converge. + +## Linear Stability Analysis + +Once a steady state $N^*$ has been found, mizer can evaluate its dynamic stability by linearising the discrete-time map around $N^*$. The time step acts as a map $N^{t+1} = G(N^t)$, and the stability of the steady state is governed by the eigenvalues of the absolute Jacobian matrix $L$: +$$ +L_{ij} = \frac{\partial G_i}{\partial N_j}(N^*) +$$ {#eq-absolute-jacobian} + +### Multiplicative Finite Differences +Because abundances in a mizer model span dozens of orders of magnitude, evaluating the Jacobian $L$ numerically using a constant absolute perturbation size (e.g., $N_j^* + \epsilon$ with $\epsilon = 10^{-4}$) is poorly scaled. The perturbation would be lost in rounding error for larvae, and catastrophically huge for large fish. + +Instead, the `getStability()` function approximates the derivative using a **multiplicative (relative)** finite-difference step. It perturbs the abundance $N_j^*$ by an amount $\Delta N_j$ proportional to the steady-state abundance: +$$ +\Delta N_j = h \cdot \max(N_j^*, \text{machine eps}) +$$ {#eq-multiplicative-step} +where $h$ is the relative step size (default $10^{-4}$). The function then evaluates the system's response: +$$ +L_{ij} \approx \frac{G_i(N^* + \Delta N_j) - G_i(N^* - \Delta N_j)}{2 \Delta N_j} +$$ {#eq-jacobian-fd} +This ensures the numerical precision is maintained across all size classes. A tiny additive floor (the machine epsilon) is kept to allow the solver to gently "poke" classes where the steady-state abundance is exactly zero (like extinct species or the truncated tail above $w_{max}$), capturing their stability against small invasions. + +### Equivalence to Log-Space Analysis +It might seem conceptually more natural to analyze stability entirely in log-space, by perturbing the relative abundances $x_j(t) = (N_j(t) - N_j^*) / N_j^*$ and computing the relative Jacobian $K$: +$$ +K_{ij} = \frac{\partial \log G_i}{\partial \log N_j} +$$ {#eq-relative-jacobian} + +Applying the chain rule, this log-space Jacobian relates to the absolute Jacobian via: +$$ +K_{ij} = \frac{\partial \log G_i}{\partial G_i} \frac{\partial G_i}{\partial N_j} \frac{\partial N_j}{\partial \log N_j} = \frac{1}{N_i^*} L_{ij} N_j^* +$$ {#eq-jacobian-similarity} +In matrix notation, this is a **similarity transform**: $K = D^{-1} L D$, where $D$ is a diagonal matrix containing the steady-state abundances $N^*$. + +A fundamental theorem of linear algebra states that similar matrices have **exactly the same eigenvalues**. Because stability is determined entirely by whether any eigenvalue's modulus exceeds 1, $L$ and $K$ yield mathematically identical stability conclusions. + +Numerically, however, explicitly computing the relative Jacobian $K$ requires dividing the elements of $L$ by the steady-state abundances $N_i^*$. For structurally zero size classes, this division by zero causes floating-point overflow (`Inf` or `NaN`), breaking standard eigenvalue solvers. Therefore, `getStability()` numerically computes the well-conditioned absolute Jacobian $L$ (using the multiplicative-scale finite difference) and relies on the similarity equivalence to guarantee the correct result. From 02e5eaba511db8efefe328c219e4013c0116e66d Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Wed, 15 Jul 2026 10:34:03 +0100 Subject: [PATCH 15/15] Restore steadyNewton docs and export after rebase onto master The rebase onto master (which removed steadyNewton) left the branch's generated NAMESPACE without the steadyNewton export and stripped the steadyNewton cross-references from helpers.R, transport.R and the numerical_details / cheatsheet vignettes. Restore them so the branch remains a complete, self-consistent feature branch for the dynamic-stability work, and regenerate the man pages (which also corrects a stale include_resource default in getStability.Rd). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01T98LXny6oRNBMy9kHubpo4 --- NAMESPACE | 2 ++ R/helpers.R | 6 +++--- R/transport.R | 3 ++- man/getStability.Rd | 2 +- man/get_steady_state_n.Rd | 6 +++--- vignettes/cheatsheet-model-setup-and-calibration.Rmd | 1 + vignettes/numerical_details.qmd | 6 +++++- 7 files changed, 17 insertions(+), 9 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 6c91ad8ce..99eee7ac8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -256,6 +256,7 @@ S3method(species_params,MizerSim) S3method(species_params,data.frame) S3method(species_params,species_params) S3method(steady,MizerParams) +S3method(steadyNewton,MizerParams) S3method(steadySingleSpecies,MizerParams) S3method(str,ArrayResourceBySize) S3method(str,ArraySpeciesBySize) @@ -595,6 +596,7 @@ export(sigmoid_length) export(sigmoid_weight) export(species_params) export(steady) +export(steadyNewton) export(steadySingleSpecies) export(truncated_lognormal_pred_kernel) export(use_predation_diffusion) diff --git a/R/helpers.R b/R/helpers.R index 983abd257..b0c4000ab 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -234,9 +234,9 @@ w2l <- function(w, species_params) { #' #' The returned abundance is held at zero above each species' `w_max`, the same #' upper boundary condition that [project()] imposes (via `zero_above_support()` -#' in `project_n()`). Without this the bottom-up solve would carry density above -#' `w_max` whenever growth is still positive there or diffusion pushes density -#' past it. +#' in `project_n()`) and that [steadyNewton()] solves on. Without this the +#' bottom-up solve would carry density above `w_max` whenever growth is still +#' positive there or diffusion pushes density past it. #' @concept helper get_steady_state_n <- function(params, g, mu, D, N0, max_iterations = 500, tol = 1e-10, diff --git a/R/transport.R b/R/transport.R index 5d26dd69b..6d9a1236a 100644 --- a/R/transport.R +++ b/R/transport.R @@ -55,7 +55,8 @@ log_dx <- function(params) { #' #' The dynamics impose this as the upper boundary condition (see #' `get_transport_coefs()`) so that abundance is held at zero above `w_max` even -#' when diffusion would otherwise carry density past it. +#' when diffusion would otherwise carry density past it, and [steadyNewton()] +#' solves on exactly this support. #' #' @param params A \linkS4class{MizerParams} object. #' @return A per-species integer vector of the top retained size-class index. diff --git a/man/getStability.Rd b/man/getStability.Rd index d9a3509c5..1f0f282bc 100644 --- a/man/getStability.Rd +++ b/man/getStability.Rd @@ -8,7 +8,7 @@ getStability( params, reproduction = c("fixed", "dynamic"), effort = params@initial_effort, - include_resource = TRUE, + include_resource = FALSE, extinction_floor = 1e-06, h = 1e-04 ) diff --git a/man/get_steady_state_n.Rd b/man/get_steady_state_n.Rd index 419af68ae..a256bdc4d 100644 --- a/man/get_steady_state_n.Rd +++ b/man/get_steady_state_n.Rd @@ -55,8 +55,8 @@ only stalls; under-relaxation makes it converge.) The returned abundance is held at zero above each species' \code{w_max}, the same upper boundary condition that \code{\link[=project]{project()}} imposes (via \code{zero_above_support()} -in \code{project_n()}). Without this the bottom-up solve would carry density above -\code{w_max} whenever growth is still positive there or diffusion pushes density -past it. +in \code{project_n()}) and that \code{\link[=steadyNewton]{steadyNewton()}} solves on. Without this the +bottom-up solve would carry density above \code{w_max} whenever growth is still +positive there or diffusion pushes density past it. } \concept{helper} diff --git a/vignettes/cheatsheet-model-setup-and-calibration.Rmd b/vignettes/cheatsheet-model-setup-and-calibration.Rmd index ce4b0c79a..1d852a56e 100644 --- a/vignettes/cheatsheet-model-setup-and-calibration.Rmd +++ b/vignettes/cheatsheet-model-setup-and-calibration.Rmd @@ -194,6 +194,7 @@ params <- newTraitParams() # or newCommunityParams(), newSingleSpeciesP # ── Steady state ────────────────────────────────────────────────────────────── params <- steady(params) params <- steadySingleSpecies(params) # fast starting spectrum +params <- steadyNewton(params) # direct solve (experimental) # ── Calibrate to data (re-run steady() after each) ──────────────────────────── params <- calibrateBiomass(params) # total biomass → kappa diff --git a/vignettes/numerical_details.qmd b/vignettes/numerical_details.qmd index d48bee59e..0d886b372 100644 --- a/vignettes/numerical_details.qmd +++ b/vignettes/numerical_details.qmd @@ -204,7 +204,9 @@ size where growth has died away, so they are already zero and imposing @eq-upper-bc-zero changes nothing. With diffusion, density would otherwise leak to sizes above $w_{max}$; the boundary condition holds it at zero there instead. Because the location of $j_{max}$ is fixed by $w_{max}$ and the scheme — not by -the densities — it is the same at every time step and at the steady state. +the densities — it is the same at every time step and at the steady state, so the +steady-state solver ([steadyNewton()]) solves on exactly this support and its +solution is a fixed point of the dynamics. ## Numerical Diffusion @@ -619,6 +621,8 @@ be at most $1$ across the grid. If the Courant number exceeds 1, the lagged limi Because the reconstruction lives entirely in the coefficients $A,B,C$, the steady-state machinery stays consistent automatically: `getRequiredRDD()` reads the same boundary coefficients and `get_steady_state_n()` solves the same system. With the van Leer reconstruction the weight depends on the solution, so the steady state is found by an *under-relaxed* fixed-point iteration (freezing $\chi$ at the current iterate, solving, repeating), because at $\Delta t=1$ the operator is not diagonally dominant. Both `getRequiredRDD()` and `steadySingleSpecies()` read the same `second_order_w` slot, so a model is automatically set up at the steady state of exactly the scheme that `project()` will use, and that state is preserved to machine precision by all three time-stepping methods. +The direct solver [steadyNewton()] reads its set of unknowns from where the supplied abundances are non-zero, rather than from $w_{max}$. This keeps it robust when $w_{max}$ is set far above the largest fish (a common choice, so that the grid need not change when a parameter update produces larger fish): the structurally-zero classes below such a $w_{max}$ are simply excluded. The growth rate alone could not be used to find the support, because the main reason fish grow past $w_{repro\_max}$ is diffusion, whose rate only grows with $w$ and never vanishes; the abundance is the only reliable indicator of where the (possibly diffusion-fed) tail has died away. With the van Leer limiter the residual is only Lipschitz, so the Newton iteration converges to a fixed point of the dynamics but not to machine precision; the unlimited centred reconstruction — which admits an undamped odd-even mode at a steady state with no physical diffusion — gives an ill-conditioned steady-state Jacobian for which the direct solver is not expected to converge. + ## Steady-State Solution