Skip to contents

This article collects the changes that may require you to update your own code or models when you upgrade to a new release of mizer. Existing model objects created with an earlier version continue to load and run — they are upgraded automatically — so the notes below are about changes in behaviour and in the functions you call, not about stored objects. The changes are grouped by the release in which they took effect, most recent first.

Only changes that can alter the behaviour of existing code are listed. The many purely additive features (new functions, new optional arguments, new plots) are described in the changelog and are not repeated here.


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, that lets mizer fill in a gamma or f0 while carrying a hand-set search volume, an extension or extra food, or that assigns a species parameter table holding only some of the model’s columns. 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 rather than the density one. If you meant the biomass density, drop the biomass argument. 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 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.

Extension packages that called the unexported plotting helpers directly should 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.

summary() of an array covers the same sizes as plot()

plot() of a species-by-size array draws each species over its own size range, from its w_min to its w_max. summary() of the same array reduced the whole size grid, so the two described different arrays. The values outside a species’ range describe an animal that does not exist, and because a rate usually grows with size they are also the extreme ones, so it was Min and Max that were reported wrongly:

summary(getEncounter(NS_params))$per_species[1, ]
#> before:  Species Sprat   Min 0.299   Mean 2929   Max 39573
#> now:     Species Sprat   Min 0.299   Mean 37.9   Max 240

40000 g/year is the encounter rate a 40 kg Sprat would have. 240 g/year is the rate over the sizes a Sprat reaches.

Both summary() methods with a size dimension — ArraySpeciesBySize and ArrayTimeBySpeciesBySize — now take all.sizes, defaulting to FALSE as plot() does. Pass all.sizes = TRUE for the whole grid:

summary(getEncounter(NS_params), all.sizes = TRUE)   # as before

A species with no values left in range now comes back as NA rather than as the -Inf/Inf and warning that min() and max() of an empty selection give.

print() and as.data.frame() are unchanged: they show the array as it is, without interpreting it. Nothing about the arrays themselves changed, so any code that indexes them directly is unaffected.

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.

An invalid w_mat25 is replaced by its default

validSpeciesParams() rejects a w_mat25 that is not smaller than w_mat. It used to say that it had corrected that by setting it to NA, which described a step the user never saw: setReproduction() fills the default in straight afterwards. It now says it is "marking it as missing so that its default will" be used.

Where the species parameters also carry l_mat25, this is a change in results and not only in wording. The length used to survive the correction and the length-to-weight conversion put the rejected value straight back, a rounding error below w_mat — small enough to pass the w_mat25 < w_mat assertion in setReproduction(), and large enough to collapse the maturity ogive to a step function. Such a model now gets the intended default, w_mat / 3^(1/10), and a smooth ogive. Maturity, and everything downstream of it, will differ (#580).

Mizer says when it defaults the weight-length parameters

A species parameter data frame with no a or b column has always been given the defaults a = 0.01 and b = 3, silently. Building a model from such a data frame now reports them among the other defaults it fills in:

i No `a` column so using a = 0.01 in w = a l^b, with w in g and l in cm.
i No `b` column so using the isometric default b = 3 in w = a l^b.

Nothing about the model changes; only the report is new, and it is at info_level = 3, so info_level = 0 silences it as it does the rest of the chatter.

It is worth reading rather than silencing. a = 0.01 with weights in grams and lengths in cm is Fulton’s condition factor K = 1, the textbook fusiform fish, and b = 3 is isometry. The exponent is a good default — it varies little across species — but a follows body shape over roughly two orders of magnitude, from around 0.001 for eel-like species to a few hundredths for deep-bodied ones, so for a specific species it is a placeholder rather than an estimate. On the twelve North Sea species shipped with mizer, whose fitted values span a = 0.001 to 0.010, the default gives lengths that are 1% to 42% short of the fitted relationship at a weight of 1 g.

That matters wherever a length reaches the model rather than a plot: the sigmoid_length, double_sigmoid_length and knife_edge_length selectivity functions convert l50, l25 and knife_edge_length from gear_params() into weights through a and b, so a defaulted relationship puts the selectivity curve at the wrong weights and changes the dynamics. The same applies to min_l and max_l in the summary and indicator functions, to getMeanLength(), and to plots drawn with size_axis = "l".

Because that case is the one that changes results, setFishing() reports it a second time, where the conversion actually happens:

i The gear selectivity for Sprat, Cod is set by length, but `a` and `b` were
  not supplied, so the conversion to weight used mizer's defaults (a = 0.01,
  b = 3). The selectivity therefore sits at weights that are unlikely to be the
  ones you intend. Supply the weight-length parameters in the species
  parameters.

