inst/skills/upgrade-mizer-code/references/mizer-3.3.md

In mizer: Dynamic Multi-Species Size Spectrum Modelling

Upgrading from mizer 3.2 to 3.3

Most of the changes in this release are corrections, so for many models nothing moves at all. Results change for a model that had opted in to second-order bin-averaging or to the van_leer flux, that sets min_w below the default, that specifies sizes as lengths, that changes the resource power law after constructing the model, or that lets mizer fill in a gamma or f0 while carrying a hand-set search volume. Three interfaces change: the spectrum plots, the names of the steady-state finders, and the "convergence" attribute they attach. Everything else is a new report, and info_level = 0 silences the lot.

biomass and per_log_size replace power

plotSpectra(), plotSpectra2(), plotCDF(), plotCDF2() and animate() now describe the plotted quantity with two independent arguments: biomass chooses a biomass density rather than a number density, and the new per_log_size chooses a density with respect to logarithmic size rather than with respect to size. The power of the weight multiplying the number density is the sum of the two:

| | per_log_size = FALSE | per_log_size = TRUE | |---|---|---| | biomass = FALSE | power = 0 | power = 1 | | biomass = TRUE | power = 1 | power = 2 |

power still works and is still the only way to ask for a power that is not the sum of the two flags, so calls that pass only power are unaffected. Two things change:

