SKILL.md
In mizer: Dynamic Multi-Species Size Spectrum Modelling

name: analyse-stability description: >- Analyse the dynamic stability of a mizer steady state and characterise limit cycles (experimental). Use whenever the user asks whether a steady state is stable or unstable, wants the leading eigenvalue, the period of an emergent oscillation, a Hopf bifurcation, the spectral radius of the numerical time step, a limit cycle to build or plot, or a bifurcation diagram over fishing effort — via getStability(), getDiscreteStability(), findSteadyState(), getOscillationModeSim() and scanModel(). This skill and calibrate-model share findSteadyState() and getSteadyResidual(): use calibrate-model to find a steady state, this skill to ask whether the state you found is stable. Assumes the standard semichemostat resource dynamics.

Analysing dynamic stability

Tools for asking whether a mizer steady state is dynamically stable — and, when it is not, characterising the limit cycle that replaces it. These are experimental: their interface may still change. findSteadyState(solver = "newton") assumes the standard semichemostat resource dynamics and solves for the resource alongside the fish, so the resource density and the feeding levels it implies are self-consistent even where consumers are satiated. getStability() and getDiscreteStability() perturb the resource alongside the fish and so work with any resource dynamics.

Distinct from calibration: the calibrate-model skill gets you onto a fixed point; what follows analyses the dynamics around it. The natural entry point is findSteadyState(params, solver = "newton"), which converges even when the fixed point is dynamically unstable, unlike the default solver = "project".

Is the steady state stable? — `getStability()`

getStability(params) computes the continuous-time eigenvalues of the model. Mizer discretises the size axis but not time, so on the size grid the model is a system of ODEs; getStability() differentiates their right-hand side directly and takes the eigenvalues of that Jacobian. No time step is involved, so the eigenvalues are a property of the model rather than of a solver. The result reports whether the state is stable or unstable, the maximum real part of the eigenvalues (max_real_part; > 0 means unstable), and the period at which the model rings (dominant_period for the dominant mode, oscillation_period for the leading oscillatory one). A complex pair is an oscillatory mode, not a bifurcation: its period is the period of an emerging limit cycle only where its real part is zero, which is a Hopf bifurcation and which a single spectrum cannot show. Watch Re(leading_oscillatory_eigenvalue) cross zero over a scanModel() before calling one. It also returns leading_eigenvectors, the top two eigenvectors of the state space split into $fish, a complex array (n_species, n_sizes, 2), and $resource, a complex matrix (n_w_full, 2), each normalised to maximum modulus 1 (the spatial shape of the oscillation).

params <- findSteadyState(params, solver = "newton")  # sit exactly on the (possibly unstable) fixed point
stab   <- getStability(params)
stab                                     # stable/unstable, growth rate, cycle period

effort sets the fishing effort used when forming the Jacobian.
h is the relative finite-difference step. Re-running with a different h is the cheapest check that the model is smooth enough for the analysis to mean anything: if the answer moves, do not trust it.

Both getStability() and getOscillationModeSim() linearise at the state stored in the object, so a model that is not on a fixed point gets eigenvalues for the neighbourhood of a point it is not sitting at. Both now warn when handed one; the fix is to run findSteadyState() (or tuneSteadyState()) first, not to ignore the warning. plot(getSteadyResidual(params)) shows how far off it is — see the calibrate-model skill.

tuneSteadyState() and findSteadyState() attach a related "convergence" attribute, which answers three questions separately: termination says why the run stopped, converged whether the solver met its own criterion, and attractor what the state reached actually is — "fixed_point", "limit_cycle" or NA. Read attractor: it is the only one set from the measured biomass drift (reported as residual), and so the only one that can be taken as a claim that the model is at a steady state.

Stopping on a fixed point takes both a distance function below distance_tol and a drift within residual_tol; a state that meets the first but not the second keeps running. Limit cycles are detected from the per-species biomass series sampled at every time step, at every check rather than only when the distance criterion has failed, because a cycle whose period divides t_check is sampled at one phase and looks converged. The relative-amplitude floor for calling an oscillation a cycle is the amplitude_tol argument (default 0.01).

getStability(), getDiscreteStability(), tuneSteadyState() and findSteadyState(solver = "newton") all treat components registered with setComponent() as fixed inputs: they are held at their stored values and have no rows in the Jacobian. The verdict is therefore about the consumer-resource subsystem with the components frozen, which is exact when a component is a fixed input and a good approximation when it is much faster or much slower than the fish. Mizer warns when it meets a component with dynamics of its own. findSteadyState(solver = "project") is the one that advances everything.

The reported residual measures that same subsystem: a component's rate of change is not folded into it, because mizer does not know what a component's state is measured in and so cannot form a biomass for it. So attractor == "fixed_point" is a claim about the consumers and the resource, not about the whole model. Mizer names any component that is moving in the same message, and attr(getSteadyResidual(params), "other") holds the rates themselves. To settle the components too, project the model with projectUntilSettled().