This one is at info_level = 1, so it survives the setting that silences the routine chatter, and it is shown even when calc_selectivity() or setFishing() is called on its own rather than through setParams(). It names only the species whose selectivity is set from a length and whose a or b mizer had to fill in, and only the parameter that was missing — supply a and b for those species and it goes away. A gear that selects on weight (knife_edge, sigmoid_weight) never triggers it, whatever the species parameters say.

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).

A column dropped from an assigned species parameter table is removed

The species parameter setters also read the absence of a column as an instruction. A column that the table you assign does not have is one you no longer supply, so it is taken out of given_species_params(). Mizer then calculates afresh the parameters it knows how to calculate, and the ones it does not know about leave the model:

species_params(params)$gamma  <- NULL   # gamma is calculated again
species_params(params)$my_col <- NULL   # my_col is gone

given_species_params(params)$… <- NULL follows the same rule. Removing a column is reported at info_level 3, with a message beginning I have removed the species parameter column.

Previously neither route removed anything. Dropping a column from the table assigned to species_params<-() did nothing at all: the column was restored from the given species parameters, still recorded as given. Dropping one from the table assigned to given_species_params<-() did take it out of the given table, but left its value standing in species_params(), where calculated_species_params() reported the user’s own number as one mizer had calculated. Nor did the removal recalculate anything, so given_species_params(params)$gamma <- NULL left the previously given gamma in place instead of handing it back to mizer; it is now the same instruction as setting that entry to NA (#578).

Two things to watch for in existing code:

  • Assigning a table with only some of the model’s columns now withdraws the rest. This used to be a way of updating a few columns and leaving the others alone. Edit the table you get from species_params(params) instead of building a new one. The scope for surprise is limited, because species_params<-() has always validated what it is given: a table without species and one of w_inf, w_max or w_repro_max is an error, not a partial update.
  • A rebuild now happens where it did not before. Removing a column from the given species parameters goes through setParams(), so a parameter that was derived from the withdrawn one moves to its calculated value. That is the point of the change, but it means numbers can shift where previously nothing did.

For extension packages this is the supported way to withdraw a species parameter column when the user switches the extension off, replacing the params@species_params[["my_col"]] <- NULL slot manipulation that was the only thing that worked before.

$ 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.

The defaults for gamma and f0 are measured on mizer’s own reference state

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. Three separate things used to leak into that measurement, and none of them does now. All three change the species parameters of the models they affected, and the new values are the right ones.

A search volume you set by hand. 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.

The extension chain. The measurement used getEncounter(), which on an object carrying an extension marker class dispatches through the extension’s projectEncounter() method. So an extension that changes the encounter rate had its change folded into mizer’s own gamma — the one quantity the function goes out of its way to derive from the species parameters alone. Both functions now measure with mizerEncounter() (#577).

The consequence was not a one-off offset. gamma is a calculated species parameter, so it is recomputed on every rebuild, and the recomputed gamma is what the search volume is then built from — so the extension’s factor was re-applied each time: an extension that halves the search volume doubled gamma on every assignment to species_params(). Where an extension’s factor is zero for some species — a therMizer species sitting exactly on its temp_min, evaluated at the arbitrary t = 0 the default calculation uses — the available energy came out as zero and the calculation failed outright. An encounter function registered with setRateFunction() no longer enters these defaults either; if your model relied on that, supply gamma or f0 explicitly.

Additive encounter contributions. The ext_encounter array and every function registered with other_encounter(), including a component’s encounter_fun, are now excluded too (#586). They used to enter the call to mizerEncounter() used for the measurement, combined with a search volume coefficient of 1 for get_gamma_default() — where they were negligible beside the predation encounter — but with the species’ real gamma for get_f0_default(), where they counted at full strength. The two therefore stopped being inverses, and a model with extra food could report a substantially higher calculated f0 than the target its gamma had been derived from. A lower recalculated f0 is the intended reference-state value; use getFeedingLevel() to inspect feeding with the extra sources present.

When no default gamma can be calculated the error now names the species concerned and reports the energy measured for them. One model that used to build now raises it: a species with no predation encounter in the reference state at all — interaction_resource = 0, or a predation kernel that does not overlap the resource — but with an additive contribution to fall back on. The gamma it used to get was derived from that contribution alone and was many orders of magnitude away from a search volume coefficient, so the model it produced was not one to keep. Supply gamma explicitly for such a species.

If you are an extension author who worked around any of this by declaring the current gamma as a given species parameter, you can drop the workaround and let gamma follow f0 again:

given_species_params(params)$gamma <- NULL

The gamma that mizer then calculates describes the species’ baseline search volume, which is what a dynamic modulation such as a temperature scalar is meant to modulate.

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 check for a column name that looks like a typo of a standard species parameter moved onto this mechanism too, and at the same time stopped repeating itself (#581). It used to run inside every validation pass, so building a model reported the same column many times over — 31 warnings from one newMultispeciesParams() call — and every later edit of an unrelated parameter reported it again. It now runs once, where a column enters the model, and looks only at the columns the model does not already have. Two things follow: code that built models quietly with info_level = 0 no longer sees the warning at all, and a test wrapping model construction in expect_warning() with a count, or in suppressWarnings() to swallow a known run of duplicates, is now looking at one warning rather than many. The message itself is unchanged: "very close to standard parameter names".

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.

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.

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 an extension package’s tests, 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.

A size class holding no fish no longer blocks convergence

Above a size where the growth rate vanishes, a species’ density decays exponentially and dN/dt decays with it, so the density in the class falls through 1e-100 and beyond while never reaching zero. distanceSSLogN() counted every class with a positive density, so log(n) in such a class fell by the same amount between every pair of states and its contribution to the distance never shrank. One trace holding 3e-92 g of fish could therefore hold the distance above any tolerance indefinitely, and the run stopped at t_max reporting converged = FALSE (#570).

distanceSSLogN() now takes a biomass_share_cutoff, defaulting to 1e-8: a size class counts only if it holds at least that share of its species’ biomass. Relevance is measured as a share of biomass rather than of density, because density falls fifteen orders or more across a healthy spectrum for entirely good reasons, so no density threshold could tell a dying trace from real large fish.

Nothing changes for a model without such a trace. There, every class the cutoff removes is one that already had no density at all and was already excluded, so the number distanceSSLogN() returns is unchanged to the last bit and any distance_tol you have tuned keeps its meaning. What changes is the model that has one: a run that used to reach t_max now converges.

Pass biomass_share_cutoff = 0, to distanceSSLogN() directly or through projectUntilSettled(), for the old behaviour. Note that this is unrelated to the biomass_cutoff species parameter used by calibrateBiomass() and matchBiomasses(), which is a size in grams. What decides whether a state is a fixed point is untouched: residual_tol is measured against a biomass drift that integrates over every size class, cut off or not, so the cutoff cannot make a drifting model look settled.

A run that stops at t_max on a state that is nonetheless a fixed point now says so, instead of reporting only the distance:

#> Simulation run did not converge after 100 years. The distance function
#> returned 2.4, which is above the distance tolerance, but the state reached is
#> a fixed point: the biomasses change at only 0.0003 per year.

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.

An extension package whose tests 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.

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.

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.

Extension marker classes are created and repaired by mizer

Two fixes to the dynamic S4 marker classes that mizer 3.2 introduced, both for extension authors.

They used to be defined in .GlobalEnv, so anything that empties the global environment destroyed them: your own rm(list = ls()), and the cleanEx() that R CMD check runs before every example. R’s class cache went on reporting the class as present, so mizer’s repair did not fire and there was no way back short of restarting R. The next mizer call that had to resolve the class failed with base R’s error that the class is not defined — which for an extension package with examples that touch a params object meant a failing R CMD check. The marker classes now live in an environment that mizer attaches to the search path, mizer:extension-classes, which clearing the workspace does not reach (#587). If you have been working around the old behaviour by defining the marker classes by hand, as in

setClass("mizerShelf", contains = "MizerParams", where = globalenv())

you can drop that.

Separately, registerExtension() and registerExtensions() now rebuild the chain’s dynamic marker classes when a class went missing during an extension-package reload — the case devtools::load_all() creates, where the reload removes the S4 classes the package’s namespace held and the next coerceToExtensionClass() failed with base R’s coercion error. The whole chain is rebuilt rather than just the missing class, because R prunes a removed class from the contains list of its subclasses (#569). An intact chain is inspected and left untouched, so the usual repeated registration does no work, and the repair never installs, version-checks or loads anything.

other_mort() and other_encounter() register contributions that have no component

getMort() and getEncounter() add the result of every function listed in params@other_mort and params@other_encounter. Until now the only exported way to write into those lists was setComponent(), which needs a dynamics_fun and an initial_value, so an extension adding a term that depends on the model state but keeps no state of its own — a starvation or senescence mortality — had to assign into the slot directly. That still works, but the supported way is now

other_mort(params)[["starvation"]] <- "starvMort"
other_encounter(params)[["scavenging"]] <- "scavengingEncounter"

which checks that the name really is a function and that it does not collide with a component’s. Names must be unique, and assigning NULL to an entry removes it. Like other_params(), the new accessors show only what does not belong to a component, so assigning a whole list through the accessor can no longer wipe a component’s entry out (#579).

Two things can break existing code:

  • Encounter contributions are now passed the current simulation time as t, as mortality contributions already were. This makes a time-dependent extra encounter possible, but an existing component encounter_fun whose exact signature accepts neither t nor ... will now error. Add one of those arguments to the function.
  • setComponent() refuses a component name that a free-standing contribution is already registered under, where it used to take the entry over silently, leaving it live but invisible to the accessors. This can only affect code that both writes into @other_mort (or @other_encounter) by hand and then creates a component of the same name.

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
creating-extension-packages guide-create-extension-package Guide: Creating a mizer extension package
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")

Two of those rows are 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. Everything that only matters once you share an extension moved into guide-create-extension-package, so the two articles split along that line: the mechanisms for changing mizer’s dynamics in one, and turning that into a package other people can install in the other. Its advice on marker classes was corrected at the same time — it still told you to define them with setClass("myExtension", contains = "MizerParams"), which mizer 3.2 made unnecessary and which actively prevents your package from being chained with another, because a sealed class cannot be re-parented into the chain.

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.


Upgrading from mizer 3.1 to 3.2

species_params<-() now detects and protects changes

Previously, modifying species parameters via species_params<-() updated the values in the model but bypassed given_species_params(). This meant that your changes were not protected, and any subsequent recalculation of defaults (for example, by a call to given_species_params<-()) would overwrite your custom values. Furthermore, changing a parameter like w_inf via species_params<-() did not automatically trigger a recalculation of downstream parameters like w_mat or w_max.

Now, species_params<-() intelligently diffs the new data frame against the old one to detect exactly which parameters you have changed. It automatically records those changed parameters in given_species_params, protecting them from future overwrites, and immediately recalculates any downstream defaults based on your changes.

How this affects existing code:

  1. If your existing code used species_params<-() to update a core parameter like w_inf and you expected w_mat or w_max to remain frozen at their old values, you will now see them automatically recalculate. If you wish to freeze downstream parameters, you must provide their frozen values explicitly in the same update.

  2. If your code computes custom parameters and saves them via species_params<-(), those parameters will now be preserved and survive future recalculations.

Setting resource parameters

Two related changes affect how you modify the resource size spectrum. Together they make the resource scalars behave like the species parameters: a scalar is an input, and the size-dependent arrays are computed from it.

Assigning to resource_params() now updates the resource arrays

Previously, assigning to resource_params() — or to one of its components, such as resource_params(params)$kappa <- ... — only stored the new scalar values. The size-dependent carrying capacity (cc_pp) and replenishment rate (rr_pp) were left unchanged until you next called setResource().

Now these assignments immediately rebuild the arrays from the scalars, exactly as species_params()<- rebuilds the species rates:

  • kappa, lambda and w_pp_cutoff rebuild the carrying capacity;
  • r_pp and n rebuild the replenishment rate.

Arrays that you have set by hand are left untouched (see Frozen arrays below).

If your code changed a resource scalar and then called setResource() to apply it, nothing breaks — you can drop the now-redundant setResource() call. If you changed a resource scalar and relied on the arrays not changing until later, review that code.

Assigning to resource_params() does not balance the resource

Balancing means adjusting the rate and capacity together so that the resource replenishes at exactly the rate at which it is consumed, keeping it at its steady state. Assigning to resource_params() rebuilds the arrays from the scalars but does not balance, so the resource steady state generally shifts.

Balancing is now solely a feature of setResource(). To change a resource coefficient and keep the resource balanced, call setResource() rather than assigning to resource_params():

# Rebuild the capacity from a new coefficient and rebalance the rate,
# so the steady state is preserved:
params <- setResource(params, resource_capacity = new_kappa)

# Likewise, set a new rate coefficient and rebalance the capacity:
params <- setResource(params, resource_rate = new_r_pp)

The resource setters gained a balance argument

resource_rate<-, resource_capacity<-, resource_level<- and resource_dynamics<- still balance by default (unchanged behaviour), but they now accept a balance argument so you can switch balancing off:

# Set the capacity but leave the rate untouched (do not rebalance):
resource_capacity(params, balance = FALSE) <- my_capacity

Frozen arrays are protected from incidental balancing

When you set the size dependence of the resource capacity or the resource rate by hand (by assigning a full vector rather than a scalar), mizer marks it “set manually” — it is frozen and will not be recomputed from the resource parameters. Previously, an operation that re-balanced the resource without being given a replacement rate or capacity — for example changing only resource_dynamics, or calling setResource() with neither a rate nor a capacity — would silently overwrite such a frozen array. It is now kept, and a warning is issued instead.

To deliberately recompute a frozen array from the resource parameters, pass reset = TRUE to setResource().

The species_params data frame is now an S3 subclass

The species_params data frame now has class c("species_params", "data.frame") (and gear_params similarly). It behaves like an ordinary data frame, but subsetting and subassignment go through class-preserving S3 methods and can trigger reactive re-validation and conversions (for example filling in a weight from a length). Code that relied on class(species_params(params)) being exactly "data.frame", or that stripped attributes with the assumption of a plain data frame, may need adjusting. When you need a plain frame, coerce explicitly with as.data.frame().

Accessing a column with $ now returns a named vector

Extracting a single column from a species_params or gear_params object with $ now returns a vector named by species (or by "species, gear" for gear_params):

species_params(params)$w_mat
#>   Sprat  Herring      Cod
#>    ...      ...      ...

The values are unchanged, but the names are new. This is convenient for identifying entries, but code that compared such a vector with identical() to an unnamed vector, or that used it as-is where names matter (for example as row/column names elsewhere), may behave differently. Strip the names with unname() if you need the old behaviour. The species column itself is returned unnamed.

Setting sel_func adds the required argument columns

Assigning a selectivity function name to a gear_params object now automatically adds the argument columns that the function needs (as NA), ready to be filled in:

gp$sel_func <- "sigmoid_length"
# gp now has l25 and l50 columns, both NA

Previously these columns had to be added by hand. Code that checks which columns are present in gear_params, or that expected setting sel_func to leave the column set unchanged, will now see the extra columns (#431).

Passing a data frame to species_params() / given_species_params() now validates it

Calling species_params() or given_species_params() on a plain data frame now runs the same validation and defaults that validSpeciesParams() and validGivenSpeciesParams() apply, rather than only checking for misspellings and converting lengths to weights. species_params(df) fills in the default columns (w_max, alpha, n, p, interaction_resource, z_ext, and the rest), and given_species_params(df) applies the consistency corrections (for example clamping w_mat below w_inf), derives w_inf from w_max/w_repro_max when it is absent, and now stops if the frame has duplicate species rows. Models built or modified through newMultispeciesParams(), setParams() and the species_params()<- / given_species_params()<- setters are unaffected, because those already ran this validation. Only code that called the two accessors directly on a bare data frame will see the extra columns and stricter checks (#432).

Printing of mizer array objects shows the values

print() on the array objects returned by the rate getters (ArraySpeciesBySize, ArrayTimeBySpecies, ArrayResourceBySize, ArrayTimeByResourceBySize and ArrayTimeBySpeciesBySize, as returned by getEncounter(), getBiomass(), getFMort(), NResource() and similar) now truncates the output instead of flooding the console with all the array entries. If your code or reports relied on the old printed format, use as.data.frame() to go back to the full output.

Upper boundary condition at w_max

The size-spectrum solver now holds the abundance at zero above each species’ maximum size w_max. Without diffusion this happens automatically and results are unchanged. With diffusion switched on this change stops a small amount of density leaking to sizes above w_max, so results there change slightly. See vignette("numerical_details").

Extension packages: dynamic marker classes

If you develop a mizer extension, an installed extension package is now recognised as a dispatching extension from the S3 methods it registers for its marker class (for example getEncounter.mizerMR), rather than only from a statically defined S4 marker class. You can now omit the static setClass("mizerFoo", contains = "MizerParams") and let mizer create the marker class dynamically. This lets two independently developed extensions be chained in either load order. See the guide to creating a mizer extension package.


Upgrading from mizer 3.0 to 3.1

Version 3.1 leaves default results unchanged from 3.0 unless you opt in to the new experimental second-order-in-size scheme. The changes below can still affect existing code in specific situations.

Maximum-size species parameters clarified

The maximum-size parameters have been given clearer, separate roles (#325):

  • w_inf, the von Bertalanffy asymptotic size, is now the primary maximum-size parameter and is used as the default for w_repro_max (the size at which a mature individual invests all its energy in reproduction) and for w_mat.
  • w_max is now purely a computational boundary — it sets the size grid and the plot range — and defaults to 1.5 * w_inf.
  • The default external mortality parameter z0 is now computed from w_inf rather than w_max, so the computational boundary w_max no longer feeds into any model parameter.

Existing models and scripts are unaffected: if w_inf is not supplied it is taken from w_repro_max or w_max, so old objects behave as before. However, new models built from the defaults may differ from 3.0.0. If you build models from scratch, check that w_inf, w_max and w_repro_max mean what you intend.

getTrophicLevel() gives the resource a size-dependent trophic level

getTrophicLevel() and getTrophicLevelBySpecies() now assign the resource a size-dependent trophic level, \(T_R(w) = \max(1,\, 1 + \log(w / w_R) / \log(\beta_R))\), instead of treating the resource as trophic level 0. The new w_R and beta_R arguments control this. Trophic levels computed with these functions will therefore be higher than before. Set the arguments explicitly if you need to reproduce old numbers.

Bug fixes that change results

Several fixes correct earlier behaviour and so change output:

  • summary() of a MizerSim now reports the fishing effort that was used during the simulation, rather than the model’s initial_effort. Gears whose effort varied over time show the mean, flagged with a note giving the range. The printed summary therefore differs for simulations run with time-varying effort.
  • MizerSim method for plotDiet() introduced in version 3.0 simply plotted the diet at the initial time of the simulation. Now plotDiet() for a MizerSim accepts a time_range argument. The diet is now computed from the simulated abundances at the requested times, defaulting to the final saved step, rather than the initial one (#357).
  • Other components and t_save. project() was advancing the abundances of other components (set via setComponent()) only once per saved time step instead of once per dt step. They are now integrated with the same dt as the consumer and resource spectra, so results for models with other components no longer depend on t_save.
  • Time-varying effort in getRDI(), getRDD(), getFlux(). On a MizerSim object these now use the simulated time-varying effort rather than the initial effort, so they change for simulations with varying effort (#370).
  • plotCDF() / plotlyCDF() bin placement. Each cumulative value is now plotted at its bin’s upper edge, correcting a one-bin offset. The curves shift by one bin compared with 3.0 (#383).
  • distanceMaxRelRDI(). Now returns Inf instead of NaN when a previous RDI is zero, so projectToSteady() no longer mistakes a NaN distance for convergence. Convergence behaviour can therefore differ in edge cases.

Second-order methods advance the resource at the midpoint

If you use project() with method = "predictor_corrector" (or the new method = "tr_bdf2"), the resource and the other components are now advanced with midpoint rates rather than the start-of-step value, so that they reach the same second-order accuracy in time as the consumer spectra. Results from these methods therefore differ slightly from 3.0. The default method = "euler" and the steady states are unchanged.

Opting in to the second-order-in-size scheme

3.1 adds an optional, experimental second-order-accurate finite-volume scheme in the size variable, controlled by the new second_order_w slot. It is off by default, so default results are byte-identical to 3.0. If you switch it on (via second_order_w()<- or the second_order_w argument of the new...Params() constructors), size-integrated diagnostics and the resource spectrum shift by \(O(\Delta w)\), so a calibrated model may need recalibrating. See ?second_order_w and the “Numerical Details” vignette.


Upgrading from mizer 2.5.4 to 3.0

Version 3.0 is a large release. Most new capabilities are additive and off by default, but there are several renamed arguments, deprecations and behavioural changes that can affect existing code.

Renamed arguments and changed defaults (breaking changes)

  • First argument of plotBiomass(), plotYield(), plotYieldGear() (and their MizerSim methods and plotly* wrappers) is renamed from sim to object, for consistency with the other plot generics. Calls that passed the simulation by name, plotBiomass(sim = my_sim), must become plotBiomass(object = my_sim). Positional calls are unaffected.
  • plotBiomassObservedVsModel() / plotlyBiomassObservedVsModel() now default to ratio = FALSE for all object types. Calls that relied on the previous ratio plot must set ratio = TRUE explicitly.
  • plotDiet() no longer accepts a time_range argument. Remove it from your calls. (In 3.1 a time_range argument returns for the MizerSim method — see above.)
  • Dimnames of getMort() and getPredRate() arrays are now sp and w (matching getFMort() and the other rate getters). Code that referred to the old dimnames by name must be updated.

Rate getters return classed array objects

Functions that return arrays of the form (species × size), (time × species) or (time × species × size) now attach extra attributes and an S3 class (ArraySpeciesBySize, ArrayTimeBySpecies or ArrayTimeBySpeciesBySize). The numeric values and ordinary matrix behaviour (arithmetic, subsetting) are unchanged, but the extra class and attributes mean that a strict comparison such as identical(getMort(params), old_value) can now report a difference where the numbers agree. Use unclass(), or compare with all.equal() on the values, if you need to ignore the class. These objects also carry print(), summary(), plot() and as.data.frame() methods, so printing them looks different from a bare matrix.

setInitialValues() is deprecated

setInitialValues() is deprecated. Replace

params <- setInitialValues(params, sim)

with

params <- finalParams(sim)

or, when averaging over a time range, with getParams(sim, time_range, geometric_mean). This reflects a shift in interpretation: a MizerParams object now represents not just the model specification but also its current state (the abundances), which can be extracted from a simulation with getParams(), finalParams() and initialParams().

Growth can no longer be negative

Growth is now forced to be non-negative, preventing unphysical shrinkage. In any model where the energy available for growth used to go negative (for example a strongly food-limited large individual), growth is now clamped at zero instead, so projected size spectra can differ from 2.5.4. No warning is issued when growth stops at or after the maturity size.

project() timing and effort handling

  • Inherited dt and method. When project() is called on an existing MizerSim object, dt and method now default to the values stored in the simulation’s new sim_params slot. If you pass values that differ from the stored ones, a warning is issued. To use different settings deliberately, pass them explicitly and expect the warning.
  • t_max / t_save with an effort array. These arguments are now respected even when an effort array is supplied (#231). With t_max the simulation extends beyond the times in the effort array using the last known effort; with t_save the save frequency is controlled independently, interpolating effort as needed. Simulations that previously derived their length or save times solely from the effort array may now produce a different set of saved steps.
  • State at t_max always saved. project() now warns when t_max is not a multiple of t_save and ensures the state at t_max is saved even if the final interval is shorter than t_save (#341). The returned simulation may therefore contain one extra saved time step compared with 3.0.

plot() and summary() are now S3 methods

The plot() and summary() methods for MizerParams, MizerSim and the mizer array classes are now registered as S3 methods rather than S4 methods, so plot() and summary() stay plain S3 generics when mizer is loaded. This avoids interfering with S4 dispatch in other packages, but code that relied on plot/summary being S4 generics (for example via selectMethod() or getMethod()) needs adjusting.

Bug fixes that change results

  • getMeanMaxWeight() now applies the species selector to the denominator as well, so its values change when a subset of species is selected.
  • plotSpectra() axis limits. It no longer forces the y-axis lower limit to 1e-20 (it auto-scales to the data) and, when resource = FALSE, it uses min(params@w) rather than min(params@w) / 100 as the default lower size limit. Plots therefore look different.
  • getFMort() on a MizerSim was silently dropping the component names from n_other, breaking rate functions that access n_other by name (e.g. n_other[["resource"]]); it now preserves them.
  • getFMort.MizerSim() now passes the time argument t to user-defined fishing-mortality functions, so a time-dependent fishing function now sees the correct time.

Predation diffusion is available but off by default

3.0 adds a diffusion term to the growth dynamics, controlled by the new use_predation_diffusion slot. It defaults to FALSE, preserving the behaviour of earlier mizer, so existing models are unchanged unless you switch it on with use_predation_diffusion(params) <- TRUE. Likewise the new species parameters z_ext, d, E_ext and D_ext for external mortality, encounter and diffusion all default to values that leave the model unchanged.