• Passing power together with biomass used to ignore biomass silently. Now the two must agree, or you get an error (#501). Where they do agree the call is honoured: plotSpectra(sim, power = 1, biomass = FALSE) used to plot the biomass density, and now plots the number density with respect to logarithmic size — the same numbers, but labelled correctly and, with size_axis = "l", converted to a length axis with the logarithmic Jacobian. If you meant the biomass density, set biomass = TRUE. The same applies to plotlySpectra(), plotlyCDF(), plotlySpectra2() and plotlyCDF2(), which passed power on internally and so ignored biomass even when you gave only biomass: those calls now plot what they were asked for.
• plotCDF() and plotCDF2() do not accept per_log_size, because integrating a density over size gives the same cumulative quantity either way. Use biomass on its own there.

Arrays say what kind of value they hold

Mizer arrays now carry a type attribute saying what kind of quantity their values are: "value" (the default) for a rate or an amount, "density" for an amount per gram of body weight, "proportion" for a fraction. Two kinds of plotting behaviour follow from it, and both used to be decided some other way.

Densities. Plotting a density against a length axis (size_axis = "l") has to multiply the values by a Jacobian, because a density per gram is not a density per centimetre. mizer used to decide which arrays those were by looking at their metadata strings, treating an array as a density if it was named "Number density" or had units "1/g". For mizer's own number spectra — initialN(), N(), finalN(), NResource(), resource_capacity() — nothing changes; they were recognised before and are tagged now. What changes is getFluxGradient(): it is a rate of change of a number density, with units g^-1/year, and neither of the old string tests recognised it, so on a length axis its values were left as densities per gram and were mislabelled as such. They are now converted with the dw/dl = b w / l Jacobian and labelled cm^-1/year. The new curve is the right one; if you were reading values off the old one, they were per gram plotted against length.

Proportions. getFeedingLevel(), getCriticalFeedingLevel(), maturity(), repro_prop(), psi() and resource_level() now declare themselves proportions, and a plot of one shows the whole of the interval from 0 to 1 on a linear y axis, so the value can be read against the scale it belongs to. Three consequences:

• plot(getFeedingLevel(params)) and the other array plots gain that y range, where they used to fit the axis to the data. This is the range plotFeedingLevel() has always shown, so the dedicated function and the array plot now agree.
• plot(resource_level(params)) gets a linear y axis instead of a logarithmic one. Pass log_y = TRUE to get the old axis back; any explicit log_y or log you already pass is respected.
• The range is only ever widened to include the data, never narrowed to less than the interval from 0 to 1. So plotFeedingLevel(include_critical = TRUE) now shows a critical feeding level above 1, which the old fixed window drew off the top of the plot. Nothing is ever hidden, and an explicit ylim still wins.

Declaring it yourself. An array of your own is taken to be a density or a proportion only if you say so, by passing type to the array constructor. If you do not pass it, the old string tests still run as a fallback, so existing code that named an array "Number density" or gave it units "1/g" keeps working, and arrays saved by earlier versions keep working when they are loaded.

If you called the unexported plotting helpers directly, note that plotComparisonDataFrame() and the internal animate_plotly() take a single density_wrt argument in place of spectrum_power and spectrum_per_log_size, and that the internal array_spectrum_power() is gone. The power-based interface of plotSpectra() and friends is unchanged.

The total is summed on the axis it is plotted against

A total can only be formed once every line sits on the same coordinate. On a weight axis they do: every species shares the model's weight grid. On a length axis they do not, because each species — and now the resource — converts weight to length with its own allometric relationship, so at a given length the lines sit at different weights. That is why the Total line used to be dropped from length-based plots.

It is now summed after the conversion, at equal length rather than at equal weight, interpolating each line onto the union of all the size coordinates (logarithmically in size, with a line contributing nothing outside its own range). Where the lines already share a grid — always on a weight axis, and on a length axis whenever the weight-length parameters agree — the union is that grid and the interpolation reproduces the values exactly. The weight-axis total is unchanged, for every power.

plotSpectra2() and plotSpectraRelative() are fixed by the same change. They used to convert the size axis after assembling the two spectra, so the total they had been handed — already summed at equal weight — reached the conversion with no species to convert it by and was silently dropped. They now let plotSpectra() do the conversion, so the total they receive is the total on the axis being plotted.

That also settles what ylim does there. plotSpectra() applies it both as the axis limits and as a filter on the data, with a hard floor at 1e-20. plotSpectra2() could not do the same on a length axis, because the values it was filtering were a Jacobian away from the ones the limits described, so it skipped the filter. Now that it converts first, the filter applies as it does everywhere else. The plot is unchanged — the axis limits hid those points anyway — but return_data = TRUE no longer hands back values outside the limits.

One thing does change on the weight axis. total = TRUE now means the same thing everywhere: the total of everything the object holds, whatever is drawn. plotSpectra() always worked that way — the resource and every species count, whether or not the resource is shown and whichever species were selected — and it still does. The array plots did not: plot(<array>, total = TRUE) summed only the species selected for display, and only the sizes inside each species' own range. It now sums the whole array, so the total no longer moves when you change species, all.sizes or background, and a plot of two species can be read against the community total. If you were relying on the total of a selection, sum the selection yourself.

Comparison plots transform each array with its own model

plot2() and plotRelative() used to convert both of their operands together, with the MizerParams attached to the first one. That is invisible while the two models agree, and wrong as soon as they do not. The weight-length relationship w = a l^b is a species parameter, so with size_axis = "l" the second array was drawn at the first model's lengths, and a density was multiplied by the first model's Jacobian. Comparing a model against a version of itself with a changed a or b — exactly what these plots are for — therefore compared it against a curve in the wrong place.

Each array now prepares its own plotting data, with its own model, and the comparison receives two series that are already on the axis they will be drawn against. Two consequences, both only on a length axis and only when the models differ:

• The two lines sit where they belong, which for plot2() means the second one may move.
• plotRelative() interpolates. Once each series is converted with its own parameters the two length grids no longer coincide, and matching them by equality of the coordinate — which is what the old join did — kept only their exact coincidences, usually none. Each series is now interpolated linearly in the logarithm of size onto the union of the two sets of coordinates, restricted to the interval both cover so that nothing is extrapolated. Where the grids already coincide, which is always the case on a weight axis, the union is that grid and the interpolation reproduces the values exactly. plotSpectraRelative() works the same way.
• The relative difference of a density on a length axis moves, because the Jacobian of the two models no longer cancels out of the ratio. The old cancellation was only ever valid when both models shared a and b.

On a weight axis, and whenever the two models agree, every one of these plots is unchanged.

Two arrays that hold different kinds of value can no longer be compared at all. The type of an array decides whether its values are multiplied by a Jacobian on a length axis and whether the y axis covers the interval from 0 to 1, so a density and a plain value have no pair of axes in common; plot2() and plotRelative() now stop instead of warning and using the first array's type for both. The warnings about a differing value_name or differing units are unchanged.

plot2() and plotRelative() also take highlight now. They used to accept it into ... and discard it, so a call that asked for a thicker line got a plot that looked right but was not.

Background species and non-positive values in time-series plots

plot() on a time-by-species array — which is what plotBiomass(), plotYield() and plotN() draw — used to select the species and then append the background species as a second, separate step. Four things follow from doing it in one pass instead.

• A background species is drawn once. With species left at its default of all of them, every background species was selected under its own name and then appended again under the "Background" legend, so it was drawn twice, in two colours.
• background = FALSE removes them. It used to skip the appending step and leave the background species on the plot under their own names, which is the opposite of what it asks for.
• A background species is drawn when the selection asks for it. Previously every background species appeared whatever species said; now species decides, as it always has for plot() on a species-by-size array. Name a background species in species to see it — it keeps the "Background" legend.
• An empty background group no longer aborts the plot with replacement has 1 row, data has 0. That happened whenever no background value survived the filtering, for instance when a background species had died out.

Separately, these plots no longer discard values of zero or less when the y axis is linear. The floor at 1e-20 belongs to a logarithmic axis, where such values cannot be shown, and is still applied there. On a linear axis they are data like any other, and a quantity that can go negative used to lose exactly the part of it that was interesting. plotRelative() on a time-by-species array is drawn on a linear axis, so it keeps them too, which is what lets it show the -2 of a species that has gone from present to absent.

Animations follow the array type and the size axis

animate() now applies the axis handling the static plots have:

• A resource animation with size_axis = "l" labels its y axis 1/cm. It used to convert the values to a density per centimetre but keep the 1/g of the weight axis in the label, which named the curve wrongly.
• An animation of a "proportion" — getFeedingLevel(sim), or a resource level through time — is drawn on a linear y axis showing the whole of the interval from 0 to 1, as a static plot of one is. Pass log_y = TRUE or an explicit ylim to override it.

getProportionOfLargeFish() on a MizerParams object was wrong

The MizerParams method multiplied the species x size abundance array by the vector of weights, which R recycles down the columns of the array rather than along the size axis, so every species but the first was weighted by the wrong sizes. Only the MizerParams method was affected; the MizerSim method was always right, and the two now agree when applied to the same state (#494). Any Large Fish Index computed from a MizerParams object in a model with more than one species needs recomputing.

yield_observed belongs to the gear parameters

plotYieldObservedVsModel() now takes the observed yield from the yield_observed column of gear_params(), where the yield is given for each gear-species pair and the plot adds it up over the gears:

gear_params(params)["Cod, Otter", "yield_observed"] <- 3e11
plotYieldObservedVsModel(params)

Nothing breaks if your model keeps yield_observed among the species parameters: a species that has no observation in the gear parameters takes its value from there. What changes is that a model following mizer's own advice — given_species_params<-() has been telling you to use gear_params()<- — now works, where before the plot stopped with "You have not provided values for the column 'yield_observed'".

Length and weight parameters follow the one you gave last

A size can be given either as a weight (w_mat, w_max, …) or as the length it converts to (l_mat, l_max, …). Mizer used to derive the weight from the length whenever both were present, so on a model specified by lengths a weight could not be set at all: the value you assigned was replaced on the spot by the one calculated from the unchanged length.

Both now follow one rule: the one you gave last wins, and if you gave both at the same time the weight wins. The other is set to match, so the two never disagree, and mizer warns, naming the species, when it changes a length to match a weight it disagrees with.

params <- newMultispeciesParams(sp)   # sp specifies l_mat, a and b

# Used to be silently undone, now it takes effect and l_mat follows
species_params(params)$w_mat[1] <- 100

The rule is applied when a data frame is assigned into a model, which is when mizer can tell which values changed. A data frame you have taken out of a model and are editing on its own is left exactly as you write it — the conversions, checks and warnings happen on assignment. One that was never in a model, for example one passed to validSpeciesParams(), carries no such history, so a length and a weight that disagree there count as given at the same time and the weight wins.

species_params<-() and given_species_params<-() apply the rule identically (#490). Previously only species_params<-() did, so the same edit made through given_species_params<-() was discarded — a maturity weight differing by up to 73% — and the given species parameters were left permanently inconsistent, which made mizer repeat the "not consistent" warning at every later parameter change. If you worked around this by setting the weight and the length together, you can now set either one on its own.

Species parameter setters distinguish edits from declarations

The setters now express two different intentions:

| Setter | Meaning | |---|---| | species_params<-() | Edit the complete table. Mizer detects and records only entries whose values changed. | | given_species_params<-() | Declare the authoritative user input. Every non-NA entry is given, even when equal to the current calculated value; NA or a removed column hands it back to mizer. |

This makes it possible to protect a calculated value without changing the current model:

given_species_params(params)$q <- species_params(params)$q

The setters rebuild through setParams() only when the model can change. Provenance-only changes, observations, direct-runtime parameters and unrelated custom columns on a base MizerParams object are stored without rebuilding all rate arrays. Dependent parameters, demotions to calculated values, arguments of the active predation kernel and unknown columns on extension objects retain the conservative rebuild path. Call setParams() explicitly if the intention is to repair an object after direct slot manipulation.

given_species_params<-() also reports instructions that cannot take effect; species_params<-() stays quiet. It warns when a parameter is overruled by another given parameter, feeds a rate array set by hand, or belongs in gear_params(). Clearing an actually given value to NA counts as a change; adding an all-NA column does not. A frozen resource array is handled the same way. These warnings expose an existing no-op rather than changing the model: reset the named array to return it to parameter control, set the array directly, or use options(mizer_info_level = 0) when the warning is not wanted (#489).

Finally, defaults added while validating a species_params<-() assignment are no longer mistaken for values the user supplied. In particular, an unrelated edit no longer freezes newly filled a and b values. Such values can now move again when their inputs change and appear in calculated_species_params() rather than given_species_params(). Set a value explicitly if it should remain fixed (#496).

$ on a parameter table no longer partially matches

$ on a species_params or gear_params table now matches column names exactly. Partial matching was silently returning the wrong parameter: in a model without the length-weight parameters a and b,

species_params(NS_params)$a   # used to return the `alpha` column
species_params(NS_params)$b   # used to return the `beta` column

complete with per-species names, so code converting weights to lengths got the assimilation efficiency and the preferred predator/prey mass ratio instead. Writing was never partially matched (sp$b <- 3 always created a new column b), so reads and writes disagreed about what $b meant.

A name that is not a column now gives NULL. If the name would have partially matched a single column, you also get a warning naming that column. So is.null(species_params(params)$foo) is now a reliable test for whether a parameter is present, and any code that was relying on the abbreviation should spell the column out in full (#487).

w_min survives a rebuild of the species parameters

w_min is now part of given_species_params, so the min_w argument to newMultispeciesParams() and emptyParams() is preserved across any operation that rebuilds the species parameters. Previously a given_species_params<- round-trip silently reset w_min to 0.001 when min_w was smaller, and emitted a spurious warning when it was larger (#460). Code that set a small min_w and worked around the reset — or that unknowingly ran on the reset grid — now gets the size grid it asked for, and results change accordingly.

Defaults for gamma and f0 ignore a hand-set search volume

get_gamma_default() works out the search volume coefficient gamma that gives a species the target feeding level f0. It does that by giving the species a search volume coefficient of 1, putting a power-law resource spectrum in front of it and measuring the energy that becomes available. get_f0_default() is its inverse. A search volume you had set by hand used to leak into that measurement, and no longer does. This changes the species parameters of the models it affected, and the new values are the right ones.

The measurement used to build its unit search volume by calling setSearchVolume(), which refuses to recalculate a search_vol array you have frozen — so mizer's own internal call was blocked along with yours, and the available energy was measured with your array. The resulting gamma was wrong by whatever factor separated your array from the unit-gamma one, which in a realistic model is many orders of magnitude. Both functions now build the search volume they need directly from the species parameters (#488).

sv <- search_vol(params)
search_vol(params) <- sv * 10          # freeze the search volume

given_species_params(params)$gamma <- NA   # ask mizer to recalculate gamma

species_params(params)$gamma
#> Used to come back ~1e9 times too large; now the same value you would
#> get without the frozen search volume.

The recalculated gamma still has no effect on the model while the search volume stays frozen — mizer warns about that separately, see Species parameter setters distinguish edits from declarations above. Call setSearchVolume(params, reset = TRUE) to put the search volume back under the control of the species parameters.

f0 is always validated

Every non-missing target feeding level f0 must now be finite and in the interval [0, 1), whether or not a search-volume coefficient gamma is also supplied. Previously f0 = 1 divided by zero when mizer calculated gamma, silently creating an infinite gamma and a non-finite search_vol; values above 1 created negative search volumes. If gamma was supplied explicitly, the same invalid f0 could instead be accepted and ignored.

If code now errors here, choose a physically attainable feeding level below 1. When gamma is the parameter you intend to control, omit the f0 value or use a valid value; gamma will still take precedence (#517).

Resource scalars refresh calculated gamma and q

The resource power law is also the reference spectrum used to calculate search volume parameters. Changing lambda through resource_params<-() or setResource() now recalculates every q and gamma entry that mizer owns; changing kappa recalculates every mizer-owned gamma. A value you supplied explicitly remains protected, including when only some species in a column were given (#497).

Previously the resource capacity was rebuilt but the calculated species parameters and search_vol were left at the values for the old resource:

params <- newMultispeciesParams(sp)
resource_params(params)$lambda <- 2.2

# These now follow the new lambda automatically
species_params(params)$q
species_params(params)$gamma

If existing code deliberately wanted to keep the old values, record them as given before changing the resource:

given <- given_species_params(params)
given$q <- species_params(params)$q
given$gamma <- species_params(params)$gamma
given_species_params(params) <- given
resource_params(params)$lambda <- 2.2

setExtMort() warns when z0pre or z0exp is ignored

The z0pre and z0exp arguments of setExtMort() are used only to calculate values of the z0 species parameter that are not present in given_species_params(). If z0 is given for every species, calls such as

params <- setExtMort(params, z0pre = 2)
params <- setParams(params, z0exp = -0.25)

were accepted but changed nothing. They now warn that the arguments were ignored. reset = TRUE does not make them applicable: it hands the external-mortality array back to the species parameters but does not remove the given z0 values.

Set z0 explicitly when changing an existing model:

given_species_params(params)$z0 <- 2 * species_params(params)$w_inf^(-0.25)

The arguments still work wherever z0 is not given. A z0 value present only in species_params() is the cached result of an earlier calculation and is recalculated. If either argument was supplied explicitly, the newly calculated z0 values are recorded in given_species_params() so that they survive later rebuilds. Values calculated from the default z0pre = 0.6 and z0exp = n - 1 remain calculated parameters and are not recorded there (#493).

One report, one switch

Nearly everything mizer says while building or changing a model now goes through the same mechanism, including the reports in steady(), projectToSteady(), validParams(), setInteraction(), setReproduction(), setResource(), newTraitParams(), newSingleSpeciesParams() and plotYieldObservedVsModel(). Two consequences for existing code:

• info_level = 0 now means silence. Reports that were plain message() calls ignored info_level altogether and appeared anyway; they no longer do. If your code relied on seeing one of them, drop the info_level = 0.
• Reports are collected and given at the end of the call, one message and one warning rather than a stream. A test doing expect_message() on an individual report inside a longer call may need adjusting, and the text now arrives with any others in the same message.

The steady-state finders have new names

Neither steady() nor projectToSteady() said what distinguished them, and projectToSteady() returned a different class depending on an argument. Both are superseded:

| Superseded | Use instead | |---|---| | steady() | tuneSteadyState() | | projectToSteady() | findSteadyState(), or projectUntilSettled() for the trajectory |

The new names say what each one keeps. tuneSteadyState() holds the inputs to the fish dynamics — the reproduction rate and the resource abundance — at the values you supply while the spectra settle, and then adjusts the parameters that generate them, erepro/R_max and the resource capacity cc_pp, so that those held values are steady too. That is what steady() always did. findSteadyState() changes no parameter and lets reproduction, the resource and the spectra settle together, which is what projectToSteady() did.

projectUntilSettled() is the run itself and always returns a MizerSim; the two finders always return a MizerParams. So the return_sim argument is gone from the new functions:

# old
sim    <- projectToSteady(params, return_sim = TRUE)
params <- projectToSteady(params)
# new
sim    <- projectUntilSettled(params)
params <- findSteadyState(params)

It is called "settled" rather than "steady" because a run may equally settle on a limit cycle.

Both finders gained a solver argument. The default solver = "project" runs the dynamics until they settle, which is the old behaviour. solver = "newton" instead solves the steady-state equation directly with a Newton-type root finder from the nleqslv package, so it converges even when the steady state is dynamically unstable, where the time-stepping solver diverges away from it. findSteadyState(solver = "newton") carries the resource densities among its unknowns and so needs the default semichemostat resource dynamics; tuneSteadyState(solver = "newton") holds the resource fixed and works with any.

Three arguments are spelled differently on the new functions, because each of them did more than one job under the old name:

| Superseded | New | | |---|---|---| | tol | distance_tol | the run now also has a tolerance on the biomass drift, residual_tol | | t_per | t_check | how often the run checks whether it has settled; defaults to 15 * dt, so it can no longer contradict a dt you chose | | — | t_save | the interval at which the returned MizerSim is saved, as in project(); it is independent of t_check |

t_save is only on projectUntilSettled(), which is the only new function that returns a trajectory. The limit-cycle detection samples the biomass at every time step and has no argument of its own.

Nothing breaks: steady() and projectToSteady() are kept as thin wrappers that reproduce the old behaviour exactly, return_sim and t_per included. They do not warn, they will not be removed, and old code and old scripts keep running untouched. Renaming is a search and replace whenever you next touch the code, plus the three arguments above.

Because the wrappers do not warn, a user running old code sees no signal at all. Do not tell them their code is about to break; it is not. Suggest the new name when you are already editing the line.

The convergence attribute has a new shape

The "convergence" attribute attached by projectToSteady() and steady() used to carry type and settled. It now carries three fields in their place, because those two were being read as one answer to three different questions:

| Field | Answers | Values | |---|---|---| | termination | Why the run stopped | "residual_tolerance", "distance_tolerance", "cycle_detected", "time_limit", "extinction", and from the Newton solver "solver_converged", "solver_failed" | | converged | Whether the solver met its own criterion | TRUE/FALSE | | attractor | What the state reached is | "fixed_point", "limit_cycle", NA |

distance, years, period and amplitude are unchanged, and there is a new residual (the largest relative rate of biomass change at the state reached, in 1/year) and extinct (a character vector naming any species that went extinct during the run, or character(0) if none).

Code that tested conv$type == "steady" or conv$settled needs updating, and so does any expect_named() or names() check on the attribute. The translation:

# old
if (conv$settled && conv$type == "below_tolerance") ...
# new — the question is almost always about the state, not the run
if (identical(conv$attractor, "fixed_point")) ...

attractor is the only one of the three that may be read as a claim that the model is at a steady state: it is set from the measured biomass drift, so it is "fixed_point" only where that drift is within tolerance. converged says the numerics went well, which is a different thing and is why a limit cycle could previously be reported as a converged fixed point.

steady() and projectToSteady() still stop as soon as the distance function drops below tol. Their stopping rule has not changed, and a script that relies on a loose tol to get a quick, rough run still gets one. What is new is that they measure the model's biomass drift at the state they stop on, and record it. Where that drift is above 0.05/year the state is not a fixed point, and the attribute says so:

conv <- attr(steady(params, tol = 1e3), "convergence")
conv$termination   # "distance_tolerance" — that criterion, and only that, was met
conv$attractor     # NA — not a fixed point
conv$residual      # how fast it is still moving, in 1/year

The successful case is termination = "residual_tolerance" with attractor = "fixed_point". Nothing errors or warns in either case; code that ignores the attribute is unaffected.

The new tuneSteadyState() and findSteadyState() take the stricter line: a run there does not stop until the drift is within residual_tol as well (default 0.05/year), and carries on to t_max otherwise, reporting

#> Simulation run did not converge after 100 years. The distance function
#> returned 0.0002, which is below the distance tolerance, but the biomasses are
#> still changing at up to 0.4 per year, which is above the residual tolerance,
#> so this state is not a fixed point.

The remedy there is nearly always to look at what is moving, with plot(getSteadyResidual(params)), rather than to loosen residual_tol. Pass t_max = t_check to stop after a single check, or residual_tol = Inf for the wrappers' rule.

If a user reports that steady() has become slow, this is not the cause — its stopping rule is untouched. Check t_max, the model size, and whether they have switched to tuneSteadyState(), which is where the stricter rule lives.

steady() reports the tolerance it reached rather than announcing convergence

steady() and projectToSteady() used to end a successful run with

#> Convergence was achieved in 12 years.

They now say what was actually tested — the distance function dropped below tol, which is not the same thing as having reached a fixed point — and they report the biomass drift every time rather than only when it is large:

#> Reached the convergence tolerance after 12 years. The biomasses change at up
#> to 3.2e-05 per year.

Where the distance function is satisfied but the drift is not negligible, the message adds Reduce the tolerance on the distance function to converge further. It is a message, not a warning, because convergence at the tol you asked for did happen.

Code that matched the old wording needs the new string: an expect_message() in a test, or a script grepping the output. In a script the "convergence" attribute is the better check anyway, since info_level = 0 suppresses the message entirely.

The steady-state run advances time like project() does

steady() and projectToSteady() broke the run into blocks of t_per years, and each block used to start its clock again at t = 0. A rate function or a component's dynamics function that reads t therefore saw the same short interval over and over instead of a clock that runs. Seasonal forcing, a time-dependent external mortality, or anything else registered with setRateFunction() or setComponent() that depends on t was affected; the result differed from project() over the same period and could not settle onto a forced cycle.

The clock now runs from 0 to the end of the run, exactly as in project(). Results move for models with time-dependent dynamics, and the new numbers are the right ones. Models whose rates do not depend on t — which is nearly all of them — are unaffected to the last digit.

projectToSteady() ignores initial transients

To decide whether a simulation has settled onto a limit cycle, projectToSteady() calculates the autocorrelation of a fine-resolution biomass series. Previously it used the entire history from the start of the run. A large initial transient could therefore dominate the autocorrelation and obscure a cycle that had settled more recently. The autocorrelation step now uses only the second half of the series (or the most recent 20 samples if the series is shorter). A cycle will be found earlier, and some cycles that were previously missed entirely will now be correctly reported.

The steady-state tools hold other components fixed

tuneSteadyState() holds the components registered with setComponent() at their stored values while it solves for the spectra, and puts their dynamics back afterwards without solving them. The stability analyses likewise give the Jacobian a row for every fish cell and every resource cell and none for any component. That has always been true and was never said, so a model built with setComponent() could come back described as being at a fixed point of dynamics that had not been solved.

Mizer now warns when it meets a component whose dynamics_fun is not constant_other:

#> The component `detritus` has dynamics of their own, and mizer's steady-state
#> and stability machinery covers the consumers and the resource only: it is
#> held at the stored value throughout and is not included in the biomass drift
#> that mizer reports. See `attr(getSteadyResidual(params), "other")` for its
#> rate of change.

Nothing is refused: the analysis is the right one whenever the component is slaved to the fish or moves on a very different timescale, and only the user knows which. What has changed is that the assumption is now visible. Use findSteadyState(solver = "project"), or projectUntilSettled() for the trajectory, where the components need to settle too: the projection advances them like everything else and is not restricted.

A model with a custom resource_dynamics and no matching balance_<dynamics>() function gets a second warning from tuneSteadyState(), because the resource capacity that would make the preserved resource abundance steady cannot be derived:

#> There is no `balance_my_resource()` function, so the resource capacity could
#> not be rebalanced and the preserved resource abundance need not be a steady
#> state of `my_resource()`.

This too used to happen silently. Supply a balance_<dynamics>() function (see balance_resource_semichemostat() for the shape of one), or read the reported residual and decide whether the drift matters. resource_constant() is exempt: it hands back the abundance it was given, so there is nothing to rebalance.

Tests that assert that these calls are silent will see them fail. That is the point of the change; suppress with options(mizer_info_level = 0) where the assumption is deliberate.

summary() reports the steady state

summary() of a MizerParams object has a new block:

Steady state:
    biomass drift:  3.2e-05 /year   (at steady state)

Code that parses the output of summary() by line position needs updating. The same number is available directly as rowSums(getSteadyResidual(params)).

The match…() functions announce that they broke the steady state

matchBiomasses(), matchNumbers() and matchGrowth() now report that they have moved the model off its steady state. This is a message, so info_level = 0 or options(mizer_info_level = 0) silences it, as does the info_level argument, which matchGrowth() gains and which matchBiomasses() and matchNumbers() previously accepted but ignored.

The calibrate…() functions and scaleModel() say nothing, because they do not break the steady state: they apply one overall scaling factor, which is an exact symmetry of the model. If your workflow re-ran steady() after every calibrate…() step, that step was never necessary.

matchNumbers() also gains the empty-selection guard that matchBiomasses() already had. Its own guard could never fire, so when it had nothing to match — no number_observed values, or none for the species you asked for — it left the abundances alone but still called setBevertonHolt(), updated time_modified and announced that it had moved the model off its steady state. It now returns the model unchanged, as matchBiomasses() always did. Code that relied on the incidental re-tuning of the reproduction parameters should call setBevertonHolt() itself.

Fixes under the second-order size scheme

These all concern a model that has opted in to second-order bin-averaging or to the van_leer flux with second_order_w(). On the default first-order scheme nothing in this section changes anything. Under the second-order scheme, several functions were hand-rolling a first-order sum over the size grid instead of using the model's own quadrature, and are now consistent with it:

• getDiet(proportion = FALSE) was applying the prey-bin quadrature twice, making it uniformly too large by a factor (1 + beta) / 2, where beta is the grid ratio — 9.7% for NS_params. Summing the diet over prey now reproduces getEncounter() * (1 - getFeedingLevel()) (#474). getDiet(proportion = TRUE), the default, was unaffected: the factor was uniform and divided out.
• getTrophicLevel() built its numerator and denominator from different quadratures, so reported trophic levels were off by up to 0.06. A predator whose prey all have trophic level 1 now comes out at exactly 2 (#474).
• getN(params, min_w = ...) now bin-averages the size-range window, so the bin straddling min_w or max_w contributes only partially — as getBiomass() already did. Numbers over a restricted size range therefore change slightly; over the full size range nothing changes (#494).
• calibrateBiomass(), calibrateNumber(), matchNumbers(), plotBiomassObservedVsModel() and plotYieldObservedVsModel() each wrote out their own sum over the size grid and cut the size range at a bin boundary. The calibration functions therefore left the abundances at values that disagreed with the getBiomass() or getN() you would check them against, and a species matched to its observed biomass was then plotted off the 1:1 line. All five now integrate the same way getBiomass(), getN() and getYield() do, so a matched species really does come out at its observation (#504, #529). plotYieldObservedVsModel() is where the size of the error matters: its model yields were 10-20% below getYield(), and the total relative error in the plot caption is computed from them, so it told you the model under-predicted the yields when it did not. If you have read a yield calibration off that plot, re-read it.
• plot() on an ArrayTimeBySpeciesBySize took its time slice by hand and lost the array's representation tag, so a bin-averaged quantity such as getFMort(sim) or getFeedingLevel(sim) was drawn at the left bin edges instead of the geometric bin centres where it lives — one bin to the left of the right place. Only the drawn location moves; the values are unchanged.
• steady() under the van_leer flux used to fall into a limit cycle instead of converging: the flux limiter weights flipped from one cell to the next between iterations, and the iteration chased itself. The limiter is now relaxed with an exponential moving average, and the run converges (#522). A workaround — a steady() call wrapped in try(), a hand-set t_max, or a fall-back to the default upwind flux — is no longer needed. The steady state it now reaches is the one the van_leer discretisation actually has, so it differs from the upwind steady state the workaround was settling on; recalibrate rather than treat the difference as a regression.

If you have published absolute diet values, trophic levels or calibrated abundances computed under second_order_w, they need recomputing.

setParams() rejects arguments it does not use

setParams() passes its ... on to the rate setters, each of which declares its own ... as unused. Any argument that none of them recognises was therefore accepted and ignored without a word. It is now an error, and the error says where the argument belongs when it belongs somewhere:

setParams(params, metabolic = 99)        # was: silently ignored
setParams(params, resource_rate = 5)     # was: silently ignored

The resource case is the one most likely to have bitten: setParams() never called setResource(), so no resource argument ever reached the model, and the deprecation warnings for setResource(r_pp) and setResource(kappa) used to recommend setParams(resource_rate) and setParams(resource_capacity), which do nothing. Use setResource() for all of these:

params <- setResource(params, resource_rate = 5)

Likewise gear_params goes to gear_params<-(), and second_order_w and use_predation_diffusion to their own assignment functions. setResource() itself now applies the same check to its own ..., so a misspelled resource argument errors there too.

If your code errors here, the argument was having no effect before, so removing it changes nothing; moving it to the right function changes the model, and that is the change you had intended all along.

Two related tidy-ups: reset is now a documented argument of setParams() (it was already forwarded through ..., undocumented) and still thaws every rate array that setParams() sets; and setExtDiffusion() is now listed among the setters that setParams() calls, which it always did.

Diagnostic: if a user reports that a resource change "did not take", check whether they went through setParams(). Before this release the call was accepted, so there is no error in their logs and the model simply kept its old resource. resource_rate(params) before and after their call is the quickest confirmation.

One name for each stored rate array

Seventeen accessors that read a parameter or rate array back out of a MizerParams object had two interchangeable names. The bare name is now the one to use — it is the one that also has a replacement function, so the pair reads the same way in both directions (catchability(params) and catchability(params) <- value, reproduction_level(params) and reproduction_level(params) <- value). The get-prefixed names are superseded:

| Superseded | Use instead | |---|---| | getCatchability() | catchability() | | getSelectivity() | selectivity() | | getInitialEffort() | initial_effort() | | getInteraction() | interaction_matrix() | | getResourceDynamics() | resource_dynamics() | | getResourceLevel() | resource_level() | | getResourceRate() | resource_rate() | | getResourceCapacity() | resource_capacity() | | getPredKernel() | pred_kernel() | | getSearchVolume() | search_vol() | | getMaxIntakeRate() | intake_max() | | getMetabolicRate() | metab() | | getExtMort() | ext_mort() | | getExtEncounter() | ext_encounter() | | getMaturityProportion() | maturity() | | getReproductionProportion() | repro_prop() | | getReproductionLevel() | reproduction_level() |

Nothing breaks: the old names are kept as plain aliases of the new ones. They do not warn, they will not be removed, and they return exactly the same value, so old code and old scripts keep running untouched. Renaming is a search and replace whenever you next touch the code.

The get prefix now means one thing — a function that calculates something from the current state of a model, like getEncounter(), getFMort() or getBiomass(). The functions above only hand back a value that is already stored in the object.

Because the aliases no longer warn, a user running old code sees no signal at all. Do not tell them their code is about to break; it is not. Suggest the new name when you are already editing the line.

When a user's bin_average diagnostic disagrees with the rate functions, check whether they reached for pred_kernel() (formerly getPredKernel()) rather than encounter_kernel() — the rename does not change that distinction, but it makes the two names look more alike than they used to.

matchYields() and calibrateYield() have been removed

Both were deprecated in mizer 2.6.0 and nobody reported a use for them. They adjusted the abundance of a species to move its yield, which is the wrong lever: the yield is what the model predicts from the abundance and the fishing. Replace matchYields() with mizerExperimental::matchYield(), which adjusts the catchability instead:

# Old
params <- calibrateYield(params)
params <- matchYields(params)
# New
params <- mizerExperimental::matchYield(params)

calibrateYield() has no replacement. It rescaled the whole model so that the total yield summed over all species matched the total observation. If you were using it to set the scale of your model, use calibrateBiomass() with observed biomasses, or scaleModel() with a factor of your own choosing.

compareParams() compares small parameters properly

compareParams() now uses a relative tolerance for species parameters, so small-magnitude parameters such as gamma (~1e-8) are no longer treated as equal when they differ by up to ~10%. Comparisons that previously reported two models as identical may now report differences — those differences were always there.

The cheatsheet articles are now called guides

The topic articles that used to be called cheatsheets are called guides. A cheatsheet reminds you of something you already know; these articles assume no prior knowledge, so the name was wrong. Each article is now named after the agent skill it is generated from, and its title is that skill's own heading:

| Old article | New article | New title | |---|---|---| | cheatsheet-size-spectrum-dynamics | guide-understand-size-spectrum-dynamics | Guide: Understanding size-spectrum dynamics | | cheatsheet-model-setup | guide-build-model | Guide: Building a mizer model | | cheatsheet-calibration | guide-calibrate-model | Guide: Reaching steady state and calibrating | | cheatsheet-changing-parameters | guide-change-parameters | Guide: Changing model parameters | | cheatsheet-fishing | guide-set-up-fishing | Guide: Setting up fishing | | cheatsheet-running-simulations | guide-run-simulation | Guide: Running a mizer simulation | | cheatsheet-analysis-and-plotting | guide-analyse-and-plot | Guide: Analysing and plotting mizer results | | cheatsheet-stability | guide-analyse-stability | Guide: Analysing dynamic stability | | cheatsheet-extending-mizer | guide-extend-mizer | Guide: Extending mizer | | using-extension-packages | guide-use-extension-packages | Guide: Using mizer extension packages | | extending-mizer | guide-extend-mizer | Guide: Extending mizer |

On the website the old addresses redirect, so a bookmark or a link in your own writing still works. In R the old name does not resolve, because a vignette is looked up by exactly its file name:

# Old
vignette("cheatsheet-fishing")
# New
vignette("guide-set-up-fishing")

One of those rows is more than a rename. "Extending mizer" and "Guide: Extending mizer" were two articles on one topic; they are now the single guide-extend-mizer, holding both the article's worked examples and the guide's rules on quadrature schemes and discontinuous rates.

The build-multispecies-model skill was renamed to build-model in the same pass: it covers newTraitParams(), newCommunityParams() and newSingleSpeciesParams() as well, so its name claimed a narrower scope than it has. If you install mizer's skills with mizerAgents::setup_mizer_agent(), re-run it to pick up the new name.

mizer documentation built on Aug. 31, 2026, 5:08 p.m.