Is mizer's time step stable? — `getDiscreteStability()`

A different question, and the one to ask when a simulation disagrees with getStability(). getDiscreteStability(params, dt = 1) linearises the one-step map that project(method = "euler") takes at that step size and returns its discrete_eigenvalues and spectral_radius. Below 1 means the numerical scheme does not amplify perturbations at that dt. It is that step and not an approximation to it: the same transport coefficients, the same Thomas solver and the model's own resource_dynamics.

The two verdicts can differ, and that is the point: the implicit transport solve damps oscillations artificially, so an unstable steady state can have a spectral radius below 1 at a large dt and a simulation will then settle onto a state the model does not hold. If getStability() says unstable but the simulation goes flat, check getDiscreteStability() at the dt you are projecting with, and reduce dt (or use method = "second_order").

The discrete eigenvalues cannot be converted into continuous ones by any exact algebraic relation, because mizer's step is not fully implicit: the rates are evaluated at the state at the start of the step and only the transport solve is implicit. That is why getStability() differentiates the rates of change directly instead of inverting the one-step map.

Visualising the oscillation — `getOscillationModeSim()`

getOscillationModeSim(params) takes the output of findSteadyState() and builds a MizerSim covering exactly one period of the leading oscillatory mode in the linear approximation,

[ x(t) = x^* + A\,\mathrm{Re}!\left[e^{i\omega t}\,\mathbf v\right], ]

over the whole state (x = (N, n_{pp})), where (\mathbf v) is getStability()'s leading_oscillatory_eigenvector — the eigenvector of the dominant oscillatory mode, which is not in general the dominant mode — and (A) is scaled so that the largest relative swing in species biomass equals the amplitude argument (default 10%): the hardest-oscillating species departs that far from its steady biomass and no species departs further. Fish and resource share one (A), so the resource oscillates with the phase the mode gives it; on NS_params at effort 1.5 it leads the fish biomass by about 0.3 years of the 5.12-year period. The growth factor (e^{\sigma t}) is deliberately omitted so the oscillation closes after one period; the run ends exactly at (T), with a shortened final interval when t_save does not divide it, so the last state is the first state again. (\sigma) is recorded as growth_rate in the result's sim_params, and it is what decides whether this is a cycle the model settles onto or a transient it rings with on the way back.

A biomass amplitude does not bound the individual size classes — a cohort trough is a far larger relative excursion than the biomass integral over it — so getOscillationModeSim() reports when the perturbation drives an abundance negative and it has to be clipped at zero. On NS_params at effort 1.5 that starts at around amplitude = 0.2. Past that point the picture is no longer the linear mode. The result is an ordinary MizerSim, so plot it with the standard tools (see the analyse-and-plot skill):

params <- findSteadyState(params, solver = "newton")
sim    <- getOscillationModeSim(params, amplitude = 0.1)
plotBiomass(sim)                         # biomass oscillation over one period
animate(plotSpectra(sim))               # the travelling wave in the spectrum

This is the linear mode — the shape near onset. For the fully nonlinear cycle, project the real dynamics with the second-order scheme (see the run-simulation skill).

Where does the cycle appear? — `scanModel()`

scanModel() sweeps any aspect of the model and, at each value, follows the attractor of the full dynamics and records the long-term range of a summary quantity. Draw that range with style = "envelope" and you have a bifurcation diagram: a stable steady state shows as a single line, a limit cycle as a band between its minimum and maximum, so a Hopf bifurcation appears as the value at which the band opens up.

scan <- scanModel(params, scan_values = seq(0, 1.5, 0.05),
                  set_func = scanEffort(), value_func = getYield)
plot(scan, style = "envelope")

set_func says what to vary — scanEffort() for the fishing effort of every gear together, scanFishingMortality() for the fishing mortality on one species, scanSpeciesParam() for any species parameter. value_func says what to measure: any of getBiomass(), getYield(), getSSB(), getN(). species restricts which series are shown. See the analyse-and-plot skill.
The settling stage runs projectUntilSettled(), whose distance_tol, residual_tol, amplitude_tol and extinction_threshold are exposed for tuning when a run is slow to settle or a species is collapsing.
The returned MizerScan object is a data frame, so the numbers behind the diagram are available without re-running the scan.

getStability() involves no time step, so the temporal numerical diffusion of the implicit solver does not enter its answer at all. But the spatial numerical diffusion from the default first-order upwind scheme does: it is part of the semi-discretised model that the eigenvalues describe. A real limit cycle can still be damped to a flat line by the spatial scheme alone!

To accurately simulate the fully nonlinear oscillation and confirm a stable cycle, build the model with second_order_w = TRUE and project with method = "second_order". See the run-simulation skill.