Package {mizer}

Title:	Dynamic Multi-Species Size Spectrum Modelling
Date:	2026-08-23
Type:	Package
Description:	A set of classes and methods to set up and run multi-species, trait based and community size spectrum ecological models, focused on the marine environment.
Maintainer:	Gustav Delius <gustav.delius@york.ac.uk>
Version:	3.3.0
License:	GPL-3
Imports:	assertthat, deSolve, dplyr, ggplot2 (≥ 3.4.0), ggrepel, grid, lubridate, methods, nleqslv, plotly, plyr, progress, Rcpp, reshape2, rlang, lifecycle, pak
LinkingTo:	Rcpp
Depends:	R (≥ 3.5)
Suggests:	testthat (≥ 3.0.0), withr, vdiffr, diffviewer, roxygen2, knitr, rmarkdown, quarto, pkgdown, covr, spelling
Collate:	'age_mat.R' 'helpers.R' 'info_signals.R' 'MizerParams-class.R' 'MizerSim-class.R' 'registerExtensions.R' 'ArraySpeciesBySize-class.R' 'ArrayTimeBySpecies-class.R' 'ArrayTimeBySpeciesBySize-class.R' 'ArrayResourceBySize-class.R' 'MizerScan-class.R' 'generic_methods.R' 'background.R' 'reproduction.R' 'saveParams.R' 'species_params.R' 'getRequiredRDD.R' 'setColours.R' 'setInteraction.R' 'setPredKernel.R' 'setSearchVolume.R' 'setMaxIntakeRate.R' 'setMetabolicRate.R' 'setMetadata.R' 'setExtMort.R' 'setExtEncounter.R' 'diffusion.R' 'second_order_w.R' 'setExtDiffusion.R' 'setReproduction.R' 'setResource.R' 'setFishing.R' 'setInitialValues.R' 'setBevertonHolt.R' 'upgrade.R' 'selectivity_funcs.R' 'pred_kernel_funcs.R' 'resource_dynamics.R' 'resource_semichemostat.R' 'resource_logistic.R' 'numerical_methods.R' 'transport.R' 'project_n.R' 'project.R' 'mizer-package.R' 'project_methods.R' 'rate_functions.R' 'sim_rates.R' 'sizeIntegral.R' 'summary_methods.R' 'indicator_functions.R' 'plots.R' 'plotBiomassObservedVsModel.R' 'plotYieldObservedVsModel.R' 'animateSpectra.R' 'newMultispeciesParams.R' 'wrapper_functions.R' 'newSingleSpeciesParams.R' 'steady.R' 'steadyState.R' 'scanModel.R' 'scan_setters.R' 'plotYieldVsF.R' 'extension.R' 'data.R' 'RcppExports.R' 'deprecated.R' 'get_initial_n.R' 'compareParams.R' 'customFunction.R' 'manipulate_species.R' 'observations.R' 'calibrate.R' 'scaleRates.R' 'match.R' 'matchGrowth.R' 'steadySingleSpecies.R' 'steadyNewton.R' 'getSteadyResidual.R' 'getOscillationModeSim.R' 'defaults_edition.R' 'validSpeciesParams.R' 'zzz.R'
Encoding:	UTF-8
LazyData:	true
URL:	https://sizespectrum.org/mizer/, https://github.com/sizespectrum/mizer
BugReports:	https://github.com/sizespectrum/mizer/issues
Language:	en-GB
RdMacros:	lifecycle
VignetteBuilder:	knitr, quarto
Config/testthat/edition:	3
Config/testthat/parallel:	true
Config/testthat/start-first:	steadyState, steadyNewton, getOscillationModeSim, rate_functions, steady, plots, manipulate_species
Config/roxygen2/version:	8.1.0
NeedsCompilation:	yes
Packaged:	2026-08-24 08:27:40 UTC; gustav
Author:	Gustav Delius [cre, aut, cph], Finlay Scott [aut, cph], Julia Blanchard [aut, cph], Ken Andersen [aut, cph], Richard Southwell [ctb, cph]
Repository:	CRAN
Date/Publication:	2026-08-24 08:50:02 UTC

mizer: Multi-species size-based modelling in R

Description

The mizer package implements multi-species size-based modelling in R. It has been designed for modelling marine ecosystems.

Details

Using mizer is relatively simple. There are three main stages:

1. Setting the model parameters. This is done by creating an object of class MizerParams. This includes model parameters such as the life history parameters of each species, and the range of the size spectrum. There are several setup functions that help to create a MizerParams objects for particular types of models:

   • newSingleSpeciesParams()

   • newCommunityParams()

   • newTraitParams()

   • newMultispeciesParams()

2. Running a simulation. This is done by calling the project() function with the model parameters. This produces an object of MizerSim that contains the results of the simulation.

3. Exploring results. After a simulation has been run, the results can be explored using a range of plotting_functions, summary_functions and indicator_functions.

See the mizer website for full details of the principles behind mizer and how the package can be used to perform size-based modelling.

Author(s)

Maintainer: Gustav Delius gustav.delius@york.ac.uk (ORCID) [copyright holder]

Authors:

• Gustav Delius gustav.delius@york.ac.uk (ORCID) [copyright holder]

• Finlay Scott drfinlayscott@gmail.com [copyright holder]

• Julia Blanchard julia.blanchard@utas.edu.au (ORCID) [copyright holder]

• Ken Andersen kha@aqua.dtu.dk (ORCID) [copyright holder]

Other contributors:

• Richard Southwell richard.southwell@york.ac.uk [contributor, copyright holder]

See Also

Useful links:

• https://sizespectrum.org/mizer/

• https://github.com/sizespectrum/mizer

• Report bugs at https://github.com/sizespectrum/mizer/issues

Check that a rate function returns the correct output dimensions

Description

Called by setRateFunction() to verify that a candidate rate function returns an array (or vector/list) of the correct dimensions for the requested rate.

Usage

.checkRateFunctionOutput(params, rate, fun)

Arguments

Value

Invisibly NULL. Called for its side-effect of stopping with an informative error if the output has the wrong shape.

S3 class for resource size spectra

Description

Several functions in mizer return a vector over the full size grid holding a resource-related quantity such as the resource number density, the resource mortality, the intrinsic resource birth rate or carrying capacity. The ArrayResourceBySize class wraps these vectors to provide convenient print(), summary(), plot(), and as.data.frame() methods.

Usage

ArrayResourceBySize(
  x,
  value_name = NULL,
  units = NULL,
  type = NULL,
  params = NULL
)

is.ArrayResourceBySize(x)

Arguments

Details

An ArrayResourceBySize object behaves just like a regular numeric vector for arithmetic operations and subsetting. It carries three lightweight attributes:

• value_name – a human-readable name for the value (e.g. "Resource mortality").

• units – the units of the value (e.g. "1/year").

• params – the MizerParams object that the value was computed from.

Value

An ArrayResourceBySize object (inherits from numeric).

is.ArrayResourceBySize() returns TRUE if x is an ArrayResourceBySize object, FALSE otherwise.

See Also

print(), summary(), as.data.frame(), plot(), plot2(), plotRelative(), addPlot()

Examples

mort <- getResourceMort(NS_params)
is.ArrayResourceBySize(mort)
summary(mort)
plot(mort)

The complete plotting data of a resource-by-size array

Description

The resource analogue of ArraySpeciesBySize_plot_data(): the weight limits, the conversion of the values and of the size coordinate onto the requested axis, and the length limits, all done with the array's own params. A resource array holds a single spectrum, so there is no selection, no background and no total to form.

Usage

ArrayResourceBySize_plot_data(
  x,
  wlim = c(NA, NA),
  llim = c(NA, NA),
  size_axis = "w",
  per_log_size = NULL
)

Arguments

Value

A data frame with the size coordinate in its first column, the values in its second, and Species and Legend columns.

S3 class for species x size rate arrays

Description

Many functions in mizer return two-dimensional arrays (species x size) holding rates like encounter rate, feeding level, growth rate, mortality etc. The ArraySpeciesBySize class wraps these arrays to provide convenient print(), summary(), plot(), and as.data.frame() methods.

Usage

ArraySpeciesBySize(
  x,
  value_name = NULL,
  units = NULL,
  type = NULL,
  params = NULL,
  representation = c("point", "average")
)

is.ArraySpeciesBySize(x)

Arguments

Details

An ArraySpeciesBySize object behaves just like a regular matrix for arithmetic operations and subsetting. It carries a few lightweight attributes:

• value_name – a human-readable name for the value (e.g. "Encounter rate").

• units – the units of the rate (e.g. "g/year").

• type – the kind of quantity the values are.

Value

An ArraySpeciesBySize object (inherits from matrix and array).

is.ArraySpeciesBySize() returns TRUE if x is an ArraySpeciesBySize object, FALSE otherwise.

See Also

print(), summary(), as.data.frame(), plot()

Examples

enc <- getEncounter(NS_params)
is.ArraySpeciesBySize(enc)
summary(enc)

The complete plotting data of a species-by-size array

Description

Everything a plot of an ArraySpeciesBySize needs, prepared once: the species selection, the background grouping, the masking of sizes outside a species' own range, the weight limits, the conversion of the values and of the size coordinate onto the requested axis, the total line, and the length limits.

Usage

ArraySpeciesBySize_plot_data(
  x,
  species = NULL,
  all.sizes = FALSE,
  wlim = c(NA, NA),
  llim = c(NA, NA),
  total = FALSE,
  background = TRUE,
  size_axis = "w",
  per_log_size = NULL
)

Arguments

Details

All of it uses the array's own params. That matters for the comparison plots, where the two operands may come from different models: a length axis and a density Jacobian are both built from the weight-length relationship of the model the values came from, so preparing the second array with the first one's parameters would put it in the wrong place on the axis. Each operand is therefore prepared here, on its own, and the comparison renderers receive data that is already on the axis it will be drawn against.

Value

A data frame with the size coordinate in its first column, the values in its second, and Species and Legend columns.

S3 class for time x resource-size arrays

Description

The NResource() function returns a two-dimensional array (time x size) holding the resource number density through time. The ArrayTimeByResourceBySize class wraps this array to provide convenient print(), summary(), plot(), and as.data.frame() methods.

Usage

ArrayTimeByResourceBySize(
  x,
  value_name = NULL,
  units = NULL,
  type = NULL,
  params = NULL
)

is.ArrayTimeByResourceBySize(x)

Arguments

Details

An ArrayTimeByResourceBySize object behaves just like a regular matrix for arithmetic operations and subsetting. It carries these lightweight attributes:

• value_name – a human-readable name for the value (e.g. "Number density").

• units – the units of the value (e.g. "1/g").

• params – the MizerParams object that the value was computed from.

Value

An ArrayTimeByResourceBySize object (inherits from matrix and array).

is.ArrayTimeByResourceBySize() returns TRUE if x is an ArrayTimeByResourceBySize object, FALSE otherwise.

See Also

print(), summary(), as.data.frame(), plot(), plot2(), plotRelative(), addPlot(), animate()

Examples

nr <- NResource(NS_sim)
is.ArrayTimeByResourceBySize(nr)
summary(nr)
plot(nr)

S3 class for time x species arrays

Description

Some functions in mizer return two-dimensional arrays (time x species) holding quantities like biomass, abundance, or yield rate through time. The ArrayTimeBySpecies class wraps these arrays to provide convenient print(), summary(), plot(), and as.data.frame() methods.

Usage

ArrayTimeBySpecies(
  x,
  value_name = NULL,
  units = NULL,
  type = NULL,
  params = NULL
)

is.ArrayTimeBySpecies(x)

Arguments

Details

An ArrayTimeBySpecies object behaves just like a regular matrix for arithmetic operations and subsetting. It carries these lightweight attributes:

• value_name – a human-readable name for the value (e.g. "Biomass").

• units – the units of the value (e.g. "g", "g/year").

• params – the MizerParams object that created the values.

Value

An ArrayTimeBySpecies object (inherits from matrix and array).

is.ArrayTimeBySpecies() returns TRUE if x is an ArrayTimeBySpecies object, FALSE otherwise.

See Also

print(), summary(), as.data.frame(), plot()

Examples

bio <- getBiomass(NS_sim)
is.ArrayTimeBySpecies(bio)
summary(bio)

S3 class for time x species x size arrays

Description

Some functions in mizer return three-dimensional arrays (time x species x size) holding quantities like fishing mortality, feeding level, or predation mortality through time. The ArrayTimeBySpeciesBySize class wraps these arrays to provide convenient print(), summary(), plot(), animate(), and as.data.frame() methods.

Usage

ArrayTimeBySpeciesBySize(
  x,
  value_name = NULL,
  units = NULL,
  type = NULL,
  params = NULL,
  representation = c("point", "average")
)

is.ArrayTimeBySpeciesBySize(x)

Arguments

Details

An ArrayTimeBySpeciesBySize object behaves just like a regular array for arithmetic operations and subsetting. It carries these lightweight attributes:

• value_name – a human-readable name for the value (e.g. "Fishing mortality").

• units – the units of the value (e.g. "1/year").

• type – the kind of quantity the values are.

• params – the MizerParams object that the value was computed from.

Value

An ArrayTimeBySpeciesBySize object (inherits from array).

is.ArrayTimeBySpeciesBySize() returns TRUE if x is an ArrayTimeBySpeciesBySize object, FALSE otherwise.

See Also

print(), summary(), as.data.frame(), plot(), animateSpectra()

Examples

fmort <- getFMort(NS_sim)
is.ArrayTimeBySpeciesBySize(fmort)
summary(fmort)
plot(fmort, time = 2007)

Beverton Holt function to calculate density-dependent reproduction rate

Description

Takes the density-independent rates R_{di} of egg production (as calculated by getRDI()) and returns reduced, density-dependent reproduction rates R_{dd} given as

R_{dd} = R_{di} \frac{R_{max}}{R_{di} + R_{max}}

where R_{max} are the maximum possible reproduction rates that must be specified in a column in the species parameter dataframe. (All quantities in the above equation are species-specific but we dropped the species index for simplicity.)

Usage

BevertonHoltRDD(rdi, species_params, ...)

Arguments

Details

This is only one example of a density-dependence. You can write your own function based on this example, returning different density-dependent reproduction rates. Three other examples provided are RickerRDD(), SheperdRDD(), noRDD() and constantRDD(). For more explanation see setReproduction().

Value

Vector of density-dependent reproduction rates.

See Also

Other functions calculating density-dependent reproduction rate: RickerRDD(), SheperdRDD(), constantEggRDI(), constantRDD(), noRDD()

Alias for set_multispecies_model()

Description

An alias provided for backward compatibility with mizer version <= 1.0

Usage

MizerParams(
  species_params,
  interaction = matrix(1, nrow = nrow(species_params), ncol = nrow(species_params)),
  min_w_pp = 1e-10,
  min_w = 0.001,
  max_w = NULL,
  no_w = 100,
  n = 2/3,
  q = 0.8,
  f0 = 0.6,
  kappa = 1e+11,
  lambda = 2 + q - n,
  r_pp = 10,
  ...
)

Arguments

Details

If species_params contains a w_inf column then it is copied to w_max. If max_w is not supplied then it is set to 1.1 * max(species_params$w_max). The supplied min_w_pp is shifted up by one grid step before being passed to newMultispeciesParams() to compensate for the fact that newer mizer versions extend the full size grid below min_w_pp.

Missing legacy columns in species_params are filled as follows: gear = species, k = 0, alpha = 0.6, erepro = 1, sel_func = "knife_edge", knife_edge_size = w_mat if needed, catchability = 1, ks = h * 0.2, and m = 1. If h is missing it is calculated from k_vb, alpha, f0 and w_max. If gamma is missing it is calculated from f0, h, beta, sigma, lambda and kappa.

