diff --git a/DESCRIPTION b/DESCRIPTION index fe73257bd..8cdc7faca 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -30,6 +30,7 @@ Imports: grid, lubridate, methods, + nleqslv, plotly, plyr, progress, @@ -102,6 +103,7 @@ Collate: 'summary_methods.R' 'indicator_functions.R' 'plots.R' + 'plotBifurcation.R' 'plotBiomassObservedVsModel.R' 'plotYieldObservedVsModel.R' 'animateSpectra.R' @@ -122,6 +124,8 @@ Collate: 'match.R' 'matchGrowth.R' 'steadySingleSpecies.R' + 'steadyNewton.R' + 'getLimitCycleSim.R' 'defaults_edition.R' 'validSpeciesParams.R' Roxygen: list(markdown = TRUE) diff --git a/NAMESPACE b/NAMESPACE index 62a468fd0..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) @@ -382,6 +383,7 @@ export(getFluxGradient) export(getGrowthCurves) export(getInitialEffort) export(getInteraction) +export(getLimitCycleSim) export(getLinetypes) export(getM2) export(getM2Background) @@ -416,6 +418,7 @@ export(getSSB) export(getSearchVolume) export(getSelectivity) export(getSimParams) +export(getStability) export(getTimes) export(getTrophicLevel) export(getTrophicLevelBySpecies) @@ -481,6 +484,7 @@ export(noRDD) export(other_params) export(parsePlotLog) export(plot2) +export(plotBifurcation) export(plotBiomass) export(plotBiomassObservedVsModel) export(plotCDF) @@ -592,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/NEWS.md b/NEWS.md index 87fd7f139..5abaf4fc4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,98 +1,144 @@ -# 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. 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. + +- `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`. 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 + 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. 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. -- 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/getLimitCycleSim.R b/R/getLimitCycleSim.R new file mode 100644 index 000000000..f7c98bf01 --- /dev/null +++ b/R/getLimitCycleSim.R @@ -0,0 +1,155 @@ +# 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`. +#' @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, t_save = 0.1, ...) { + # ------------------------------------------------------------------ + # 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 + # ------------------------------------------------------------------ + 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 + 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 + + # ------------------------------------------------------------------ + # 3. Time grid + # ------------------------------------------------------------------ + n_steps <- ceiling(T_period / t_save) + t_seq <- seq(0, by = t_save, length.out = n_steps + 1) + + # ------------------------------------------------------------------ + # 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/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/plotBifurcation.R b/R/plotBifurcation.R new file mode 100644 index 000000000..103bc7a68 --- /dev/null +++ b/R/plotBifurcation.R @@ -0,0 +1,191 @@ +# 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 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 +#' 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, + tol = 0.01, + amplitude_tol = 0.01, + extinction_threshold = 1e-6, + 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, + 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, + 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, + tol = tol, amplitude_tol = amplitude_tol, + extinction_threshold = extinction_threshold, + 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/R/steady.R b/R/steady.R index fff46c24b..5df704046 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,7 +114,10 @@ projectToSteady <- function(params, t_per = 1.5, t_max = 100, 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, @@ -113,7 +131,10 @@ projectToSteady.MizerParams <- function(params, t_per = 1.5, t_max = 100, 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, @@ -123,10 +144,19 @@ 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") } + 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 +191,47 @@ 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`. + 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 @@ -182,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.") @@ -199,21 +259,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, amplitude_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 +314,100 @@ 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 `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, +#' 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 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, amplitude_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 <= 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) + # ... 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,8 +426,26 @@ 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 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. @@ -268,7 +462,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,8 +474,9 @@ projectToSteady.MizerParams <- function(params, #' params <- steady(params) #' plotSpectra(params) #' } -steady <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, - tol = 0.1 * dt, return_sim = FALSE, +steady <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, t_save = dt, + 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, @@ -288,7 +486,9 @@ 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, - tol = 0.1 * dt, return_sim = FALSE, + t_save = dt, + 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, @@ -322,11 +522,17 @@ 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, + amplitude_tol = amplitude_tol, + extinction_threshold = extinction_threshold, 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 +560,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/R/steadyNewton.R b/R/steadyNewton.R new file mode 100644 index 000000000..1845f893b --- /dev/null +++ b/R/steadyNewton.R @@ -0,0 +1,797 @@ +# 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) + + 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 + for (i in seq_len(nrow(N))) { + lo <- params@w_min_idx[i] + # 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 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") { + 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}. +#' 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. +#' +#' 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, 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()]. +#' @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`.} +#' \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 +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) + # ------------------------------------------------------------------------- + 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 + + 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 + } + + # ------------------------------------------------------------------------- + # 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, + leading_eigenvectors = leading_eigenvectors + ) +} 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/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/getLimitCycleSim.Rd b/man/getLimitCycleSim.Rd new file mode 100644 index 000000000..274917742 --- /dev/null +++ b/man/getLimitCycleSim.Rd @@ -0,0 +1,62 @@ +% 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, t_save = 0.1, ...) +} +\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}.} + +\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}.} +} +\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/man/getStability.Rd b/man/getStability.Rd new file mode 100644 index 000000000..1f0f282bc --- /dev/null +++ b/man/getStability.Rd @@ -0,0 +1,125 @@ +% 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}.} +\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{ +\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}. +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. + +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, 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{ +\code{\link[=steadyNewton]{steadyNewton()}} +} 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/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..3c346f439 --- /dev/null +++ b/man/plotBifurcation.Rd @@ -0,0 +1,143 @@ +% 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, + tol = 0.01, + amplitude_tol = 0.01, + extinction_threshold = 1e-06, + 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{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.} + +\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/man/projectToSteady.Rd b/man/projectToSteady.Rd index 641110e0f..5a2214dbb 100644 --- a/man/projectToSteady.Rd +++ b/man/projectToSteady.Rd @@ -11,7 +11,10 @@ projectToSteady( t_per = 1.5, t_max = 100, 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, @@ -42,9 +45,30 @@ 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.} +\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.} @@ -65,6 +89,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 +114,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..2f1c5baad 100644 --- a/man/steady.Rd +++ b/man/steady.Rd @@ -9,7 +9,10 @@ steady( t_max = 100, t_per = 1.5, 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, @@ -29,9 +32,30 @@ 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.} +\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.} @@ -54,7 +78,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/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..0632422f9 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 @@ -142,6 +143,9 @@ reference: walks through this calibration workflow. contents: - steady + - steadyNewton + - getStability + - getLimitCycleSim - steadySingleSpecies - matchGrowth - plotBiomassObservedVsModel @@ -229,6 +233,7 @@ reference: - plotFeedingLevel - plotYield - plotYieldGear + - plotBifurcation - plotFMort - plotMizerParams - plotMizerSim diff --git a/tests/testthat/test-cycle_detection.R b/tests/testthat/test-cycle_detection.R new file mode 100644 index 000000000..bb9af0288 --- /dev/null +++ b/tests/testthat/test-cycle_detection.R @@ -0,0 +1,131 @@ +# 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, 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. + 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, amplitude_tol = 0.01)) + # A flat series is not a cycle. + expect_null(detect_limit_cycle(matrix(1000, nrow = 200, ncol = 1), + t_save, amplitude_tol = 0.01)) + # Too few samples. + expect_null(detect_limit_cycle(matrix(1000, nrow = 10, ncol = 1), + 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 ----------------------------------------- + +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) +}) diff --git a/tests/testthat/test-getLimitCycleSim.R b/tests/testthat/test-getLimitCycleSim.R new file mode 100644 index 000000000..a05d1d94c --- /dev/null +++ b/tests/testthat/test-getLimitCycleSim.R @@ -0,0 +1,92 @@ +# 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 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) + + 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", { + 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, 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])) / + 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) +}) 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" + ) +}) diff --git a/tests/testthat/test-steadyNewton.R b/tests/testthat/test-steadyNewton.R new file mode 100644 index 000000000..50edb3f7a --- /dev/null +++ b/tests/testthat/test-steadyNewton.R @@ -0,0 +1,341 @@ +# 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 zeroes out the species + expect_warning(pn_ext <- steadyNewton(p_extinct, reproduction = "dynamic", extinction_floor = 1e-6), + "went extinct and were set to zero") + + # Verify Cod abundance is exactly 0 + lo <- p_extinct@w_min_idx[3] + expect_equal(pn_ext@initial_n[3, lo], 0) +}) + +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", + "leading_eigenvectors")) + 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", + "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) + 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) +}) + +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)) +}) 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/dynamic_stability.Rmd b/vignettes/dynamic_stability.Rmd new file mode 100644 index 000000000..c8be7e0c5 --- /dev/null +++ b/vignettes/dynamic_stability.Rmd @@ -0,0 +1,154 @@ +--- +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, 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 + +Mizer projects the state of the ecosystem forward in discrete time steps. We can write this discrete-time map abstractly as: + +$$ 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^*)$. + +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)$: + +$$ N_i(t) = N_i^* \left(1 + x_i(t)\right) $$ + +Linearising the dynamics for these small relative perturbations $x(t)$ gives: + +$$ 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 + +Let's load mizer and start with the standard North Sea model. + +```{r setup} +suppressMessages(devtools::load_all("..", quiet = TRUE)) +params <- NS_params +``` + +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. + +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. + +```{r find_steady} +params_f2 <- steadyNewton(params, effort = 2) +``` + +## Diagnosing the instability + +Is this steady state stable? We compute the spectral radius with `getStability()`. + +```{r default_stability} +getStability(params_f2, effort = 2)$spectral_radius +``` + +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: + +```{r coupled_stability} +stab <- getStability(params_f2, effort = 2, include_resource = TRUE) +stab$spectral_radius +stab$stable +``` + +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). + +```{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} = \frac{2\pi}{|\arg(\lambda)|} \text{ time steps}, $$ + +which mizer reports as `dominant_period`: + +```{r dominant_period} +stab$dominant_period +``` + +So the linear analysis predicts a growing oscillation with a period of roughly `r round(stab$dominant_period, 1)` years. + +## Watching the limit cycle emerge + +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. + +```{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) + +plotBiomass(sim_cycle) +``` + +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.) + +## 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. + +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") +``` + +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: + +```{r pts_steady} +attr(projectToSteady(params, t_max = 100), "convergence")$type +``` + +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. + +## Visualising the oscillation shape + +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. + +```{r limit_cycle} +lcs <- getLimitCycleSim(params_f2, include_resource = TRUE) +``` + +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: + +```{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. diff --git a/vignettes/numerical_details.qmd b/vignettes/numerical_details.qmd index 4c2f3563b..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 @@ -684,3 +688,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.