Value

A MizerParams object

A class to hold the parameters for a size based model.

Description

Although it is possible to build a MizerParams object by hand it is not recommended and several constructors are available. Dynamic simulations are performed using project() function on objects of this class. As a user you should never need to access the slots inside a MizerParams object directly.

Details

The MizerParams class is fairly complex with a large number of slots, many of which are multidimensional arrays. The dimensions of these arrays is strictly enforced so that MizerParams objects are consistent in terms of number of species and number of size classes.

The MizerParams class does not hold any dynamic information, e.g. abundances or harvest effort through time. These are held in MizerSim objects.

Slots

metadata

A list with metadata information. See setMetadata().

mizer_version

The package version of mizer (as returned by packageVersion("mizer")) that created or upgraded the model.

extensions

Describes the extension chain needed to run the model. The entries are named by extension identifier (also the S4 marker class name) and ordered in S3 dispatch order, from outermost to innermost extension. It is either a named character vector whose values are requirement strings (version strings, installation specifications, or NA_character_), or a named list whose entries are length-2 character vectors c(requirement = ..., version = ...). The version records the version of the extension package that last upgraded the object (NA if unknown) and is used by needs_upgrading(). Use recordExtension() to write entries rather than modifying the slot directly. Extension subclasses are marker classes only and must not add slots.

time_created

A POSIXct date-time object with the creation time.

time_modified

A POSIXct date-time object with the last modified time.

w

The size grid for the fish part of the spectrum. An increasing vector of weights (in grams) running from the smallest egg size to the largest maximum size.

dw

The widths (in grams) of the size bins

w_full

The size grid for the full size range including the resource spectrum. An increasing vector of weights (in grams) running from the smallest resource size to the largest maximum size of fish. The last entries of the vector have to be equal to the content of the w slot.

dw_full

The width of the size bins for the full spectrum. The last entries have to be equal to the content of the dw slot.

w_min_idx

A vector holding the index of the weight of the egg size of each species

maturity

An array (species x size) that holds the proportion of individuals of each species at size that are mature. This enters in the calculation of the spawning stock biomass with getSSB(). Set with setReproduction().

psi

An array (species x size) that holds the allocation to reproduction for each species at size, \psi_i(w). Changed with setReproduction().

intake_max

An array (species x size) that holds the maximum intake for each species at size. Changed with setMaxIntakeRate().

search_vol

An array (species x size) that holds the search volume for each species at size. Changed with setSearchVolume().

metab

An array (species x size) that holds the metabolism for each species at size. Changed with setMetabolicRate().

mu_b

An array (species x size) that holds the external mortality rate \mu_{ext.i}(w). Changed with setExtMort().

ext_encounter

An array (species x size) that holds the external encounter rate E_{ext.i}(w). Changed with setExtEncounter().

ext_diffusion

An array (species x size) that holds the external rate at which the abundance density is redistributed over body size due to mixing, beyond the deterministic growth dynamics. Changed with ext_diffusion().

pred_kernel

An array (species x predator size x prey size) that holds the predation coefficient of each predator at size on each prey size. If this is NA then the following two slots will be used. Changed with setPredKernel().

ft_pred_kernel_e

An array (species x log of predator/prey size ratio) that holds the Fourier transform of the feeding kernel in a form appropriate for evaluating the encounter rate integral. If this is NA then the pred_kernel will be used to calculate the available energy integral. Changed with setPredKernel().

ft_pred_kernel_p

An array (species x log of predator/prey size ratio) that holds the Fourier transform of the feeding kernel in a form appropriate for evaluating the predation mortality integral. If this is NA then the pred_kernel will be used to calculate the integral. Changed with setPredKernel().

ft_pred_kernel_d

An array (species x log of predator/prey size ratio) that holds the Fourier transform of the feeding kernel in a form appropriate for evaluating the predation-diffusion integral (used when use_predation_diffusion is TRUE). It differs from ft_pred_kernel_e only when second_order_w[["bin_average"]] is TRUE, where it carries the extra power of prey size that the diffusion integrand needs; otherwise it equals ft_pred_kernel_e. Changed with setPredKernel().

rr_pp

A vector the same length as the w_full slot. The size specific growth rate of the resource spectrum.

cc_pp

A vector the same length as the w_full slot. The size specific carrying capacity of the resource spectrum.

resource_dynamics

Name of the function for projecting the resource abundance density by one timestep.

other_dynamics

A named list of functions for projecting the values of other dynamical components of the ecosystem that may be modelled by a mizer extensions you have installed. The names of the list entries are the names of those components.

other_encounter

A named list of functions for calculating the contribution to the encounter rate from each other dynamical component.

other_mort

A named list of functions for calculating the contribution to the mortality rate from each other dynamical components.

other_params

A list containing the parameters needed by any mizer extensions you may have installed to model other dynamical components of the ecosystem.

rates_funcs

A named list with the names of the functions that should be used to calculate the rates needed by project(). By default this will be set to the names of the built-in rate functions.

sc

The community abundance of the scaling community

species_params

A data.frame to hold the species specific parameters. See species_params() for details.

given_species_params

A data.frame to hold the species parameters that were given explicitly rather than obtained by default calculations.

gear_params

Data frame with parameters for gear selectivity. See setFishing() for details.

interaction

The species specific interaction matrix, \theta_{ij}. Changed with setInteraction().

selectivity

An array (gear x species x w) that holds the selectivity of each gear for species and size, S_{g,i,w}. Changed with setFishing().

catchability

An array (gear x species) that holds the catchability of each species by each gear, Q_{g,i}. Changed with setFishing().

initial_effort

A vector containing the initial fishing effort for each gear. Changed with setFishing().

initial_n

An array (species x size) that holds the initial abundance of each species at each weight.

initial_n_pp

A vector the same length as the w_full slot that describes the initial resource abundance at each weight.

initial_n_other

A list with the initial abundances of all other ecosystem components. Has length zero if there are no other components.

resource_params

List with parameters for resource.

A

Formerly used to flag background species via NA values. Replaced by the is_background column in species_params. Will be removed in a future version.

linecolour

A named vector of colour values, named by species. Used to give consistent colours in plots.

linetype

A named vector of linetypes, named by species. Used to give consistent line types in plots.

ft_mask

An array (species x w_full) with zeros for weights larger than the maximum weight of each species. Used to efficiently minimize wrap-around errors in Fourier transform calculations.

use_predation_diffusion

A logical flag controlling whether predation-induced diffusion is included when calculating rates with mizerDiffusion(). Defaults to FALSE to preserve the behaviour of previous mizer versions. Set to TRUE to enable the diffusion term from the jump-growth equation.

second_order_w

A named list with entry flux (the advective flux scheme: "upwind", "van_leer", or "centred") and logical entry bin_average (controls whether bin-averaging is used for rates). Both default to the first-order setting to preserve the behaviour of previous mizer versions.

See Also

project() MizerSim() emptyParams() newMultispeciesParams() newCommunityParams() newTraitParams()

S3 class for the result of a parameter scan

Description

scanModel() varies one aspect of a model over a range of values and measures a quantity on the attractor the model settles on at each of them. It returns its result as a MizerScan object, which is a data frame carrying, in addition, everything that plot() needs to draw it.

Usage

MizerScan(
  x,
  scan_name = NULL,
  scan_units = NULL,
  value_name = NULL,
  value_units = NULL,
  type = NULL,
  params = NULL,
  reference_lines = NULL,
  settings = NULL
)

is.MizerScan(x)

Arguments

Details

A MizerScan object behaves like an ordinary data frame with one row for each combination of scanned value and series. Its columns are, in order:

1

The scanned value. The column is named after the scan_name, so for example a scan over fishing effort has a column called "Fishing effort". Use names(scan)[[1]] rather than hard-coding it.

2

The measured quantity, averaged over the attractor. Named after the value_name.

Species

The series the row refers to. Named Species whatever the series are, because that is the column that mizer's colour and line-type machinery reads.

ymin, ymax

The smallest and largest value over the sampling window. On a fixed point these both equal the value; on a limit cycle they give the range of the oscillation.

attractor

What the state reached at this scan value is: "fixed_point", "limit_cycle" or NA for neither. This is the column that says whether the value in this row can be read as an equilibrium.

termination, converged

Why the run at this scan value stopped, and whether the solver met its own criterion. Both come from the "convergence" attribute that projectUntilSettled() attaches to its result, and neither is a claim about the state — see attractor for that.

period

The period of the limit cycle in years, or NA.

residual

How far the state still is from a fixed point, as a per-capita rate in 1/year, see getSteadyResidual().

The first three columns are the x, y and grouping variable in that order, which is the layout plotDataFrame() expects.

It also carries these attributes:

• scan_name, scan_units – name and units of the quantity that was varied, used for the x-axis label.

• value_name, value_units – name and units of the quantity that was measured, used for the y-axis label.

• type – the kind of quantity the values are, see array_types.

• params – the MizerParams object the scan started from, used for species colours and line types.

• reference_lines – an optional named numeric vector of positions on the x axis to mark with vertical lines, for example c(F_MSY = 0.32).

• at_max, max_value – for each series, the scanned value at which the measured quantity is largest, and the value it takes there. See the section below.

• settings – a list recording the settings the scan was run with.

Value

A MizerScan object, which inherits from data.frame.

is.MizerScan() returns TRUE if x is a MizerScan object, FALSE otherwise.

Where the maximum is

The at_max attribute holds, for each series, the scanned value at which that series' measured quantity is largest. On a yield-versus-fishing-mortality scan that is F_{MSY}; on a scan over fishing effort it is the effort that maximises the quantity being plotted. max_value holds the value attained there.

This is the largest value among those that were scanned, not the maximum of the underlying curve. It is therefore only as good as the grid you gave in scan_values, and the way to sharpen it is to scan a finer grid near the maximum, not to interpolate a coarse one. Subsetting a MizerScan with [ recomputes both attributes from the rows that remain, so they never go stale.

Limitations

Because the object is a data frame subclass, its attributes survive base R subsetting with [ but are dropped by functions that rebuild the data frame, including dplyr::filter(), dplyr::mutate(), subset() and transform(). A scan that has lost its attributes can no longer be plotted. Subset with [, or rebuild the object with MizerScan().

See Also

scanModel(), plot.MizerScan()

Other scan functions: plot.MizerScan(), plotYieldVsF(), scanEffort(), scanModel()

Examples

scan <- scanModel(NS_params, scan_values = c(0, 0.5, 1),
                  set_func = scanEffort(), species = "Cod")
scan
summary(scan)
attr(scan, "at_max")

Constructor for the MizerSim class

Description

A constructor for the MizerSim class. This is used by project() to create MizerSim objects of the right dimensions. It is not necessary for users to use this constructor.

Usage

MizerSim(params, t_dimnames = NA, t_max = 100, t_save = 1)

Arguments

Value

An object of type MizerSim

A class to hold the results of a simulation

Description

A class that holds the results of projecting a MizerParams object through time using project().

Details

A new MizerSim object can be created with the MizerSim() constructor, but you will never have to do that because the object is created automatically by project() when needed.

As a user you should never have to access the slots of a MizerSim object directly. Instead there are a range of functions to extract the information. N() and NResource() return arrays with the saved abundances of the species and the resource population at size respectively. getEffort() returns the fishing effort of each gear through time. getTimes() returns the vector of times at which simulation results were stored and idxFinalT() returns the index with which to access specifically the value at the final time in the arrays returned by the other functions. getParams() extracts the ecosystem state as a MizerParams object with initial abundances set to values from the simulation; finalParams() and initialParams() are convenient shorthands for the final and initial time steps. There are also several summary_functions and plotting_functions available to explore the contents of a MizerSim object.

The arrays all have named dimensions. The names of the time dimension denote the time in years. The names of the w dimension are weights in grams rounded to three significant figures. The names of the sp dimension are the same as the species name in the order specified in the species_params data frame. The names of the gear dimension are the names of the gears, in the same order as specified when setting up the MizerParams object.

Extensions of mizer can use the n_other slot to store the abundances of other ecosystem components and these extensions should provide their own functions for accessing that information.

The MizerSim class has changed since previous versions of mizer. To use a MizerSim object created by a previous version, you need to upgrade it with validSim().

Slots

params

An object of type MizerParams. If this params object uses extensions, the MizerSim object uses the same extension chain via params@extensions; MizerSim has no separate extension slot.

n

Three-dimensional array (time x species x size) that stores the projected community number densities.

n_pp

An array (time x size) that stores the projected resource number densities.

n_other

A list array (time x component) that stores the projected values for other ecosystem components.

effort

An array (time x gear) that stores the fishing effort by time and gear.

sim_params

A named list of the parameters passed to project() or projectUntilSettled() to produce this simulation, such as method and dt.

Time series of size spectra

Description

Fetch the simulation results for the size spectra over time.

Usage

N(sim)

NResource(sim)

Arguments

Value

For N(): An ArrayTimeBySpeciesBySize object (time x species x size) with the number density of consumers.

For NResource(): An array (time x size) with the number density of resource

Examples

str(N(NS_sim))
str(NResource(NS_sim))

Time series of other components

Description

Fetch the simulation results for other components over time.

Usage

NOther(sim)

finalNOther(sim)

Arguments

Value

For NOther: A list array indexed by time and component that stores the projected values for other ecosystem components.

For finalNOther: A named list holding the values of other ecosystem components at the end of the simulation

See Also

Other extension tools: clearExtensionChain(), coerceToExtensionClass(), getRegisteredExtensions(), initialNOther<-(), recordExtension(), registerExtension(), registerExtensions(), setComponent(), setRateFunction()

Example interaction matrix for the North Sea example

Description

The interaction coefficient between predator and prey species in the North Sea.

Usage

NS_interaction

Format

A 12 x 12 matrix.

Source

Blanchard et al.

Examples

params <- newMultispeciesParams(NS_species_params_gears,
                                interaction = NS_interaction)

Example MizerParams object for the North Sea example

Description

A MizerParams object created from the NS_species_params_gears species parameters and the inter interaction matrix together with an initial condition corresponding to the steady state obtained from fishing with an effort effort = c(Industrial = 0, Pelagic = 1, Beam = 0.5, Otter = 0.5).

Usage

NS_params

Format

A MizerParams object

Source

Blanchard et al.

See Also

Other example parameter objects: NS_sim

Examples

sim = project(NS_params, effort = c(Industrial = 0, Pelagic = 1, 
                                    Beam = 0.5, Otter = 0.5))
plot(sim)

Example MizerSim object for the North Sea example

Description

A MizerSim object containing a simulation with historical fishing mortalities from the North Sea, as created in the tutorial "A Multi-Species Model of the North Sea".

Usage

NS_sim

Format

A MizerSim object

Source

https://sizespectrum.org/mizer/articles/a_multispecies_model_of_the_north_sea.html

See Also

Other example parameter objects: NS_params

Examples

plotBiomass(NS_sim)

Example species parameter set based on the North Sea

Description

This data set is based on species in the North Sea (Blanchard et al.). It is a data.frame that contains all the necessary information to be used by the MizerParams() constructor. As there is no gear column, each species is assumed to be fished by a separate gear.

Usage

NS_species_params

Format

A data frame with 12 rows and 7 columns. Each row is a species.

species

Name of the species

w_max

Computational upper size boundary, defaulting to 1.5 * w_inf.

w_mat

Size at maturity

beta

Size preference ratio

sigma

Width of the size-preference

R_max

Maximum reproduction rate

k_vb

The von Bertalanffy k parameter

w_inf

The von Bertalanffy asymptotic size. This is the required maximum-size parameter.

Source

Blanchard et al.

Examples

params <- newMultispeciesParams(NS_species_params)

Example species parameter set based on the North Sea with different gears

Description

This data set is based on species in the North Sea (Blanchard et al.). It is similar to the data set NS_species_params except that this one has an additional column specifying the fishing gear that operates on each species.

Usage

NS_species_params_gears

Format

A data frame with 12 rows and 8 columns. Each row is a species.

species

Name of the species

w_max

Computational upper size boundary, defaulting to 1.5 * w_inf.

w_mat

Size at maturity

beta

Size preference ratio

sigma

Width of the size-preference

R_max

Maximum reproduction rate

k_vb

The von Bertalanffy k parameter

w_inf

The von Bertalanffy asymptotic size. This is the required maximum-size parameter.

gear

Name of the fishing gear

Source

Blanchard et al.

Examples

params <- MizerParams(NS_species_params_gears)

Ricker function to calculate density-dependent reproduction rate

Description

Takes the density-independent rates R_{di} of egg production and returns reduced, density-dependent rates R_{dd} given as

R_{dd} = R_{di} \exp(- b R_{di})

Usage

RickerRDD(rdi, species_params, ...)

Arguments

Value

Vector of density-dependent reproduction rates.

See Also

Other functions calculating density-dependent reproduction rate: BevertonHoltRDD(), SheperdRDD(), constantEggRDI(), constantRDD(), noRDD()

Sheperd function to calculate density-dependent reproduction rate

Description

Takes the density-independent rates R_{di} of egg production and returns reduced, density-dependent rates R_{dd} given as

R_{dd} = \frac{R_{di}}{1+(b\ R_{di})^c}

Usage

SheperdRDD(rdi, species_params, ...)

Arguments

Details

With b = 1/R_{max} and c = 1 this reduces to the Beverton-Holt reproduction rate, see BevertonHoltRDD().

Value

Vector of density-dependent reproduction rates.

See Also

Other functions calculating density-dependent reproduction rate: BevertonHoltRDD(), RickerRDD(), constantEggRDI(), constantRDD(), noRDD()

Add lines to an existing plot

Description

addPlot() adds another set of values to an existing ggplot, for example to compare the same rate before and after a model change. There are methods for all the mizer array classes. Each checks whether the existing plot uses a compatible x variable, and warns if the y variable or y-axis units appear to differ.

Usage

addPlot(
  plot,
  x,
  species = NULL,
  total = FALSE,
  background = TRUE,
  colour = NULL,
  linetype = "dashed",
  linewidth = 0.8,
  alpha = 1,
  ...
)

Arguments

Value

A ggplot2 object.

See Also

Other plotting functions: animate(), plot, plot2(), plotBiomass(), plotCDF(), plotCDF2(), plotDiet(), plotFMort(), plotFeedingLevel(), plotGrowthCurves(), plotMizerParams, plotMizerSim, plotPredMort(), plotRelative(), plotSpectra(), plotSpectra2(), plotSpectraRelative(), plotYield(), plotYieldGear(), plotYieldVsF(), plotting_functions

Examples

p <- plot(getEncounter(NS_params), species = "Cod")
addPlot(p, getEncounter(NS_params), species = "Cod")

pr <- plot(getResourceMort(NS_params))
addPlot(pr, getResourceMort(NS_params))

Add new species

Description

Takes a MizerParams object and adds additional species with given parameters to the ecosystem. It sets the initial values for these new species to their steady-state solution in the given initial state of the existing ecosystem. This will be close to the true steady state if the abundances of the new species are sufficiently low. Hence the abundances of the new species are set so that they are at most 1/100th of the resource power law. Their reproductive efficiencies are set so as to keep them at that low level.

Usage

addSpecies(
  params,
  species_params,
  gear_params = data.frame(),
  initial_effort,
  interaction,
  steady = TRUE,
  info_level = 3,
  ...
)

Arguments

Details

The resulting MizerParams object will use the same size grid where possible, but if one of the new species needs a larger range of w (either because a new species has an egg size smaller than those of existing species or a maximum size larger than those of existing species) then the grid will be expanded and all arrays will be enlarged accordingly.

If any of the rate arrays of the existing species had been set by the user to values other than those calculated as default from the species parameters, then these will be preserved. Only the rates for the new species will be calculated from their species parameters.

After adding the new species, the background species are not retuned and the system is not run to steady state. This could be done with tuneSteadyState(). The new species will have a reproduction level of 1/4, this can then be changed with reproduction_level<-().

Value

An object of type MizerParams

See Also

removeSpecies(), renameSpecies()

Examples

params <- newTraitParams()
species_params <- data.frame(
    species = "Mullet",
    w_max = 173,
    w_mat = 15,
    beta = 283,
    sigma = 1.8,
    h = 30,
    a = 0.0085,
    b = 3.11
)
params <- addSpecies(params, species_params)
plotSpectra(params)

Add the annotation layers to a MizerScan plot

Description

Add the annotation layers to a MizerScan plot

Usage

add_scan_annotations(
  p,
  x,
  plot_dat,
  reference_lines = TRUE,
  mark_max = FALSE,
  show_unsettled = TRUE
)

Arguments

Value

The ggplot object with the extra layers, still a mizer_plot.

Add a total line to plotting data by summing over its series

Description

The total has to be formed after the size coordinate has been converted, not before. On a weight axis every series shares the model's weight grid, so summing at equal weight and summing at equal position are the same thing. On a length axis they are not: each species, and the resource, converts weight to length with its own allometric relationship, so at a given length the series sit at different weights and their grids no longer coincide. The sum that means something there is the sum at equal length — the number of organisms per unit length, whatever they are — which is what this computes.

Usage

add_total_line(
  plot_dat,
  x_var = names(plot_dat)[[1]],
  value_col = 2,
  by = NULL
)

Arguments

Details

Each series is interpolated onto the sorted union of all the size coordinates, linearly in the logarithm of size, since the grid is logarithmic. A series contributes nothing outside its own range. When the series 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, so nothing is approximated in the cases where nothing needs to be.

Value

plot_dat with the total appended as a series named "Total".

Adjust the size grid

Description

This function adjusts the size grid in a MizerParams object to the desired minimum and maximum size. It can both expand and truncate the grid. If the grid is truncated, any data outside the new grid is discarded. A warning is issued if there is non-negligible biomass in the discarded size bins.

Usage

adjustSizeGrid(params, ...)

## S3 method for class 'MizerParams'
adjustSizeGrid(
  params,
  new_min_w = min(params@species_params$w_min),
  new_max_w = max(params@species_params$w_max),
  new_min_w_pp = min(params@w_full),
  preserve_species = params@species_params$species,
  tol = 1e-06,
  ...
)

Arguments

Value

A new MizerParams object with the updated size grid.

Calculate age at maturity

Description

Uses the size-dependent growth rate and the size at maturity to calculate the age at maturity.

Usage

age_mat(params, ...)

Arguments

Details

Using that by definition of the growth rate g(w) = dw/dt we have that

\mathrm{age_{mat}} = \int_0^{w_{mat}.}\frac{dw}{g(w)}

In the implementation this integral is approximated on the model size grid by summing dw / g(w) over all size bins with w < w_mat.

Value

A named vector. The names are the species names and the values are the ages at maturity.

Examples

age_mat(NS_params)

Calculate age at maturity from von Bertalanffy growth parameters

Description

This is not a good way to determine the age at maturity because the von Bertalanffy growth curve is not reliable for larvae and juveniles. However this was used in previous versions of mizer and is supplied for backwards compatibility.

Usage

age_mat_vB(object, ...)

Arguments

Details

Uses the age at maturity that is implied by the von Bertalanffy growth curve specified by the w_inf, k_vb, t0, a and b parameters in the species_params data frame.

If any of k_vb is missing for a species, the function returns NA for that species. Default values of b = 3 and t0 = 0 are used if these are missing. If w_inf is missing, w_max is used instead.

Value

A named vector. The names are the species names and the values are the ages at maturity.

Animate size-dependent quantities through time

Description

Creates an interactive plotly animation in which a play button steps through time, drawing one line per species at each frame.

Usage

animate(
  x,
  species = NULL,
  log_x = TRUE,
  log_y = TRUE,
  log = NULL,
  wlim = c(NA, NA),
  llim = c(NA, NA),
  ylim = c(NA, NA),
  tlim = c(NA, NA),
  size_axis = c("w", "l"),
  per_log_size = NULL,
  total = FALSE,
  background = TRUE,
  frame_duration = 500,
  transition_duration = frame_duration,
  easing = "linear",
  ...
)

animateSpectra(sim, ...)

Arguments

Details

The function dispatches on the class of x:

• MizerSim — animates the community abundance spectra (number density or biomass density vs body size). Resource, background species, and a community total can be added via the resource, background, and total arguments. The total is the total of everything the model holds, whatever is drawn. The biomass and per_log_size arguments choose the plotted quantity in the same way as in plotSpectra(), and the power argument is available as the same alternative to them. Both axes are log10 by default and can each be switched to linear with log_x = FALSE or log_y = FALSE.

• ArrayTimeBySpeciesBySize — animates any per-species, size-resolved quantity returned by a MizerSim accessor, such as getFMort(), getFeedingLevel(), or getPredMort(). Both axes are log10 by default and can each be switched to linear with log_x = FALSE or log_y = FALSE. Background species and a total can be added via the background and total arguments. The total is the total over every species the array holds, whatever is drawn.

• ArrayTimeByResourceBySize — animates the size-resolved resource quantity returned by NResource() on a MizerSim object. There is only a single resource spectrum, so species, total and background do nothing and warn if set.

Species linecolours and linetypes follow params@linecolour and params@linetype.

animateSpectra() is retained as a backward-compatible alias.

Value

A plotly object with one animated line trace per plotted group. Use the play button or the slider to step through time.

See Also

Other plotting functions: addPlot(), plot, plot2(), plotBiomass(), plotCDF(), plotCDF2(), plotDiet(), plotFMort(), plotFeedingLevel(), plotGrowthCurves(), plotMizerParams, plotMizerSim, plotPredMort(), plotRelative(), plotSpectra(), plotSpectra2(), plotSpectraRelative(), plotYield(), plotYieldGear(), plotYieldVsF(), plotting_functions

Examples

# Animate biomass density spectra, showing only sizes above 0.1 g
animate(NS_sim, power = 2, wlim = c(0.1, NA), tlim = c(1997, 2007))

# Animate fishing mortality through time
animate(getFMort(NS_sim))

# Animate feeding level for two species only
animate(getFeedingLevel(NS_sim), species = c("Cod", "Herring"))

# Animate the resource spectrum
animate(NResource(NS_sim))

Convert the contributors to a total and append the total line

Description

Wraps the two steps that every plot of an array with total = TRUE needs: the contributors are put on the same size axis as the data they will join, and then summed there, see add_total_line().

Usage

append_total_line(
  plot_dat,
  total_dat,
  params,
  size_axis,
  x,
  per_log_size = NULL
)

Arguments

Value

plot_dat with the total appended as a series named "Total".

Restrict plot data to a range of weights

Description

Internal helper that filters a plot data frame to the weight range given by wlim. It is exported so that extension packages (such as mizerMR) can reuse it in their own array plot() methods.

Usage

apply_wlim(data, wlim)

Arguments

Value

The subset of data with w inside wlim.

The density measure of a mizer array

Description

The bridge from the array metadata into the density machinery of the plots. Mizer arrays are indexed by the model's weight grid, so a stored density is always a density with respect to weight; the other measures in density_measures arise only for quantities that the spectrum plots compute on the fly, such as a density per logarithmic weight.

Usage

array_density_wrt(x)

Arguments

Value

"w" if the array holds a density, otherwise NA_character_.

The logarithmic y axis an array's type calls for

Description

A proportion belongs on a linear axis, so a plot of one turns log_y off unless the caller asked for a particular axis. Called before parsePlotLog(), which is why it also has to check log.

Usage

array_log_y(x, log_y, log, given)

Arguments

Value

The log_y to use.

The type of a mizer array

Description

The type of a mizer array

Usage

array_type(x)

Arguments

Value

One of array_types.

Kinds of quantity a mizer array can hold

Description

Mizer arrays record what kind of quantity their values are in their type attribute, because some kinds need handling that the numbers alone do not reveal:

"value"

the default: a rate, an amount, anything that needs no special handling.

"density"

an amount per gram of body weight, like a number density. Plotting a density against a length axis restates it per centimetre, which changes the values and not just the axis.

"proportion"

a fraction, like the feeding level. Plotted on a linear y axis showing the whole of the interval from 0 to 1, so that the value can be read against the scale it belongs to.

Usage

array_types

Format

A character vector of the three types.

Details

A "proportion" is not restricted to the interval from 0 to 1: the critical feeding level and the resource level can both exceed 1, and their plots show it. The type is a statement about what the number means, not a bound that mizer enforces.

Y-axis limits an array's type calls for

Description

Only a "proportion" has an opinion. A "density" is handled where the size axis is converted, and a "value" needs nothing.

Usage

array_ylim(x, ylim, log_y, values)

Arguments

Value

A numeric vector of length two.

Convert mizer arrays to data frames

Description

The as.data.frame() methods for mizer array classes turn matrix- and array-like results into tidy long-form data frames, with one row per observed combination of species, size and/or time. The numeric result is always stored in a column called value.

Usage

## S3 method for class 'ArraySpeciesBySize'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)
## S3 method for class 'ArrayTimeBySpecies'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)
## S3 method for class 'ArrayTimeBySpeciesBySize'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Details

The returned columns are:

• ArraySpeciesBySize: w, value, Species.

• ArrayTimeBySpecies: time, value, Species.

• ArrayTimeBySpeciesBySize: time, Species, w, value.

If the original object has non-numeric or missing dimension names, sequential indices are used for the time or w columns. Species names are taken from the row, column or dimension names of the original object.

Value

A data frame in long format.

See Also

print(), summary(), str(), plot(), ArraySpeciesBySize(), ArrayTimeBySpecies(), ArrayTimeBySpeciesBySize()

Examples

enc <- getEncounter(NS_params)
head(as.data.frame(enc))

biomass <- getBiomass(NS_sim)
head(as.data.frame(biomass))

Bring what a value function returned into a time by series matrix

Description

Bring what a value function returned into a time by series matrix

Usage

as_series_matrix(x, default_name = "Value")

Arguments

Value

A matrix with time in the rows and the series in the columns.

Assert that an object's extension chain is compatible with the session

Description

Stops with an informative error if the object's extension chain is not a suffix of the session's registered maximal chain, or (when check_class is TRUE) if the object does not inherit from the expected S4 marker class.

Usage

assertExtensionChain(
  object,
  extensions = objectExtensions(object),
  check_class = TRUE
)

Arguments

Value

Invisibly TRUE. Called for its side-effect of stopping on incompatibility.

Strip extension classes from a mizer object

Description

Coerces a MizerParams or MizerSim object back to its plain base class, removing any S4 extension marker classes. For MizerSim, also strips the extension class from the embedded params slot.

Usage

baseMizerClass(object)

Arguments

Value

The same object coerced to MizerParams or MizerSim.

Bin-average the weight of a size-spectrum integral

Description

Prepares the weight K(w) of an integral over the size spectrum so that the integral is evaluated with the quadrature scheme the model is actually using. Use this when writing your own indicator or diagnostic function; the built-in summary and indicator functions call it for you.

Usage

bin_average_weight(K, params)

Arguments

Details

An integral \int N(w) K(w)\, dw is discretised on mizer's finite-volume grid as \sum_j N_j \bar K_j \Delta w_j, where N_j is the cell average of the density over bin [w_j, w_{j+1}]. Only the weight is approximated: N_j is already a cell average and \Delta w_j is exact, so neither the abundance nor the bin widths should ever be passed through this function.

Whether the point weight K(w_j) is replaced by the bin average

\bar K_j = \frac{1}{\Delta w_j}\int_{w_j}^{w_{j+1}} K(w)\,dw \approx \tfrac12\big(K(w_j) + K(w_{j+1})\big)

is controlled by the bin_average entry of the model's second_order_w() slot. When it is FALSE (the default) K is returned unchanged, so an indicator written with this function reproduces the left-edge Riemann sums of previous mizer versions byte-for-byte. When it is TRUE the trapezoidal bin average is returned, which is uniformly second order and exact whenever K is linear in w (e.g. the first moment K = w, for which it equals (w_{j+1}^2 - w_j^2)/(2\Delta w_j)).

Because the gating happens inside, always call this rather than averaging unconditionally: a hard-coded bin average silently changes the results of models that are on the default scheme.

If K is a product of several size-dependent factors, average the product and not the individual factors — the average of a product is not the product of the averages. Spawning stock biomass, for example, averages maturity * w as a single weight.

The top bin has no right-hand neighbour on the grid, so its weight is left unaveraged (one-sided); the density there is negligible, so this does not affect the second-order accuracy of the totals.

Value

The weight K, bin-averaged when params@second_order_w[["bin_average"]] is TRUE, otherwise returned unchanged.

See Also

second_order_w(), get_size_range_array(), encounter_kernel()

Examples

# Biomass of each species above 10g -- what getBiomass() does internally.
params <- NS_params
K <- get_size_range_array(params, min_w = 10)   # species x size, 0/1
K <- sweep(K, 2, params@w, "*")                 # weight by w to get biomass
K <- bin_average_weight(K, params)              # gated on second_order_w
rowSums(sweep(initialN(params) * K, 2, params@dw, "*"))

Geometric bin centres of the size grid

Description

Internal helper for the second-order plotting code. A finite-volume cell average N_j = (1/\Delta w_j)\int_{w_j}^{w_{j+1}} N\,dw does not live at the left bin edge w_j but at the geometric bin centre

w^*_j = \sqrt{w_j\,w_{j+1}} = w_j\sqrt{\beta},

where \beta = w_{j+1}/w_j is the (constant) bin ratio of the logarithmic grid. This is the log-symmetric, second-order-correct location at which to plot a bin-averaged quantity (it is exact for the community spectrum N\propto w^{-2}). It is a uniform half-bin shift to the right on the log axis, the same for the consumer grid w and the full prey/resource grid w_full.

Usage

bin_midpoints(params, w_full = FALSE)

Arguments

Value

A numeric vector of geometric bin centres, one per grid node.

Box predation kernel

Description

A predation kernel where the predator/prey mass ratio is uniformly distributed on an interval.

Usage

box_pred_kernel(ppmr, ppmr_min, ppmr_max)

Arguments

Details

Writing the predator mass as w and the prey mass as w_p, the feeding kernel is 1 if w/w_p is between ppmr_min and ppmr_max inclusive and zero otherwise. ppmr_min must be strictly smaller than ppmr_max. The parameters need to be given in the species parameter dataframe in the columns ppmr_min and ppmr_max.

Value

A vector giving the value of the predation kernel at each of the predator/prey mass ratios in the ppmr argument.

See Also

setPredKernel()

Other predation kernel: gaussian_mixture_pred_kernel(), lognormal_pred_kernel(), power_law_pred_kernel(), truncated_lognormal_pred_kernel()

Examples

params <- NS_params
# Set all required paramters before changing kernel type
species_params(params)$ppmr_max <- 4000
species_params(params)$ppmr_min <- 200
species_params(params)$pred_kernel_type <- "box"
plot(w_full(params), pred_kernel(params)["Cod", 10, ], type="l", log="x")

Expand an array to a larger set of labelled dimensions

Description

Internal helper for sizeIntegral(). Replicates x along the dimensions it does not have and permutes its dimensions into the order given by target_labels, so that arrays with different dimensions can be multiplied together elementwise.

Usage

broadcast_dims(x, labels, target_labels, extent)

Arguments

Value

An array with dimensions extent[target_labels].

Calculate selectivity from gear parameters

Description

This function calculates the selectivity for each gear, species and size from the gear parameters. It is called by setFishing() when the selectivity is not set by the user. The returned array is initialised to zero, so gear-species combinations that are not listed in gear_params(params) remain zero. For each listed combination the function named in sel_func is called with w = params@w, the corresponding species parameters, and the selectivity parameters from the matching row in gear_params(params).

Usage

calc_selectivity(params)

Arguments

Value

An array (gear x species x size) with the selectivity values

Bin-averaged selectivity

By default the selectivity is point-sampled at the grid nodes params@w, i.e. at the left edge of each size bin. This is only first-order accurate in the bin size when the selectivity is used in the finite-volume update of the size spectrum. When the bin_average entry of the second_order_w() slot is TRUE, each selectivity function is instead integrated over its size bin, so that selectivity[g, i, j] holds the bin average

\bar S_{g,i,j} = \frac{1}{\Delta w_j} \int_{w_j}^{w_{j+1}} S_{g,i}(w)\, dw.

The integral is evaluated with a composite-midpoint rule on a log-spaced sub-grid of each bin, mirroring the bin-integrated predation kernel. This lifts the fishing mortality towards second order at no extra runtime cost (the integration happens once here, the rate functions are unchanged). A welcome side effect is that a knife-edge gear then gets the exact fraction of the straddling bin that lies above the knife edge, removing a grid artefact.

Examples

params <- NS_params
str(calc_selectivity(params))
calc_selectivity(params)["Pelagic", "Herring", ]

Calibrate the model scale to match total observed biomass

Description

Given a MizerParams object params for which biomass observations are available for at least some species via the biomass_observed column in the species_params data frame, this function returns an updated MizerParams object which is rescaled with scaleModel() so that the total biomass in the model agrees with the total observed biomass.

Usage

calibrateBiomass(params, ...)

Arguments

Details

Biomass observations usually only include individuals above a certain size. This size should be specified in a biomass_cutoff column of the species parameter data frame. If this is missing, it is assumed that all sizes are included in the observed biomass, i.e., it includes larval biomass.

After using this function the total biomass in the model will match the total biomass, summed over all species. However the biomasses of the individual species will not match observations yet, with some species having biomasses that are too high and others too low. So after this function you may want to use matchBiomasses(). This is described in the blog post at https://blog.mizer.sizespectrum.org/posts/2021-08-20-a-5-step-recipe-for-tuning-the-model-steady-state/.

Value

A MizerParams object. If no non-missing observed biomass values are provided, the original object is returned unchanged.

Examples

params <- NS_params
species_params(params)$biomass_observed <- 
    c(0.8, 61, 12, 35, 1.6, 20, 10, 7.6, 135, 60, 30, 78)
species_params(params)$biomass_cutoff <- 10
params2 <- calibrateBiomass(params)
plotBiomassObservedVsModel(params2)

Calibrate the model scale to match total observed number

Description

Given a MizerParams object params for which number observations are available for at least some species via the number_observed column in the species_params data frame, this function returns an updated MizerParams object which is rescaled with scaleModel() so that the total number in the model agrees with the total observed number.

Usage

calibrateNumber(params, ...)

Arguments

Details

Number observations usually only include individuals above a certain size. This size should be specified in a number_cutoff column of the species parameter data frame. If this is missing, it is assumed that all sizes are included in the observed number, i.e., it includes larval number.

After using this function the total number in the model will match the total number, summed over all species. However the numbers of the individual species will not match observations yet, with some species having numbers that are too high and others too low. So after this function you may want to use matchNumbers(). This is described in the blog post at https://blog.mizer.sizespectrum.org/posts/2021-08-20-a-5-step-recipe-for-tuning-the-model-steady-state/.

Value

A MizerParams object. If no non-missing observed number values are provided, the original object is returned unchanged.

Examples

params <- NS_params
species_params(params)$number_observed <-
    c(0.8, 61, 12, 35, 1.6, 20, 10, 7.6, 135, 60, 30, 78)
species_params(params)$number_cutoff <- 10
params2 <- calibrateNumber(params)

Calibrate the model scale to match a total observation

Description

Internal implementation shared by calibrateBiomass() and calibrateNumber(). Rescales the model with scaleModel() so that the total over all observed species of the modelled quantity agrees with the total of the observations. Species with no observation are left out of both totals.

Usage

calibrate_to(params, to = c("biomass", "number"))

Arguments

Value

A MizerParams object. If no non-missing observations are provided, the original object is returned unchanged.

Y-axis label for a cumulative-distribution plot

Description

Y-axis label for a cumulative-distribution plot

Usage

cdf_y_label(power, normalise)

Arguments

Value

A character string for the y-axis label.

Check that a data frame has the variables a band style needs

Description

Check that a data frame has the variables a band style needs

Usage

check_band_vars(var_names, style)

Arguments

Value

Nothing; called for the error.

Check that the rate arrays hold only finite values

Description

Gives an error if any of the arrays in the MizerParams object contains non-finite values, with the exception of the maximum intake rate, which is allowed to be infinite.

Usage

check_finite(params)

Arguments

Details

This check is cheap and is run on every call to validParams(), also when the repair work is skipped, because it catches exactly the kind of damage that the fingerprint in validation_key() does not cover: a bad value written into an array whose shape is unchanged.

Value

TRUE, invisibly.

Check that per_log_size applies to a mizer array

Description

Expressing values per logarithmic size only means anything for a density, so asking for it on anything else is an argument error rather than something to be quietly ignored — which is what ... used to do with it.

Usage

check_per_log_size(x, per_log_size)

Arguments

Value

per_log_size, invisibly, if it applies.

Clear the registered extension chain

Description

Clears the session's extension registry. You can then create a new extension chain with registerExtensions().

Usage

clearExtensionChain()

Value

Invisibly, an empty character vector.

See Also

The guide to using mizer extension packages

Other extension tools: NOther(), coerceToExtensionClass(), getRegisteredExtensions(), initialNOther<-(), recordExtension(), registerExtension(), registerExtensions(), setComponent(), setRateFunction()

Coerce a mizer object to its registered extension class

Description

Coerces a MizerParams or MizerSim object to the S4 marker class corresponding to the object's own extension chain. For MizerSim, the extension chain is read from sim@params@extensions.

Usage

coerceToExtensionClass(object, extensions = objectExtensions(object))

Arguments

Value

The same object coerced to the appropriate marker class, or to the base class for an empty extension chain.

See Also

"Creating a mizer extension package": Creating a mizer extension package

Other extension tools: NOther(), clearExtensionChain(), getRegisteredExtensions(), initialNOther<-(), recordExtension(), registerExtension(), registerExtensions(), setComponent(), setRateFunction()

Assemble the dimnames of a broadcast array

Description

Internal helper for sizeIntegral() that takes the dimnames of each labelled dimension from the first of the given arrays that has them.

Usage

collect_dimnames(target_labels, arrays, labels)

Arguments

Value

A named list of dimnames, or NULL if none of the arrays has any.

Compare two extension chains

Description

Compare two extension chains

Usage

compareExtensionChains(old, new)

Arguments

Value

One of "identical", "new_is_suffix", "old_is_suffix", or "incompatible".

Compare two MizerParams objects and print out differences

Description

Compare two MizerParams objects and print out differences

Usage

compareParams(params1, params2, ...)

Arguments

Value

Invisibly returns a character vector of difference messages, one element per difference. As a side effect, prints the differences in a human-readable format.

Examples

params1 <- NS_params
params2 <- params1
species_params(params2)$w_mat[1] <- 10
# Keep this example focused on the model parameter.
params2@time_modified <- params1@time_modified
compareParams(params1, params2)

Alias for validSpeciesParams()

Description

An alias provided for backward compatibility with mizer version <= 2.5.2

Usage

completeSpeciesParams(species_params)

Arguments

Details

validGivenSpeciesParams() checks the validity of the given species parameters. It throws an error if

• the species column does not exist or contains duplicates

• the asymptotic size w_inf is not specified for all species (but see the backwards-compatibility note below)

If a weight-based parameter is missing but the corresponding length-based parameter is given, as well as the a and b parameters for length-weight conversion, then the weight-based parameters are added. If both length and weight are given, then weight is used and an info_about_default condition is signalled if the two are inconsistent.

The required maximum-size parameter is w_inf, the von Bertalanffy asymptotic size of an average individual. For backwards compatibility, if no w_inf column is given, its values are taken from the w_repro_max column if that is present, or otherwise from the w_max column, and an informational message is issued. (w_repro_max is preferred over w_max because in earlier versions of mizer it was the size at which growth stopped and is therefore the closest analogue to the asymptotic size.)

Some inconsistencies in the size parameters are resolved as follows:

• Any w_mat that is not smaller than w_inf is set to w_inf / 4.

• Any w_mat25 that is not smaller than w_mat is set to NA.

• Any w_min that is not smaller than w_mat is set to 0.001 or w_mat /10, whichever is smaller.

• Any w_repro_max that is not larger than w_mat is set to 4 * w_mat.

The row names of the returned data frame will be the species names. If species_params was provided as a tibble it is converted back to an ordinary data frame.

The function tests for some typical misspellings of parameter names, like wrong capitalisation or missing underscores and issues a warning if it detects such a name.

validSpeciesParams() first calls validGivenSpeciesParams() but then goes further by adding default values for species parameters that were not provided. It only sets defaults for those species parameters that are not owned by a single rate-setting function, namely those that are read by several of them (n), that are used only when projecting (alpha), that determine the size grid (w_min, w_max), that are needed for the length-weight conversion (a, b) or that are used only for reporting (is_background). The function sets default values if any of the following species parameters are missing or NA:

• w_max is set to 1.5 * w_inf (it is only a computational boundary)

• w_repro_max is set to w_inf

• w_mat is set to w_inf/4

• w_min is set to 0.001

• alpha is set to 0.6

• n is set to 3/4

• a is set to 0.01

• b is set to 3

• is_background is set to FALSE

All other species parameters are given their default values by the rate-setting function that uses them, so that each default has a single home. For example p and k are set by setMetabolicRate(), z_ext, d and z0 by setExtMort(), E_ext by setExtEncounter(), D_ext by setExtDiffusion(), interaction_resource by setInteraction(), beta and sigma by setPredKernel(), q and gamma by setSearchVolume(), and erepro, m, w_mat25 and R_max by setReproduction(). These columns are therefore absent from the data frame returned by validSpeciesParams() but present in the species parameters of a MizerParams object, because setParams() calls all the rate-setting functions.

Note that the species parameters returned by these functions are not guaranteed to produce a viable model. More checks of the parameters are performed by the individual rate-setting functions (see setParams() for the list of these functions).

Value

For validSpeciesParams(): A valid species parameter data frame with additional parameters with default values.

For validGivenSpeciesParams(): A valid species parameter data frame without additional parameters.

See Also

species_params(), validGearParams(), validParams(), validSim()

Choose egg production to keep egg density constant

Description

The new egg production is set to compensate for the loss of individuals from the smallest size class through growth and mortality. The result should not be modified by density dependence, so this should be used together with the noRDD() function, see example.

Usage

constantEggRDI(params, n, e_growth, mort, diffusion, ...)

Arguments

Value

Vector with the value for each species

See Also

Other functions calculating density-dependent reproduction rate: BevertonHoltRDD(), RickerRDD(), SheperdRDD(), constantRDD(), noRDD()

Examples

# choose an example params object
params <- NS_params
# We set the reproduction rate functions
params <- setRateFunction(params, "RDI", "constantEggRDI")
params <- setRateFunction(params, "RDD", "noRDD")
# Now the egg density should stay fixed no matter how we fish
sim <- project(params, effort = 10, progress_bar = FALSE)
# To check that indeed the egg densities have not changed, we first construct
# the indices for addressing the egg densities
no_sp <- nrow(params@species_params)
idx <- (params@w_min_idx - 1) * no_sp + (1:no_sp)
# Now we can check equality between egg densities at the start and the end
all.equal(finalN(sim)[idx], initialN(params)[idx])

Give constant reproduction rate

Description

Simply returns the value from species_params$constant_reproduction.

Usage

constantRDD(rdi, species_params, ...)

Arguments

Value

Vector species_params$constant_reproduction

See Also

Other functions calculating density-dependent reproduction rate: BevertonHoltRDD(), RickerRDD(), SheperdRDD(), constantEggRDI(), noRDD()

Helper function to keep other components constant

Description

Helper function to keep other components constant

Usage

constant_other(params, n_other, component, ...)

Arguments

Value

The current value of the component

Restate the units of a density in a different density measure

Description

The per-size factor is found in either of the two spellings mizer uses for it, 1/g and g^-1, and is then swapped for the size unit of the target measure. A density per logarithmic size carries no size unit at all — a number per log weight interval is just a number — so converting to one removes the factor instead: 1/g becomes dimensionless and g^-1/year becomes 1/year.

Usage

convert_density_units(units, from, to)

Arguments

Details

Units that state no per-size factor are returned unchanged, since there is then nothing to identify as the size unit.

Value

The units expressed in the to measure.

Express plotting data on the requested size axis

Description

Converts the size coordinate of the plotting data to the requested axis and, when the values are a density, multiplies them by the Jacobian that restates them in the density measure that axis calls for (see density_target_measure()). Values that are not a density are left alone.

Usage

convert_plot_density_axis(
  plot_dat,
  params,
  size_axis,
  density_wrt = NA_character_,
  per_log_size = NULL,
  species_col = "Species",
  value_col = 2
)

Arguments

Details

Anything involving a length is a per-species quantity, because the weight-length relationship is, so rows whose species is not one of the model's species — the "Total" row, for instance — are dropped when a length is needed. Going from a density per size to one per logarithmic size needs no length, and keeps those rows.

Value

The plotting data with its size coordinate, and where called for its values, expressed for the requested axis. The size coordinate is the first column.

Convert plotting data from weight to length

Description

When size_axis = "l", adds a length column l computed from the weight column w using the weight-length relationship of each line, see plot_length_params(). Rows with no such relationship are dropped. For size_axis = "w" the data is returned unchanged.

Usage

convert_plot_size_axis(
  plot_dat,
  params,
  size_axis,
  species_col = "Species",
  drop_w = TRUE
)

Arguments

Value

The plotting data, with a length column l added (and w optionally dropped) when size_axis = "l".

Convert a weight-based spectrum to a length-based spectrum

Description

A density with respect to weight is converted to a density with respect to length with the Jacobian dw/dl = b * w / l. A density with respect to logarithmic weight is instead converted with ⁠d log(w) / d log(l) = b⁠. This is the interface used by the power-based spectrum plots; arrays carry their density measure explicitly and use convert_plot_density_axis() instead.

Usage

convert_plot_spectrum_axis(
  plot_dat,
  params,
  size_axis,
  power,
  per_log_size = power == 2,
  species_col = "Species",
  value_col = 2
)

Arguments

Value

The plotting data with both its size coordinate and spectrum values expressed for the requested axis.

Copy the metadata of a MizerScan onto another object

Description

Copy the metadata of a MizerScan onto another object

Usage

copy_scan_attributes(to, from)

Arguments

Value

to with the attributes and class of a MizerScan.

Replace a mizer function with a custom version

Description

This function allows you to make arbitrary changes to how mizer works by allowing you to replace any mizer function with your own version. You should do this only as a last resort, when you find that you can not use the standard mizer extension mechanism to achieve your goal.

Usage

customFunction(name, fun)

Arguments

Details

If the function you need to overwrite is one of the mizer rate functions, then you should use setRateFunction() instead of this function. Similarly you should use ⁠resource_dynamics()<-⁠ to change the resource dynamics and setReproduction() to change the density-dependence in reproduction. You should also investigate whether you can achieve your goal by introducing additional ecosystem components with setComponent().

If you find that your goal really does require you to overwrite a mizer function, please also create an issue on the mizer issue tracker at https://github.com/sizespectrum/mizer/issues to describe your goal, because it will be interesting to the mizer community and may motivate future improvements to the mizer functionality.

Note that customFunction() only overwrites the function used by the mizer code. It does not overwrite the function that is exported by mizer. This will become clear when you run the code in the Examples section.

This function does not in any way check that your replacement function is compatible with mizer. Calling this function can totally break mizer. However you can always undo the effect by reloading mizer with

detach(package:mizer, unload = TRUE)
library(mizer)

Value

No return value, called for side effects

See Also

"Extending mizer": guide to extending mizer

Examples

## Not run: 
fake_project <- function(...) "Fake"
customFunction("project", fake_project)
mizer::project(NS_params) # This will print "Fake"
project(NS_params) # This will still use the old project() function
# To undo the effect:
customFunction("project", project)
mizer::project(NS_params) # This will again use the old project()

## End(Not run)

The minimum weights given by a cutoff species parameter

Description

Internal helper for getBiomass() and for the calibration and matching functions. Returns the ⁠<to>_cutoff⁠ column of the species parameters, with any NAs replaced by the smallest weight in the model. If the model has no such column at all, the smallest weight is used for every species, so that the whole size range is counted.

Usage

cutoff_min_w(params, to = c("biomass", "number"))

Arguments

Value

A numeric vector with one minimum weight for each species.

The default level of information that mizer gives

Description

Returns the mizer_info_level option if it is set and fallback otherwise. This is the default of the info_level argument of the functions that report information, so that options(mizer_info_level = 0) quietens mizer as a whole, including the functions that have no info_level argument of their own, such as species_params<-() and the rate setters.

Usage

default_info_level(fallback = 3)

Arguments

Details

Extension packages should use this as the default of their own info_level argument, so that a constructor or setter of theirs follows the option like mizer's own do:

newFooParams <- function(species_params, ...,
                         info_level = default_info_level()) {
    newMultispeciesParams(species_params, info_level = info_level, ...)
}

Take the argument explicitly like that rather than hard-coding a value in the call, which would make a user's own info_level collide with it.

Value

A single number, or NA to leave the reporting to a handler further out.

Examples

default_info_level()

# Setting the option changes what every reporting function defaults to.
old <- options(mizer_info_level = 1)
default_info_level()
options(old)

Set defaults for predation kernel parameters

Description

If the predation kernel type has not been specified for a species, then it is set to "lognormal" and the default values are set for the parameters beta and sigma.

Usage

default_pred_kernel_params(object)

Arguments

Value

The object with updated columns in the species params data frame.

Default editions

Description

Function to set and get which edition of default choices is being used.

Usage

defaults_edition(edition = NULL)

Arguments

Details

The mizer functions for creating new models make a lot of choices for default values for parameters that are not provided by the user. Sometimes we find better ways to choose the defaults and update mizer accordingly. When we do this, we will increase the edition number.

If you call defaults_edition() without an argument it returns the currently active edition. Otherwise it sets the active edition to the given value.

Users who want their existing code for creating models not to change behaviour when run with future versions of mizer should explicitly set the desired defaults edition at the top of their code.

The most recent edition is edition 2. It will become the default in the next release. The current default is edition 1. The following defaults are changed in edition 2:

• catchability = 0.3 instead of 1

• ⁠initial effort⁠ = 1 instead of 0

• gamma is set to ensure a feeding level of f0 for larvae with the current value of interaction_resource instead of interaction_resource = 1'.

• initial_n is set using get_steady_state_n() instead of the rather arbitrary old choice.

• In setReproduction(), psi is no longer forced to 1 above w_repro_max; its value is determined entirely by the maturity ogive and the reproductive proportion.

Value

If edition is NULL, the currently active edition number. If edition is supplied, the function sets the global mizer_defaults_edition option, emits a message, and returns the supplied value invisibly.

Define S4 marker classes for a set of dispatch extensions

Description

Creates a linear inheritance chain of S4 classes: the outermost extension extends the next, which extends the next, down to the base MizerParams / MizerSim class. Existing classes are checked for compatibility instead of being redefined.

Usage

defineExtensionClasses(extensions)

Arguments

Value

Invisibly, the named character vector of dispatch extensions.

Define an S4 class or verify it extends the expected parent

Description

If class does not yet exist, defines it as a virtual-free S4 class that contains parent, registered in .GlobalEnv. If class already exists, stops with an error unless it already extends parent.

Usage

defineOrCheckClass(class, parent)

Arguments

Value

Invisibly, class.

Jacobian converting between two density measures

Description

Jacobian converting between two density measures

Usage

density_measure_jacobian(from, to, w, l, b)

Arguments

Value

A numeric vector by which to multiply a density with respect to from to obtain the density with respect to to.

Factor relating a density measure to a density with respect to weight

Description

Writing N_w for the density with respect to weight, the density with respect to measure m is N_w times the factor returned here. With the allometric weight-length relationship w = a l^b these factors are 1 for "w", w for "log_w", dw/dl = b w / l for "l" and l\,dw/dl = b w for "log_l".

Usage

density_measure_weight(measure, w, l, b)

Arguments

Value

A numeric vector of factors.

Density measures a spectrum can be expressed in

Description

A size spectrum is a density, and a density only has a meaning together with the variable it is a density with respect to. That variable is one of

"w"

a density with respect to weight, e.g. numbers per gram.

"log_w"

a density with respect to logarithmic weight, e.g. numbers per log weight interval.

"l"

a density with respect to length, e.g. numbers per cm.

"log_l"

a density with respect to logarithmic length.

NA

not a density, e.g. a rate or a dimensionless quantity. Such values are left alone when the size axis changes.

Usage

density_measures

Format

A character vector of the four density measures.

Details

Mizer arrays are indexed by the model's weight grid, so an array that holds a density (type = "density", see array_types) always holds one with respect to weight. The other measures arise for quantities the spectrum plots compute on the fly: plotSpectra(per_log_size = TRUE) shows a density with respect to logarithmic weight, and either can be restated per unit length by size_axis = "l".

The size unit appearing in the units of a density

Description

The size unit appearing in the units of a density

Usage

density_size_unit(measure)

Arguments

Value

"g" or "cm", or NA_character_ for a density with respect to a logarithmic size, whose units carry no size unit.

The density measure a plot calls for

Description

A density is expressed with respect to two independent choices: the size variable, which follows size_axis, and whether it is per size or per logarithmic size, which follows per_log_size. Plotting against a length axis therefore turns a density with respect to weight into one with respect to length, and a density with respect to logarithmic weight into one with respect to logarithmic length.

Usage

density_target_measure(density_wrt, size_axis, per_log_size = NULL)

Arguments

Value

The density measure to express the values in, or NA_character_ if the values are not a density.

Check whether two objects are different

Description

Check whether two objects are numerically different, ignoring all attributes.

Usage

different(a, b)

Arguments

Details

We use this helper function in particular to see if a new value for a slot in MizerParams is different from the existing value in order to give the appropriate messages.

Value

TRUE or FALSE

Collect the extent of each labelled dimension

Description

Internal helper for sizeIntegral() that checks that arrays sharing a dimension label agree on its extent.

Usage

dim_extents(arrays, labels)

Arguments

Value

A named integer vector giving the extent of each label.

The dimensions of an array, given its labels

Description

Internal helper for sizeIntegral(). A scalar has no labels and no dimensions, a vector has one.

Usage

dim_from_labels(x, labels)

Arguments

Value

An integer vector of dimensions, of the same length as labels.

Filter an extension vector to those that participate in S3/S4 dispatch

Description

An extension participates in dispatch if its requirement is NA_character_ (in-development), if an S4 class with its name already exists, or if its loaded package registers S3 dispatch methods for its marker class (see providesDispatchMethods()). The last case lets an installed extension package participate without defining its marker class statically; mizer creates the class dynamically in defineExtensionClasses().

Usage

dispatchExtensions(extensions)

Arguments

Value

A named character vector containing only the dispatch extensions, preserving order.

Measure distance between current and previous state in terms of RDI

Description

This function can be used in projectUntilSettled() to decide when sufficient convergence to steady state has been achieved.

Usage

distanceMaxRelRDI(params, current, previous)

Arguments

Value

The largest absolute relative change in rdi: max(abs((current_rdi - previous_rdi) / previous_rdi)). If any entry of previous_rdi is zero, the result can be infinite.

See Also

Other distance functions: distanceSSLogN()

Measure distance between current and previous state in terms of fish abundances

Description

Calculates the sum squared difference between log(N) in current and previous state. This function can be used in projectUntilSettled() to decide when sufficient convergence to steady state has been achieved.

Usage

distanceSSLogN(params, current, previous)

Arguments

Value

The sum of squares of the difference in the logs of the (nonzero) fish abundances n, ignoring entries where either state has zero abundance: sum((log(current$n) - log(previous$n))^2)

See Also

Other distance functions: distanceMaxRelRDI()

Length based double-sigmoid selectivity function

Description

A hump-shaped selectivity function with a sigmoidal rise and an independent sigmoidal drop-off. This drop-off is what distinguishes this from the function sigmoid_length() and it is intended to model the escape of large individuals from the fishing gear.

Usage

double_sigmoid_length(w, l25, l50, l50_right, l25_right, species_params, ...)

Arguments

Details

You would not usually call this function directly. Instead, set the sel_func column in gear_params() to "double_sigmoid_length" and provide the l25, l50, l50_right and l25_right values as additional columns. setFishing() will then call this function automatically when calculating the selectivity array.

The selectivity is obtained as the product of two sigmoidal curves, one rising and one dropping. The sigmoidal rise is based on the two parameters l25 and l50 which determine the length at which 25% and 50% of the stock is selected respectively. The sigmoidal drop-off is based on the two parameters l50_right and l25_right which determine the length at which the selectivity curve has dropped back to 50% and 25% respectively. The selectivity is given by the function

S(l) = \frac{1}{1 + \exp\left(\log(3)\frac{l50 -l}{l50 - l25}\right)}\frac{1}{1 + \exp\left(\log(3)\frac{l50_{right} -l}{l50_{right} - l25_{right}}\right)}

As the size-based model is weight based, and this selectivity function is length based, it uses the length-weight parameters a and b to convert between length and weight.

l = \left(\frac{w}{a}\right)^{1/b}

Value

Vector of selectivities at the given sizes. Requires ⁠l25 < l50 < l50_right < l25_right⁠.

See Also

gear_params() for setting the selectivity parameters.

Other selectivity functions: knife_edge(), knife_edge_length(), sigmoid_length(), sigmoid_weight()

Examples

# Hump-shaped selectivity: rises from l25=10 to l50=15,
# then drops back to 50% at l50_right=40 and 25% at l25_right=50
sp <- list(a = 0.01, b = 3)
w <- c(1, 10, 100, 500, 1000)
double_sigmoid_length(w, l25 = 10, l50 = 15,
                      l50_right = 40, l25_right = 50,
                      species_params = sp)

Create empty MizerParams object of the right size

Description

An internal function. Sets up a valid MizerParams object with all the slots initialised and given dimension names, but with some slots left empty. This function is to be used by other functions to set up full parameter objects.

Usage

emptyParams(
  species_params,
  gear_params = data.frame(),
  no_w = 100,
  min_w = 0.001,
  max_w = NA,
  min_w_pp = 1e-12
)

Arguments

Value

An empty but valid MizerParams object

Size grid

A size grid is created so that the log-sizes are equally spaced. The spacing is chosen so that there will be no_w fish size bins, with the smallest starting at min_w and the largest starting at max_w. For the resource spectrum there is a larger set of bins containing additional bins below min_w, with the same log size. The number of extra bins is such that min_w_pp comes to lie within the smallest bin.

Changes to species params

The species_params slot of the returned MizerParams object may differ from the data frame supplied as argument to this function because default values are set for missing parameters.

See Also

See newMultispeciesParams() for a function that fills the slots left empty by this function.

The predation kernel as used by the encounter quadrature

Description

Returns the kernel array \Phi_i(w_k, w_p) for which

E_i(w_k) = \gamma_i(w_k) \sum_p \Phi_i(w_k, w_p) N^{eff}_i(w_p) w_p \Delta w_p

reproduces exactly the available energy computed by mizerEncounter(), where N^{eff} is the interaction-weighted prey density. It is the kernel that any summary function must use if its result is to be consistent with getEncounter().

Usage

encounter_kernel(params)

Arguments

Details

On the default first-order path this is just the point-sampled kernel returned by pred_kernel(). When second-order bin-averaging is switched on (see second_order_w()) the two differ: setPredKernel() then builds the Fourier-transformed kernel from the kernel integrated over the prey bin, divided by \beta - 1 so that the plain point weight w_p \Delta w_p carried by the prey vector is cancelled. Those bin-integrated weights are recovered here from params@ft_pred_kernel_e by an inverse Fourier transform, which costs one FFT and keeps this helper automatically in step with whatever quadrature setPredKernel() used.

Pair it with the plain point prey weight params@w_full * params@dw_full. That weight is a normalisation which the kernel construction is built to cancel, not a first-order quadrature weight, so it must not be passed through bin_average_weight(): doing so applies the prey-bin integral twice. A summary function that instead pairs the point-sampled pred_kernel() with a bin-averaged prey weight double-counts that quadrature; that was the bug behind issue #474.

Value

An array (predator species x predator size x prey size).

See Also

pred_kernel() for the point-sampled kernel used for plotting and for supplying a custom kernel, second_order_w(), bin_average_weight()

Load (and optionally install) namespaces for all non-NA extensions

Description

For each extension whose requirement is not NA_character_, checks that the package is installed and up-to-date, installs or upgrades via pak::pkg_install() if install = TRUE, then calls loadNamespace().

Usage

ensureExtensionNamespaces(extensions, install = FALSE)

Arguments

Value

Invisibly TRUE.

Expand the size grid

Description

This function expands the size grid in a MizerParams object to the desired min and max size, preserving all existing species. The function is deprecated because you can achieve the same more flexibly with adjustSizeGrid().

Usage

expandSizeGrid(params, ...)

## S3 method for class 'MizerParams'
expandSizeGrid(
  params,
  new_min_w = min(params@w),
  new_max_w = max(params@w),
  preserve_species = params@species_params$species,
  ...
)

Arguments

Value

A new MizerParams object with the updated size grid.

Expand kernel weights indexed by grid offset into a full kernel array

Description

The predation kernel depends only on the predator/prey mass ratio, so on the geometric grid it is a function of the offset m between the predator and prey grid indices alone. phis[i, m + 1] holds the weight of species i at offset m. This helper writes those weights into the (predator species x predator size x prey size) array that the non-FFT code paths work with.

Usage

expand_kernel_offsets(phis, params, species)

Arguments

Value

An array (predator species x predator size x prey size).

Extract the requirement view of an extension chain

Description

The ⁠@extensions⁠ slot may be stored either as a named character vector of requirement strings (the legacy/unversioned form) or as a named list whose entries are length-2 character vectors c(requirement = ..., version = ...). This helper returns the requirement strings as a plain named character vector, which is the form used for dispatch and suffix comparison.

Usage

extensionRequirements(ext)

Arguments

Value

A named character vector of requirement strings.

Get the recorded version stamp for one extension on an object

Description

Get the recorded version stamp for one extension on an object

Usage

extensionVersion(params, name)

Arguments

Value

The recorded version string, or NA_character_ if none.

Extract the version stamps of an extension chain

Description

Returns the version of the extension package that last upgraded the object for each extension, or NA_character_ where no stamp is recorded (including the legacy character-vector form, which carries no versions).

Usage

extensionVersions(ext)

Arguments

Value

A named character vector of version strings (NA where unknown).

Whether any extension recorded on an object needs upgrading

Description

Returns TRUE if, for any extension recorded in the object's ⁠@extensions⁠ slot whose package is installed, the recorded version stamp is missing (NA) or older than the installed version of that package. A missing stamp counts as needing an upgrade so that objects created before extension-version tracking are brought up to date (and stamped) on first use; this is safe because runExtensionUpgrades() only calls an upgrade method if one is registered and such methods are written to be idempotent.

Usage

extension_needs_upgrading(params)

Arguments

Value

TRUE or FALSE.

Filter plotting data to the requested length limits

Description

Filter plotting data to the requested length limits

Usage

filter_plot_length_limits(plot_dat, llim)

Arguments

Value

plot_dat filtered to the length limits, or unchanged if it has no l column.

Size spectra at end of simulation

Description

Size spectra at end of simulation

Usage

finalN(sim)

finalNResource(sim)

idxFinalT(sim)

Arguments

Value

For finalN(): An ArraySpeciesBySize object (species x size) holding the consumer number densities at the end of the simulation

For finalNResource(): A vector holding the resource number densities at the end of the simulation for all size classes

For idxFinalT(): An integer giving the index for extracting the results for the final time step

Examples

str(finalN(NS_sim))

# This could also be obtained using `N()` and `idxFinalT()`
identical(N(NS_sim)[idxFinalT(NS_sim), , ], finalN(NS_sim))
str(finalNResource(NS_sim))
idx <- idxFinalT(NS_sim)
idx
# This coincides with
length(getTimes(NS_sim))
# and corresponds to the final time
getTimes(NS_sim)[idx]
# We can use this index to extract the result at the final time
identical(N(NS_sim)[idx, , ], finalN(NS_sim))
identical(NResource(NS_sim)[idx, ], finalNResource(NS_sim))

Find the steady state of a model

Description

Puts the model onto a steady state of its own dynamics, changing no parameter: the reproduction rate, the resource and the consumer spectra all settle together at whatever the parameters you already have imply.

Usage

findSteadyState(
  params,
  solver = c("project", "newton"),
  effort = params@initial_effort,
  info_level = default_info_level(),
  ...
)

Arguments

Details

This is the counterpart of tuneSteadyState(), which instead holds the reproduction rate and the resource at the values you supplied and adjusts erepro/R_max and cc_pp to make those values steady. Use this function when the parameters are the thing you want to keep — when asking what state a given model settles into, for instance under a changed fishing effort — and tuneSteadyState() while calibrating.

Nothing being held fixed means the search has more ways to end badly. With solver = "project" the run can settle on a limit cycle or drive a species extinct rather than reach a fixed point, and with either solver reproduction can collapse. Check the "convergence" attribute rather than assuming a fixed point was reached; see the section below.

Value

A MizerParams object with the initial state replaced by the steady state found and no parameter changed. It carries a "convergence" attribute describing the solution found; see projectUntilSettled().

Choosing a solver

solver = "project" (the default) runs the dynamics until they settle, via projectUntilSettled(), and takes the final state. It is exactly that function with the trajectory thrown away; call projectUntilSettled() instead if you want to watch the approach.

solver = "newton" 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. This is the natural entry point for a stability analysis with getStability().

The Newton solver treats the resource densities as unknowns alongside the fish and appends the resource steady-state equation to the system, so the resource density and the feeding levels it implies are self-consistent even where consumers are satiated. That equation is written for the default semichemostat resource dynamics, so solver = "newton" stops with an error for any other resource_dynamics; use solver = "project" there.

The Newton iteration also needs the residual F(N) to be continuous. A custom rate function registered with setRateFunction() that jumps as a function of the abundances makes F discontinuous, and where the equilibrium lies on the switching threshold there is no root at all, because neither branch is in equilibrium there. The solver then stalls (nleqslv termination code 3) and returns an iterate pinned to the threshold. See Discontinuous rate functions.

It also respects the active transport scheme: if the experimental second-order scheme is enabled (see second_order_w()) it solves the steady-state equation of that scheme. With the van Leer reconstruction the residual is only Lipschitz, so the iteration converges to a fixed point of the dynamics but not to machine precision. The unlimited "centred" reconstruction admits an undamped odd-even mode at a steady state with no physical diffusion, giving an ill-conditioned steady-state Jacobian for which the solver is not expected to converge.

What you get back may not be a steady state

The stopping criterion is a proxy. It says that two states t_per years apart differ by less than distance_tol on whatever scale the criterion is measured on; it does not say that the state reached is a fixed point. There are four ways the returned object can fail to be one:

• the run converged on its own scale while the biomasses are still visibly drifting (termination = "distance_tolerance");

• the run reached t_max without converging (termination = "time_limit");

• the run settled on a limit cycle (termination = "cycle_detected"), in which case the state stored is one point on that cycle;

• the run stopped because a species was going extinct (termination = "extinction").

So treat the result as a claim to be checked rather than as a guarantee:

attr(params, "convergence")$attractor  # "fixed_point", "limit_cycle" or NA
attr(params, "convergence")$residual   # largest biomass drift, in 1/year
isSteady(params)                       # TRUE if within tolerance
summary(params)                        # includes the biomass-drift verdict
plot(getSteadyResidual(params))        # which species, and at which sizes

attractor is the field that answers the question: it is "fixed_point" only where the measured biomass drift is within residual_tol, so it cannot be satisfied by a distance function that has merely gone quiet. termination says how the run ended and converged whether the solver met its own criterion; neither is a claim about the state. The last line says where the model is not steady, which is the one to reach for when it is not: a model that is off steady state is usually off in one species or one part of the size range, and the plot names it. See getSteadyResidual() for why the verdict is phrased in terms of biomass drift rather than the largest per-capita rate.

The messages this function prints say the same thing — a converged run whose biomasses are still moving reports the drift and adds "Reduce the tolerance on the distance function to converge further." — but they are suppressed by info_level = 0, so in a script the "convergence" attribute is the reliable check.

Finally, a genuine fixed point need not be a stable one. Use getStability() to find out, and solver = "newton" to converge onto a fixed point that the dynamics themselves would run away from.

See Also

tuneSteadyState(), projectUntilSettled(), isSteady(), getSteadyResidual(), getStability()

Examples

params <- findSteadyState(NS_params, solver = "newton")
plotSpectra(params)

Assemble the flux matrix from growth, diffusion and recruitment rates

Description

Internal helper holding the arithmetic shared by getFlux.MizerParams and getFlux.MizerSim. Keeping it separate lets the MizerSim method resolve the rate functions once and reuse them across all saved time steps.

Usage

flux_from_rates(params, n, g, d, rdd, power = 0, flux_limiter = "none")

Arguments

Value

A plain species x size matrix of fluxes (no mizer array class).

Units string for the flux returned by getFlux()

Description

Units string for the flux returned by getFlux()

Usage

flux_units(power)

Arguments

Value

A character string with the units of the flux.

Format an extension chain as a human-readable string

Description

Format an extension chain as a human-readable string

Usage

formatExtensionChain(extensions)

Arguments

Value

A character string such as "mizerExtB -> mizerExtA", or "<empty>" for a zero-length chain.

A gear name that the model is not already using

Description

A gear name that the model is not already using

Usage

free_gear_name(params)

Arguments

Value

A string naming a gear that does not exist in params.

Which parameters feed which frozen array

Description

A lookup table used by signal_frozen_changes() to decide whether a change the user made can take effect. Each entry is named after a slot of MizerParams that can be frozen and gives the quantity as the user knows it, the call that unfreezes it, the parameters that the setter and its default calculations read, and what kind of parameters those are.

Usage

frozen_rate_params()

Details

The list of parameters does not have to be exhaustive, and deliberately is not: it names the parameters that the setter reads directly together with the main inputs of the default calculations for those parameters. A parameter that is missing simply means that the user is not warned, and is left with the message that the setter itself gives, see signal_not_recalculated(). Listing a parameter that in fact has no influence is the worse mistake, because it warns about a change that did take effect.

Value

A named list of lists with entries quantity, reset_call, params and derived_from.

Gaussian-mixture predation kernel

Description

A predation kernel for which the log predator/prey mass ratio follows a mixture of Gaussian distributions.

Usage

gaussian_mixture_pred_kernel(ppmr, kernel_p, kernel_mean, kernel_sd)

Arguments

Details

Writing the predator mass as w, the prey mass as w_p, and x = \ln(w / w_p), the feeding kernel is

\phi_i(w, w_p) = \sum_j a_{ij} \exp\left[-\frac{(x - \mu_{ij})^2}{2\sigma_{ij}^2}\right], \qquad a_{ij} = \frac{p_{ij}/\sigma_{ij}} {\sum_k p_{ik}/\sigma_{ik}}.

for predator/prey mass ratios greater than or equal to one, and zero for smaller ratios.

This is proportional to the Gaussian-mixture probability density with mixing proportions p_{ij}, means \mu_{ij}, and standard deviations \sigma_{ij}. The scaling makes the sum of the component peak heights equal to one. Consequently the kernel is at most one, and a one-component mixture is identical to lognormal_pred_kernel() with beta = exp(kernel_mean) and sigma = kernel_sd.

The three component parameters are vectors of equal length. When this function is selected in a species parameter data frame, they should be held in the list-columns kernel_p, kernel_mean, and kernel_sd. The values in kernel_p must be non-negative with at least one positive value, but they do not need to sum to one because they are normalised by the function.

Value

A vector giving the value of the predation kernel at each of the predator/prey mass ratios in the ppmr argument.

See Also

setPredKernel()

Other predation kernel: box_pred_kernel(), lognormal_pred_kernel(), power_law_pred_kernel(), truncated_lognormal_pred_kernel()

Examples

ppmr <- exp(seq(0, 12, length.out = 200))
phi <- gaussian_mixture_pred_kernel(
    ppmr,
    kernel_p = c(0.3, 0.7),
    kernel_mean = c(4, 8),
    kernel_sd = c(0.8, 1.5)
)
plot(ppmr, phi, type = "l", log = "x")

Every gear name a model uses

Description

Every gear name a model uses

Usage

gear_names(params)

Arguments

Value

A character vector of gear names.

Gear parameters

Description

These functions allow you to get or set the gear parameters stored in a MizerParams object. These are used by setFishing() to set up the selectivity and catchability and thus together with the fishing effort determine the fishing mortality.

Usage

gear_params(object)

gear_params(object) <- value

is.gear_params(x)

Arguments

Details

The gear_params data has one row for each gear-species pair and one column for each parameter that determines how that gear interacts with that species. The columns are:

• species The name of the species

• gear The name of the gear

• catchability A number specifying how strongly this gear selects this species.

• sel_func The name of the function that calculates the selectivity curve.

• One column for each selectivity parameter needed by the selectivity functions.

For the details see setFishing().

There can optionally also be a column yield_observed that allows you to specify for each gear and species the total annual fisheries yield in grams per year. This is used by plotYieldObservedVsModel(), which adds the yields up over the gears to get the observed yield of each species, see get_yield_observed().

The fishing effort, which is also needed to determine the fishing mortality exerted by a gear is not set via the gear_params data frame but is set with initial_effort() or is specified when calling project().

If you change a gear parameter, this will be used to recalculate the selectivity and catchability arrays by calling setFishing(), unless you have previously set these by hand.

⁠gear_params<-⁠ automatically sets the row names to contain the species name and the gear name, separated by a comma and a space. The last example below illustrates how this facilitates changing an individual gear parameter.

Value

Data frame with gear parameters

is.gear_params() returns TRUE if x is a gear_params object, FALSE otherwise.

See Also

Other functions for setting parameters: setExtDiffusion(), setExtEncounter(), setExtMort(), setFishing(), setInteraction(), setMaxIntakeRate(), setMetabolicRate(), setParams(), setPredKernel(), setReproduction(), setSearchVolume(), species_params(), use_predation_diffusion()

Examples

params <- NS_params

# gears set up in example
gear_params(params)

# setting totally different gears
gear_params(params) <- data.frame(
    gear = c("gear1", "gear2", "gear1"),
    species = c("Cod", "Cod", "Haddock"),
    catchability = c(0.5, 2, 1),
    sel_func = c("sigmoid_weight", "knife_edge", "sigmoid_weight"),
    sigmoidal_weight = c(1000, NA, 800),
    sigmoidal_sigma = c(100, NA, 100),
    knife_edge_size = c(NA, 1000, NA)
    )
gear_params(params)

# changing an individual entry
gear_params(params)["Cod, gear1", "catchability"] <- 0.8

Calculate the total biomass of each species within a size range at each time step.

Description

Calculates the total biomass through time within user defined size limits. The default option is to use the size range starting at the size specified by the biomass_cutoff species parameter, if it is set, or else the full size range of each species. You can specify minimum and maximum weight or length range for the species. Lengths take precedence over weights (i.e. if both min_l and min_w are supplied, only min_l will be used).

Usage

getBiomass(object, use_cutoff = FALSE, ...)

Arguments

Details

When no size range arguments are provided, the function checks if the biomass_cutoff column exists in the species parameters. If it does, those values are used as the minimum weight for each species. For species with NA values in biomass_cutoff, the default minimum weight (smallest weight in the model) is used.

Value

If called with a MizerParams object, a named vector with the biomass in grams for each species in the model. If called with a MizerSim object, an ArrayTimeBySpecies object (time x species) containing the biomass in grams at each time step for all species.

See Also

Other summary functions: getDiet(), getGrowthCurves(), getN(), getSSB(), getSteadyResidual(), getTrophicLevel(), getTrophicLevelBySpecies(), getYield(), getYieldGear()

Examples

biomass <- getBiomass(NS_sim)
biomass["1972", "Herring"]
biomass <- getBiomass(NS_sim, min_w = 10, max_w = 1000)
biomass["1972", "Herring"]

# If species_params contains a `biomass_cutoff`` column, it can be used
# as the minimum weight when use_cutoff = TRUE
species_params(NS_sim@params)$biomass_cutoff <- 10
biomass <- getBiomass(NS_sim, use_cutoff = TRUE)  # Uses biomass_cutoff as min_w
biomass["1972", "Herring"]

Calculate the slope of the community abundance

Description

Calculates the slope of the community abundance by performing a linear regression on the logged total numerical abundance at weight and logged weights (natural logs, not log to base 10, are used). You can specify minimum and maximum weight or length range for the species. Lengths take precedence over weights (i.e. if both min_l and min_w are supplied, only min_l will be used). You can also specify the species to be used in the calculation.

Usage

getCommunitySlope(object, species = NULL, biomass = TRUE, ...)

Arguments

Value

A data.frame with columns slope, intercept and the coefficient of determination R^2 (and a time step column when called with a MizerSim object).

See Also

Other functions for calculating indicators: getMeanMaxWeight(), getMeanWeight(), getProportionOfLargeFish()

Examples

# Slope based on biomass, using all species and sizes
slope_biomass <- getCommunitySlope(NS_sim)
slope_biomass[1, ] # in 1976
slope_biomass[idxFinalT(NS_sim), ] # in 2010

# Slope based on numbers, using all species and sizes
slope_numbers <- getCommunitySlope(NS_sim, biomass = FALSE)
slope_numbers[1, ] # in 1976

# Slope based on biomass, using all species and sizes between 10g and 1000g
slope_biomass <- getCommunitySlope(NS_sim, min_w = 10, max_w = 1000)
slope_biomass[1, ] # in 1976

# Slope based on biomass, using only demersal species and
# sizes between 10g and 1000g
dem_species <- c("Dab","Whiting", "Sole", "Gurnard", "Plaice",
                 "Haddock", "Cod", "Saithe")
slope_biomass <- getCommunitySlope(NS_sim, species = dem_species,
                                   min_w = 10, max_w = 1000)
slope_biomass[1, ] # in 1976

getCommunitySlope(NS_params)

Get critical feeding level

Description

The critical feeding level is the feeding level at which the food intake is just high enough to cover the metabolic costs, with nothing left over for growth or reproduction.

Usage

getCriticalFeedingLevel(params)

Arguments

Value

An ArraySpeciesBySize object (species x size) with the critical feeding level

Examples

str(getFeedingLevel(NS_params))

Get diet of predator at size, resolved by prey species

Description

Calculates the rate at which a predator of a particular species and size consumes biomass of each prey species, resource, and other components of the ecosystem. Returns either the rates in grams/year or the proportion of the total consumption rate.

Usage

getDiet(object, proportion = TRUE, ...)

Arguments

Details

The rates D_{ij}(w) at which a predator of species i and size w consumes biomass from prey species j are calculated from the predation kernel \phi_i(w, w_p), the search volume \gamma_i(w), the feeding level f_i(w), the species interaction matrix \theta_{ij} and the prey abundance density N_j(w_p):

D_{ij}(w, w_p) = (1-f_i(w)) \gamma_i(w) \theta_{ij} \int N_j(w_p) \phi_i(w, w_p) w_p dw_p.

The prey index j runs over all species and the resource.

Extra columns are added for the external encounter rate and for any extra ecosystem components in your model for which you have defined an encounter rate function. These encounter rates are multiplied by 1-f_i(w) to give the rate of consumption of biomass from these extra components.

This function performs the same integration as getEncounter() but does not aggregate over prey species, and multiplies by 1-f_i(w) to get the consumed biomass rather than the available biomass. Outside the range of sizes for a predator species the returned rate is zero. Summing the result of getDiet(proportion = FALSE) over prey therefore reproduces getEncounter(params) * (1 - getFeedingLevel(params)), whichever quadrature scheme the model uses (see second_order_w()).

Value

• MizerParams: An array (predator species x predator size x (prey species + resource + other components)). Dimnames are "predator", "w", and "prey".

• MizerSim: A four-dimensional array (time x predator species x predator size x prey) with the diet at each selected saved time step. If drop = TRUE then dimensions of length 1 are removed.

See Also

plotDiet()

Other summary functions: getBiomass(), getGrowthCurves(), getN(), getSSB(), getSteadyResidual(), getTrophicLevel(), getTrophicLevelBySpecies(), getYield(), getYieldGear()

Examples

diet <- getDiet(NS_params)
str(diet)

# For a MizerSim the diet is returned at each saved time step
sim <- project(NS_params, t_max = 20, effort = 0.5)
# Diet at the saved time steps over years 15 - 20
diet <- getDiet(sim, time_range = c(15, 20))
str(diet)

Get diffusion rate from predation

Description

Calculates the diffusion rate D_i(w) (grams^2/year) for each species. This diffusion rate has two components:

1. The diffusion due due to the variability in prey sizes. This is the diffusion term from the jump-growth equation.

2. Any externally specified diffusion, which is added via setExtDiffusion()

Usage

getDiffusion(object, ...)

Arguments

Details

The diffusion due due to the variability in prey sizes is determined by summing over all prey species and the resource spectrum and then integrating over all prey sizes w_p, weighted by predation kernel \phi(w,w_p):

d_i(w) = (1-f_i(w))(\alpha_i(1-\psi_i(w)))^2\gamma_i(w) \int \left( \theta_{ip} N_R(w_p) + \sum_{j} \theta_{ij} N_j(w_p) \right) \phi_i(w,w_p) w_p^2 \, dw_p.

Here N_j(w) is the abundance density of species j and N_R(w) is the abundance density of resource. The overall prefactor \gamma_i(w) determines the predation power of the predator. It could be interpreted as a search volume and is set with the setSearchVolume() function. The predation kernel \phi(w,w_p) is set with the setPredKernel() function. The species interaction matrix \theta_{ij} is set with setInteraction() and the resource interaction vector \theta_{ip} is taken from the interaction_resource column in species_params(). f(w) is the feeding level calculated with getFeedingLevel(). \psi(w) is the proportion of the available energy that is invested in reproduction instead of growth, obtained with psi().

The diffusion integral is normally evaluated efficiently with a fast Fourier transform, which assumes that the predation kernel depends only on the ratio of predator to prey size. If a custom predation kernel that depends on predator and prey size separately has been set with setPredKernel(), the integral is instead evaluated by direct summation over the full predation kernel, as in getEncounter().

Value

• MizerParams: An ArraySpeciesBySize object (predator species x predator size) with the diffusion rates.

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x predator species x predator size) with the diffusion rates at every time step. If drop = TRUE then dimensions of length 1 will be removed.

References

Datta, S., Delius, G. W. and Law, R. (2010). A jump-growth model for predator-prey dynamics: derivation and application to marine ecosystems. Bulletin of Mathematical Biology, 72(6):1361–1382

See Also

Other rate functions: getEGrowth(), getERepro(), getEReproAndGrowth(), getEncounter(), getFMort(), getFMortGear(), getFeedingLevel(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Analyse the stability of mizer's numerical time step

Description

Computes the eigenvalues \mu_i of the linearised one-step-ahead map at the steady state stored in params@initial_n, for a given step size dt. These describe how mizer's numerical scheme, rather than the model, behaves near the steady state: the map does not amplify perturbations when the spectral radius \max_i|\mu_i| is less than 1.

Usage

getDiscreteStability(params, effort = params@initial_effort, h = 1e-04, dt = 1)

Arguments

Details

This is the numerical counterpart of getStability(), which analyses the model itself and involves no time step at all. Use getStability() to ask whether the steady state of the model is stable, and this function to ask what mizer's solver does at a particular dt. The two can disagree, and that disagreement is the point: the implicit transport solve damps oscillations artificially, so a physically unstable steady state can have a spectral radius below 1 at a large dt, and the simulation then sits at a state the model does not actually hold.

The map that is linearised

One step is what project() takes with method = "euler": the rates are evaluated at the state at the start of the step, and the resulting transport problem is solved implicitly,

A(N^t, n_{pp}^t)\,N^{t+1} = S(N^t, n_{pp}^t),

with the same project_n_loop() C++ Thomas solver and the same spatial scheme (second_order_w()) as the regular dynamics.

Because the rates are evaluated at the input state, the step is not fully implicit, and the discrete eigenvalues therefore cannot be converted into continuous-time eigenvalues by any exact algebraic relation. That conversion is what getStability() avoids by differentiating the rates of change themselves.

The resource is advanced by the model's own resource_dynamics function, the one project() calls. Nothing is substituted for it: the map that is differentiated here reproduces a single project(method = "euler") step exactly, which is what makes the spectral radius a statement about mizer's solver rather than about a nearby scheme.

Value

A named list with the following components:

discrete_eigenvalues

Complex vector of the eigenvalues \mu_i of the one-step map, sorted by decreasing modulus.

spectral_radius

\max_i|\mu_i|. Less than 1 means the numerical scheme is stable at this dt.

stable

Logical: TRUE when spectral_radius < 1.

dt

The step size the map was evaluated at.

n_active

Dimension of the Jacobian.

leading_eigenvectors

The eigenvectors of the two largest-modulus eigenvalues, in the same shape as for getStability().

params

The validated params object the analysis was made at.

Requires smooth dynamics

The finite-difference Jacobian is only meaningful if the rates of change are differentiable at N^*. A custom rate function registered with setRateFunction() that jumps as a function of the abundances breaks this in two ways. If the state sits on the switching threshold, some perturbations straddle it and pick up the jump, and the reported eigenvalues then vary wildly with h. If the state is near but not on the threshold, no perturbation crosses it, and the function silently returns the stability of the single branch the state happens to lie on — which can read as stable for a model whose simulations never settle.

Re-running with a different h is the cheapest check: if the answer moves, do not trust it. See Discontinuous rate functions.

See Also

getStability(), findSteadyState()

Get energy rate available for growth

Description

Calculates the energy rate g_i(w) (grams/year) available by species and size for growth after metabolism, movement and reproduction have been accounted for.

Usage

getEGrowth(object, ...)

Arguments

Details

The growth rate is calculated as the difference between the energy available for reproduction and growth (obtainable with getEReproAndGrowth()) and the energy used for reproduction (obtainable with getERepro()), but is set to 0 if the result would be negative.

Value

• MizerParams: An ArraySpeciesBySize object (species x size) with the somatic growth rates (grams/year).

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x species x size) with the growth rates at every time step. If drop = TRUE then dimensions of length 1 will be removed.

Your own growth rate function

By default getEGrowth() calls mizerEGrowth(). However you can replace this with your own alternative growth rate function. If your function is called "myEGrowth" then you register it in a MizerParams object params with

params <- setRateFunction(params, "EGrowth", "myEGrowth")

Your function will then be called instead of mizerEGrowth(), with the same arguments.

See Also

getERepro(), getEReproAndGrowth()

Other rate functions: getDiffusion(), getERepro(), getEReproAndGrowth(), getEncounter(), getFMort(), getFMortGear(), getFeedingLevel(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <- NS_params
# Project with constant fishing effort for all gears for 20 time steps
sim <- project(params, t_max = 20, effort = 0.5)
# Get the energy at a particular time step
growth <- getEGrowth(params, n = N(sim)[15, , ], n_pp = NResource(sim)[15, ], t = 15)
# Growth rate at this time for Sprat of size 2g
growth["Sprat", "2"]

Get energy rate available for reproduction

Description

Calculates the energy rate (grams/year) available for reproduction after growth and metabolism have been accounted for.

Usage

getERepro(object, ...)

Arguments

Value

• MizerParams: An ArraySpeciesBySize object (species x size) holding

  \psi_i(w)\max(0, E_{r.i}(w))

  where E_{r.i}(w) is the rate at which energy becomes available for growth and reproduction, calculated with getEReproAndGrowth(), and \psi_i(w) is the proportion of this energy that is used for reproduction. Negative values of E_{r.i}(w) are clipped to 0 before multiplying by \psi_i(w). This proportion is taken from the params object and is set with setReproduction().

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x species x size) with the energy for reproduction at every time step. If drop = TRUE then dimensions of length 1 will be removed.

Your own reproduction rate function

By default getERepro() calls mizerERepro(). However you can replace this with your own alternative reproduction rate function. If your function is called "myERepro" then you register it in a MizerParams object params with

params <- setRateFunction(params, "ERepro", "myERepro")

Your function will then be called instead of mizerERepro(), with the same arguments.

See Also

Other rate functions: getDiffusion(), getEGrowth(), getEReproAndGrowth(), getEncounter(), getFMort(), getFMortGear(), getFeedingLevel(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <- NS_params
# Project with constant fishing effort for all gears for 20 time steps
sim <- project(params, t_max = 20, effort = 0.5)
# Get the rate at a particular time step
erepro <- getERepro(params, n = N(sim)[15, , ], n_pp = NResource(sim)[15, ], t = 15)
# Rate at this time for Sprat of size 2g
erepro["Sprat", "2"]

Get energy rate available for reproduction and growth

Description

Calculates the energy rate E_{r.i}(w) (grams/year) available for reproduction and growth after metabolism and movement have been accounted for.

Usage

getEReproAndGrowth(object, ...)

Arguments

Value

• MizerParams: An ArraySpeciesBySize object (species x size) with the energy rate E_{r.i}(w) available for growth and reproduction (grams/year).

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x species x size) with the energy rate at every time step. If drop = TRUE then dimensions of length 1 will be removed.

Your own energy rate function

By default getEReproAndGrowth() calls mizerEReproAndGrowth(). However you can replace this with your own alternative energy rate function. If your function is called "myEReproAndGrowth" then you register it in a MizerParams object params with

params <- setRateFunction(params, "EReproAndGrowth", "myEReproAndGrowth")

Your function will then be called instead of mizerEReproAndGrowth(), with the same arguments.

See Also

The part of this energy rate that is invested into growth is calculated with getEGrowth() and the part that is invested into reproduction is calculated with getERepro().

Other rate functions: getDiffusion(), getEGrowth(), getERepro(), getEncounter(), getFMort(), getFMortGear(), getFeedingLevel(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <- NS_params
# Project with constant fishing effort for all gears for 20 time steps
sim <- project(params, t_max = 20, effort = 0.5)
# Get the energy at a particular time step
e <- getEReproAndGrowth(params, n = N(sim)[15, , ],
                        n_pp = NResource(sim)[15, ], t = 15)
# Rate at this time for Sprat of size 2g
e["Sprat", "2"]

Alias for getERepro()

Description

An alias provided for backward compatibility with mizer version <= 1.0

Usage

getESpawning(object, ...)

Arguments

Value

• MizerParams: An ArraySpeciesBySize object (species x size) holding

  \psi_i(w)\max(0, E_{r.i}(w))

  where E_{r.i}(w) is the rate at which energy becomes available for growth and reproduction, calculated with getEReproAndGrowth(), and \psi_i(w) is the proportion of this energy that is used for reproduction. Negative values of E_{r.i}(w) are clipped to 0 before multiplying by \psi_i(w). This proportion is taken from the params object and is set with setReproduction().

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x species x size) with the energy for reproduction at every time step. If drop = TRUE then dimensions of length 1 will be removed.

Your own reproduction rate function

By default getERepro() calls mizerERepro(). However you can replace this with your own alternative reproduction rate function. If your function is called "myERepro" then you register it in a MizerParams object params with

params <- setRateFunction(params, "ERepro", "myERepro")

Your function will then be called instead of mizerERepro(), with the same arguments.

See Also

Other rate functions: getDiffusion(), getEGrowth(), getEReproAndGrowth(), getEncounter(), getFMort(), getFMortGear(), getFeedingLevel(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <- NS_params
# Project with constant fishing effort for all gears for 20 time steps
sim <- project(params, t_max = 20, effort = 0.5)
# Get the rate at a particular time step
erepro <- getERepro(params, n = N(sim)[15, , ], n_pp = NResource(sim)[15, ], t = 15)
# Rate at this time for Sprat of size 2g
erepro["Sprat", "2"]

Fishing effort used in simulation

Description

Note that the array returned may not be exactly the same as the effort argument that was passed in to project(). This is because only the saved effort is stored (the frequency of saving is determined by the argument t_save).

Usage

getEffort(sim)

Arguments

Value

An array (time x gear) that contains the fishing effort by time and gear.

Examples

str(getEffort(NS_sim))

Get encounter rate

Description

Returns the rate at which a predator of species i and weight w encounters food (grams/year).

Usage

getEncounter(object, ...)

Arguments

Value

• MizerParams: An ArraySpeciesBySize object (predator species x predator size) with the encounter rates.

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x predator species x predator size) with the encounter rates at every time step. If drop = TRUE then dimensions of length 1 will be removed.

Predation encounter

The encounter rate E_i(w) at which a predator of species i and weight w encounters food has contributions from the encounter of fish prey and of resource. This is determined by summing over all prey species and the resource spectrum and then integrating over all prey sizes w_p, weighted by predation kernel \phi(w,w_p):

E_i(w) = \gamma_i(w) \int \left( \theta_{ip} N_R(w_p) + \sum_{j} \theta_{ij} N_j(w_p) \right) \phi_i(w,w_p) w_p \, dw_p.

Here N_j(w) is the abundance density of species j and N_R(w) is the abundance density of resource. The overall prefactor \gamma_i(w) determines the predation power of the predator. It could be interpreted as a search volume and is set with the setSearchVolume() function. The predation kernel \phi(w,w_p) is set with the setPredKernel() function. The species interaction matrix \theta_{ij} is set with setInteraction() and the resource interaction vector \theta_{ip} is taken from the interaction_resource column in params@species_params.

Details

The encounter rate is multiplied by 1-f_0 to obtain the consumption rate, where f_0 is the feeding level calculated with getFeedingLevel(). This is used by the project() function for performing simulations.

The function returns values also for sizes outside the size-range of the species. These values should not be used, as they are meaningless.

If your model contains additional components that you added with setComponent() and for which you specified an encounter_fun function then the encounters of these components will be included in the returned value.

Extension hook

projectEncounter() is the S3 generic used by extension-aware projections. Extension packages can add methods for their marker classes and call NextMethod() to compose encounter-rate changes. The MizerParams method contains the standard mizer calculation and is also exported as mizerEncounter() for compatibility.

Your own encounter function

By default getEncounter() calls mizerEncounter() on models without extensions. However you can replace this with your own alternative encounter function. If your function is called "myEncounter" then you register it in a MizerParams object params with

params <- setRateFunction(params, "Encounter", "myEncounter")

Your function will then be called instead of mizerEncounter(), with the same arguments.

See Also

Other rate functions: getDiffusion(), getEGrowth(), getERepro(), getEReproAndGrowth(), getFMort(), getFMortGear(), getFeedingLevel(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

encounter <- getEncounter(NS_params)
str(encounter)

Get the total fishing mortality rate from all fishing gears by time, species and size.

Description

Calculates the total fishing mortality (in units 1/year) from all gears by species and size and possibly time. See setFishing() for details of how fishing gears are set up.

Usage

getFMort(object, ...)

Arguments

Details

The total fishing mortality is just the sum of the fishing mortalities imposed by each gear, F_i(w)=\sum_g F_{g,i,w}. The fishing mortality for each gear is obtained as catchability x selectivity x effort.

Value

• MizerParams with vector effort: An ArraySpeciesBySize object (species x size) with the fishing mortality rates.

• MizerParams with time-dimensioned effort or MizerSim: An ArrayTimeBySpeciesBySize object (time x species x size).

The effort argument is only used if a MizerParams object is passed in. The effort argument can be a two dimensional array (time x gear), a vector of length equal to the number of gears (each gear has a different effort that is constant in time), or a single numeric value (each gear has the same effort that is constant in time). The order of gears in the effort argument must be the same as in the MizerParams object.

If the object argument is of class MizerSim then the effort slot of the MizerSim object is used and the effort argument is not used.

Your own fishing mortality function

By default getFMort() calls mizerFMort(). However you can replace this with your own alternative fishing mortality function. If your function is called "myFMort" then you register it in a MizerParams object params with

params <- setRateFunction(params, "FMort", "myFMort")

Your function will then be called instead of mizerFMort(), with the same arguments.

See Also

Other rate functions: getDiffusion(), getEGrowth(), getERepro(), getEReproAndGrowth(), getEncounter(), getFMortGear(), getFeedingLevel(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <- NS_params
# Get the total fishing mortality in the initial state
F <- getFMort(params, effort = 1)
str(F)
# Get the initial total fishing mortality when effort is different
# between the four gears:
F <- getFMort(params, effort = c(0.5,1,1.5,0.75))
# Get the total fishing mortality when effort is different
# between the four gears and changes with time:
effort <- array(NA, dim = c(20,4))
effort[, 1] <- seq(from = 0, to = 1, length = 20)
effort[, 2] <- seq(from = 1, to = 0.5, length = 20)
effort[, 3] <- seq(from = 1, to = 2, length = 20)
effort[, 4] <- seq(from = 2, to = 1, length = 20)
F <- getFMort(params, effort = effort)
str(F)
# Get the total fishing mortality using the effort already held in a
# MizerSim object.
sim <- project(params, t_max = 20, effort = 0.5)
F <- getFMort(sim)
F <- getFMort(sim, time_range = c(10, 20))

Get the fishing mortality by time, gear, species and size

Description

Calculates the fishing mortality rate F_{g,i,w} by gear, species and size and possibly time (in units 1/year).

Usage

getFMortGear(object, ...)

Arguments

Value

An array. If the effort argument has a time dimension, or a MizerSim is passed in, the output array has four dimensions (time x gear x species x size). If the effort argument does not have a time dimension (i.e. it is a vector or a single numeric), the output array has three dimensions (gear x species x size).

Note

Here: fishing mortality = catchability x selectivity x effort.

The effort argument is only used if a MizerParams object is passed in. The effort argument can be a two dimensional array (time x gear), a vector of length equal to the number of gears (each gear has a different effort that is constant in time), or a single numeric value (each gear has the same effort that is constant in time). The order of gears in the effort argument must be the same the same as in the MizerParams object. If the effort argument is not supplied, its value is taken from the ⁠@initial_effort⁠ slot in the params object.

If the object argument is of class MizerSim then the effort slot of the MizerSim object is used and the effort argument is not used.

See Also

Other rate functions: getDiffusion(), getEGrowth(), getERepro(), getEReproAndGrowth(), getEncounter(), getFMort(), getFeedingLevel(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <-NS_params
# Get the fishing mortality in initial state
F <- getFMortGear(params, effort = 1)
str(F)
# Get the initial fishing mortality when effort is different
# between the four gears:
F <- getFMortGear(params, effort = c(0.5, 1, 1.5, 0.75))
# Get the fishing mortality when effort is different
# between the four gears and changes with time:
effort <- array(NA, dim = c(20, 4))
effort[, 1] <- seq(from=0, to = 1, length = 20)
effort[, 2] <- seq(from=1, to = 0.5, length = 20)
effort[, 3] <- seq(from=1, to = 2, length = 20)
effort[, 4] <- seq(from=2, to = 1, length = 20)
F <- getFMortGear(params, effort = effort)
str(F)
# Get the fishing mortality using the effort already held in a MizerSim object.
sim <- project(params, t_max = 20, effort = 0.5)
F <- getFMortGear(sim)
F <- getFMortGear(sim, time_range = c(10, 20))

Get feeding level

Description

Returns the feeding level. By default this function uses mizerFeedingLevel() to calculate the feeding level, but this can be overruled via setRateFunction().

Usage

getFeedingLevel(object, ...)

Arguments

Value

• MizerParams: An ArraySpeciesBySize object (predator species x predator size) with the feeding level.

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x predator species x predator size) with the feeding level at every time step. If drop = TRUE then dimensions of length 1 will be removed.

Feeding level

The feeding level f_i(w) is the proportion of its maximum intake rate at which the predator is actually taking in fish. It is calculated from the encounter rate E_i and the maximum intake rate h_i(w) as

f_i(w) = \frac{E_i(w)}{E_i(w)+h_i(w)}.

The encounter rate E_i is passed as an argument or calculated with getEncounter(). The maximum intake rate h_i(w) is taken from the params object, and is set with setMaxIntakeRate(). As a consequence of the above expression for the feeding level, 1-f_i(w) is the proportion of the food available to it that the predator actually consumes.

Your own feeding level function

By default getFeedingLevel() calls mizerFeedingLevel(). However you can replace this with your own alternative feeding level function. If your function is called "myFeedingLevel" then you register it in a MizerParams object params with

params <- setRateFunction(params, "FeedingLevel", "myFeedingLevel")

Your function will then be called instead of mizerFeedingLevel(), with the same arguments.

See Also

Other rate functions: getDiffusion(), getEGrowth(), getERepro(), getEReproAndGrowth(), getEncounter(), getFMort(), getFMortGear(), getFlux(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <- NS_params
# Get initial feeding level
fl <- getFeedingLevel(params)
# Project with constant fishing effort for all gears for 20 time steps
sim <- project(params, t_max = 20, effort = 0.5)
# Get the feeding level at all saved time steps
fl <- getFeedingLevel(sim)
# Get the feeding level for years 15 - 20
fl <- getFeedingLevel(sim, time_range = c(15, 20))

Get flux into size bins

Description

Calculates the flux J_i(w) (numbers/year) entering each size class from the one below it. This is composed of an advective flux from somatic growth and a diffusive flux from the redistribution of individuals.

Usage

getFlux(object, ..., power = 0)

Arguments

Details

At the recruitment size, the flux is simply the recruitment rate R_{dd,i} (see getRDD()). For sizes below the recruitment size the flux is zero.

The flux at weight w is multiplied by w raised to the power given by the power argument, similar to the power argument of plotSpectra(). The default power = 0 returns the flux of individuals (numbers/year). With power = 1 the result is the flux of biomass (grams/year).

Value

• MizerParams: An ArraySpeciesBySize object (species x size) with the flux entering each size class. The units are numbers/year when power = 0 and g^power/year otherwise.

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x species x size) with the flux at every time step. If drop = TRUE then dimensions of length 1 will be removed.

See Also

getEGrowth(), getRDD()

Other rate functions: getDiffusion(), getEGrowth(), getERepro(), getEReproAndGrowth(), getEncounter(), getFMort(), getFMortGear(), getFeedingLevel(), getFluxGradient(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <- NS_params
# Project with constant fishing effort for all gears for 20 time steps
sim <- project(params, t_max = 20, effort = 0.5)
# Get the flux at a particular time step
flux <- getFlux(params, n = N(sim)[15, , ], n_pp = NResource(sim)[15, ], t = 15)
# Flux for Sprat of size 2g
flux["Sprat", "2"]

Get flux gradient

Description

Calculates the flux divergence (J_{j+1} - J_j)/\Delta w_j that appears as the second term in the discretised size-spectrum transport equation

\frac{\partial N_j}{\partial t} + \frac{J_{j+1} - J_j}{\Delta w_j} = -\mu_j N_j.

The bin-boundary fluxes J_j are obtained from getFlux(), which uses the advective-flux scheme stored in the flux entry of the second_order_w slot of params. The flux leaving the largest size class through the upper boundary (J_{K+1}) is evaluated with the same scheme using the boundary condition N_{K+1} = 0.

Usage

getFluxGradient(object, ...)

Arguments

Value

• MizerParams: An ArraySpeciesBySize object (species x size) giving the flux divergence in each size bin, in units of g^{-1} \, \text{year}^{-1}.

• MizerSim: An ArrayTimeBySpeciesBySize object (time step x species x size) with the flux divergence at every saved time step. If drop = TRUE then dimensions of length 1 will be removed.

See Also

getFlux(), second_order_w()

Other rate functions: getDiffusion(), getEGrowth(), getERepro(), getEReproAndGrowth(), getEncounter(), getFMort(), getFMortGear(), getFeedingLevel(), getFlux(), getMort(), getPredMort(), getPredRate(), getRDD(), getRDI(), getRates(), getResourceMort()

Examples

params <- NS_params
fg <- getFluxGradient(params)
sim <- project(params, t_max = 5)
fg_sim <- getFluxGradient(sim)

Get growth curves giving weight as a function of age

Description

Get growth curves giving weight as a function of age

Usage

getGrowthCurves(object, species = NULL, max_age = 20, percentage = FALSE)

Arguments

Value

An array (species x age) containing the weight in grams.

See Also

Other summary functions: getBiomass(), getDiet(), getN(), getSSB(), getSteadyResidual(), getTrophicLevel(), getTrophicLevelBySpecies(), getYield(), getYieldGear()