diff --git a/.claude/skills/upgrade-mizer-data.md b/.claude/skills/upgrade-mizer-data.md new file mode 100644 index 000000000..e413ec9f4 --- /dev/null +++ b/.claude/skills/upgrade-mizer-data.md @@ -0,0 +1,66 @@ +# Upgrade MizerParams / MizerSim data objects + +Use this skill after any change to the `MizerParams` or `MizerSim` S4 class +definition — new slots, removed slots, changed slot types, or new entries in +named lists such as `@rates_funcs`. + +## Steps + +### 1. Update `upgradeParams()` in `R/upgrade.R` + +Add code that patches old objects to the new state. Place it just before the +final block: + +```r +params@mizer_version <- packageVersion("mizer") +params <- validParams(params, info_level = 0) +``` + +Example — adding a new `rates_funcs` entry: + +```r +# Add Diffusion rate function if missing (added in 2.5.4.9122) +if (is.null(params@rates_funcs[["Diffusion"]])) { + params@rates_funcs[["Diffusion"]] <- "mizerDiffusion" +} +``` + +### 2. Bump the version threshold in `needs_upgrading()` in `R/upgrade.R` + +Change the comparison to the **new** version so that existing objects stored +with the old version are detected as needing an upgrade: + +```r +!.hasSlot(params, "mizer_version") || + params@mizer_version < "2.5.4.9122" # ← new version +``` + +### 3. Bump the version in `DESCRIPTION` + +``` +Version: 2.5.4.9122 +``` + +### 4. Reload, upgrade, and save the data objects + +```r +devtools::load_all() +NS_params <- upgradeParams(NS_params) +NS_sim <- upgradeSim(NS_sim) +save(NS_params, file = "data/NS_params.rda", compress = "xz") +save(NS_sim, file = "data/NS_sim.rda", compress = "xz") +``` + +### 5. Reload again before running tests + +```r +devtools::load_all() +devtools::test() +``` + +## Why this matters + +Tests load `NS_params` and `NS_sim` from the `data/` directory. If those +objects predate the class change they will be missing new slots or list +entries, causing errors that look like rate-function or projection failures +rather than the real cause (stale stored objects). diff --git a/.gitignore b/.gitignore index 29b17b045..53e0cab7f 100644 --- a/.gitignore +++ b/.gitignore @@ -9,4 +9,7 @@ docs/dev .vscode .positai .claude +!.claude/ +!.claude/skills/ +!.claude/skills/** .codex diff --git a/AGENTS.md b/AGENTS.md index f0c382495..8f02fd89e 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -42,7 +42,9 @@ devtools::clean_dll(); devtools::load_all() - Use snapshot tests for complex outputs - Run `devtools::document()` after adding or changing exports - Run `devtools::load_all()` before running tests +- After modifying the `MizerParams` or `MizerSim` class (new/removed slots, changes to `@rates_funcs`, etc.), follow the steps in `.claude/skills/upgrade-mizer-data.md`. Claude Code users can invoke this as `/upgrade-mizer-data`. ## Before Submitting -Update `NEWS.md` when adding features or fixing bugs. +- After adding a new file under `R/`, add it to the `Collate:` field in `DESCRIPTION` (roxygen2 does not manage this automatically in this package). +- Update `NEWS.md` when adding features or fixing bugs. diff --git a/DESCRIPTION b/DESCRIPTION index de4b1a520..d4c14bcf3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,7 +19,7 @@ Authors@R: c(person("Gustav", "Delius", email="gustav.delius@york.ac.uk", comment = c(ORCID = "0000-0002-8478-3430")), person("Richard", "Southwell", email="richard.southwell@york.ac.uk", role=c("ctb", "cph"))) -Version: 2.5.4.9113 +Version: 2.5.4.9123 License: GPL-3 Imports: assertthat, @@ -47,6 +47,7 @@ Suggests: roxygen2, knitr, rmarkdown, + quarto, pkgdown, covr, spelling @@ -71,7 +72,8 @@ Collate: 'setMetadata.R' 'setExtMort.R' 'setExtEncounter.R' - 'setDiffusion.R' + 'diffusion.R' + 'setExtDiffusion.R' 'setReproduction.R' 'setResource.R' 'setFishing.R' @@ -114,7 +116,6 @@ Collate: 'steadySingleSpecies.R' 'defaults_edition.R' 'validSpeciesParams.R' -RoxygenNote: 7.3.3 Roxygen: list(markdown = TRUE) Encoding: UTF-8 LazyData: true @@ -122,5 +123,6 @@ URL: https://sizespectrum.org/mizer/, https://github.com/sizespectrum/mizer BugReports: https://github.com/sizespectrum/mizer/issues Language: en-GB RdMacros: lifecycle -VignetteBuilder: knitr +VignetteBuilder: knitr, quarto Config/testthat/edition: 3 +Config/roxygen2/version: 8.0.0 diff --git a/NAMESPACE b/NAMESPACE index 7351e68db..d644b4367 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,7 +3,7 @@ S3method("[",ArraySpeciesBySize) S3method("[",ArraySpeciesByTime) S3method("catchability<-",MizerParams) -S3method("diffusion<-",MizerParams) +S3method("ext_diffusion<-",MizerParams) S3method("ext_encounter<-",MizerParams) S3method("ext_mort<-",MizerParams) S3method("gear_params<-",MizerParams) @@ -43,12 +43,12 @@ S3method(calibrateNumber,MizerParams) S3method(calibrateYield,MizerParams) S3method(catchability,MizerParams) S3method(compareParams,MizerParams) -S3method(diffusion,MizerParams) S3method(distanceMaxRelRDI,MizerParams) S3method(distanceSSLogN,MizerParams) S3method(dw,MizerParams) S3method(dw_full,MizerParams) S3method(expandSizeGrid,MizerParams) +S3method(ext_diffusion,MizerParams) S3method(ext_encounter,MizerParams) S3method(ext_mort,MizerParams) S3method(finalN,MizerSim) @@ -62,6 +62,7 @@ S3method(getCommunitySlope,MizerParams) S3method(getCommunitySlope,MizerSim) S3method(getCriticalFeedingLevel,MizerParams) S3method(getDiet,MizerParams) +S3method(getDiffusion,MizerParams) S3method(getEGrowth,MizerParams) S3method(getERepro,MizerParams) S3method(getEReproAndGrowth,MizerParams) @@ -181,7 +182,7 @@ S3method(search_vol,MizerParams) S3method(selectivity,MizerParams) S3method(setBevertonHolt,MizerParams) S3method(setColours,MizerParams) -S3method(setDiffusion,MizerParams) +S3method(setExtDiffusion,MizerParams) S3method(setExtEncounter,MizerParams) S3method(setExtMort,MizerParams) S3method(setFishing,MizerParams) @@ -207,7 +208,7 @@ S3method(validSim,MizerSim) S3method(w,MizerParams) S3method(w_full,MizerParams) export("catchability<-") -export("diffusion<-") +export("ext_diffusion<-") export("ext_encounter<-") export("ext_mort<-") export("gear_params<-") @@ -231,6 +232,7 @@ export("resource_rate<-") export("search_vol<-") export("selectivity<-") export("species_params<-") +export("use_predation_diffusion<-") export(ArraySpeciesBySize) export(ArraySpeciesByTime) export(BevertonHoltRDD) @@ -262,7 +264,6 @@ export(constant_other) export(customFunction) export(default_pred_kernel_params) export(defaults_edition) -export(diffusion) export(distanceMaxRelRDI) export(distanceSSLogN) export(double_sigmoid_length) @@ -270,6 +271,7 @@ export(dw) export(dw_full) export(emptyParams) export(expandSizeGrid) +export(ext_diffusion) export(ext_encounter) export(ext_mort) export(finalN) @@ -283,6 +285,7 @@ export(getCommunitySlope) export(getComponent) export(getCriticalFeedingLevel) export(getDiet) +export(getDiffusion) export(getEGrowth) export(getERepro) export(getEReproAndGrowth) @@ -367,6 +370,7 @@ export(matchYields) export(maturity) export(melt) export(metab) +export(mizerDiffusion) export(mizerEGrowth) export(mizerERepro) export(mizerEReproAndGrowth) @@ -437,7 +441,7 @@ export(selectivity) export(setBevertonHolt) export(setColours) export(setComponent) -export(setDiffusion) +export(setExtDiffusion) export(setExtEncounter) export(setExtMort) export(setFishing) @@ -464,6 +468,7 @@ export(species_params) export(steady) export(steadySingleSpecies) export(truncated_lognormal_pred_kernel) +export(use_predation_diffusion) export(validEffortVector) export(validGearParams) export(validGivenSpeciesParams) diff --git a/NEWS.md b/NEWS.md index 5286e8859..2528ef35d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,14 @@ # Development version 2.5.4.9102 +- `MizerParams` gains a `use_predation_diffusion` slot (logical, default `FALSE`). + When `FALSE` (the default), [mizerDiffusion()] omits the predation-induced + diffusion term, preserving the behaviour of previous mizer versions. Set to + `TRUE` to enable the jump-growth diffusion term. + +- `diffusion()` / `diffusion<-()` / `setDiffusion()` have been renamed to + `ext_diffusion()` / `ext_diffusion<-()` / `setExtDiffusion()` to follow the + same naming convention as `setExtMort()` and `setExtEncounter()`. + - The `MizerSim` accessors defined in `R/MizerSim-class.R` (`validSim()`, `N()`, `NResource()`, `finalN()`, `finalNResource()`, `idxFinalT()`, `getTimes()`, `getEffort()`, and `getParams()`) are now registered as S3 @@ -115,8 +124,8 @@ - The numerical scheme now supports diffusion in the McKendrick-von Foerster equation, allowing individual variability in growth to be modelled. A new `diffusion` slot in `MizerParams` holds the diffusion coefficient (species x - size). Use `setDiffusion()` / `diffusion()` / `diffusion<-()` to set and - retrieve it. + size). Use `setExtDiffusion()` / `ext_diffusion()` / `ext_diffusion<-()` to + set and retrieve it. - New `getFlux()` function calculates the flux of individuals entering each size class, combining the advective flux from somatic growth and the diffusive flux. - `getRequiredRDD()` is now exported. It calculates the recruitment rate needed diff --git a/R/MizerParams-class.R b/R/MizerParams-class.R index bd453e6a9..04ba9ccd7 100644 --- a/R/MizerParams-class.R +++ b/R/MizerParams-class.R @@ -65,8 +65,8 @@ validMizerParams <- function(object) { length(dim(params@maturity)), length(dim(params@ft_mask)), length(dim(params@catchability)), - length(dim(params@diffusion))) == 2)) { - msg <- "initial_n, psi, intake_max, search_vol, metab, mu_b, ext_encounter, interaction, maturity, ft_mask, catchability and diffusion must all be two dimensions" + length(dim(params@ext_diffusion))) == 2)) { + msg <- "initial_n, psi, intake_max, search_vol, metab, mu_b, ext_encounter, interaction, maturity, ft_mask, catchability and ext_diffusion must all be two dimensions" errors <- c(errors, msg) } # 3D arrays @@ -89,9 +89,9 @@ validMizerParams <- function(object) { dim(params@interaction)[1], dim(params@interaction)[2], dim(params@ft_mask)[1], - dim(params@diffusion)[1]) == + dim(params@ext_diffusion)[1]) == dim(params@species_params)[1])) { - msg <- "The number of species in the model must be consistent across the species_params, psi, intake_max, search_vol, mu_b, ext_encounter, interaction (dim 1), selectivity, catchability, interaction (dim 2) and diffusion slots" + msg <- "The number of species in the model must be consistent across the species_params, psi, intake_max, search_vol, mu_b, ext_encounter, interaction (dim 1), selectivity, catchability, interaction (dim 2) and ext_diffusion slots" errors <- c(errors, msg) } # Check number of size groups @@ -105,9 +105,9 @@ validMizerParams <- function(object) { dim(params@search_vol)[2], dim(params@metab)[2], dim(params@selectivity)[3], - dim(params@diffusion)[2]) == + dim(params@ext_diffusion)[2]) == no_w)) { - msg <- "The number of size bins in the model must be consistent across the w, initial_n, maturity, psi, mu_b, ext_encounter, intake_max, search_vol, metab, selectivity (dim 3) and diffusion slots" + msg <- "The number of size bins in the model must be consistent across the w, initial_n, maturity, psi, mu_b, ext_encounter, intake_max, search_vol, metab, selectivity (dim 3) and ext_diffusion slots" errors <- c(errors, msg) } # Check number of gears @@ -128,8 +128,8 @@ validMizerParams <- function(object) { names(dimnames(params@ext_encounter))[1], names(dimnames(params@selectivity))[2], names(dimnames(params@catchability))[2], - names(dimnames(params@diffusion))[1]) == "sp")) { - msg <- "Name of first dimension of initial_n, maturity, psi, intake_max, search_vol, metab, mu_b, ext_encounter, diffusion and the second dimension of selectivity and catchability must be 'sp'" + names(dimnames(params@ext_diffusion))[1]) == "sp")) { + msg <- "Name of first dimension of initial_n, maturity, psi, intake_max, search_vol, metab, mu_b, ext_encounter, ext_diffusion and the second dimension of selectivity and catchability must be 'sp'" errors <- c(errors, msg) } #interaction dimension names @@ -148,8 +148,8 @@ validMizerParams <- function(object) { names(dimnames(params@search_vol))[2], names(dimnames(params@metab))[2], names(dimnames(params@selectivity))[3], - names(dimnames(params@diffusion))[2]) == "w")) { - msg <- "Name of second dimension of psi, intake_max, search_vol, metab, diffusion and third dimension of selectivity must be 'w'" + names(dimnames(params@ext_diffusion))[2]) == "w")) { + msg <- "Name of second dimension of psi, intake_max, search_vol, metab, ext_diffusion and third dimension of selectivity must be 'w'" errors <- c(errors, msg) } if (!all(c( @@ -171,9 +171,9 @@ validMizerParams <- function(object) { dimnames(params@catchability)[[2]], dimnames(params@interaction)[[1]], dimnames(params@interaction)[[2]], - dimnames(params@diffusion)[[1]]) == + dimnames(params@ext_diffusion)[[1]]) == params@species_params$species)) { - msg <- "The species names of species_params, psi, intake_max, search_vol, metab, mu_b, ext_encounter, selectivity, catchability, interaction and diffusion must all be the same" + msg <- "The species names of species_params, psi, intake_max, search_vol, metab, mu_b, ext_encounter, selectivity, catchability, interaction and ext_diffusion must all be the same" errors <- c(errors, msg) } # Check dimnames of w @@ -229,6 +229,12 @@ validMizerParams <- function(object) { # } # } + # use_predation_diffusion must be a single logical value + if (!is.logical(params@use_predation_diffusion) || length(params@use_predation_diffusion) != 1) { + msg <- "use_predation_diffusion must be a single logical value (TRUE or FALSE)" + errors <- c(errors, msg) + } + # Should not have legacy r_max column (has been renamed to R_max) if ("r_max" %in% names(params@species_params)) { msg <- "The 'r_max' column in species_params should be called 'R_max'. You can use 'validParams()' to upgrade your params object." @@ -307,9 +313,10 @@ validMizerParams <- function(object) { #' \eqn{\mu_{ext.i}(w)}. Changed with [setExtMort()]. #' @slot ext_encounter An array (species x size) that holds the external encounter rate #' \eqn{E_{ext.i}(w)}. Changed with [setExtEncounter()]. -#' @slot diffusion An array (species x size) that holds the rate at which the +#' @slot 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 [diffusion()]. +#' deterministic growth dynamics. Changed with [ext_diffusion()]. #' @slot 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 @@ -380,6 +387,10 @@ validMizerParams <- function(object) { #' @slot 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. +#' @slot 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. #' #' @seealso [project()] [MizerSim()] #' [emptyParams()] [newMultispeciesParams()] @@ -410,7 +421,7 @@ setClass( ft_pred_kernel_p = "array", mu_b = "array", ext_encounter = "array", - diffusion = "array", + ext_diffusion = "array", rr_pp = "numeric", cc_pp = "numeric", resource_dynamics = "character", @@ -433,7 +444,8 @@ setClass( A = "numeric", linecolour = "character", linetype = "character", - ft_mask = "array" + ft_mask = "array", + use_predation_diffusion = "logical" ), ) @@ -688,7 +700,7 @@ emptyParams <- function(species_params, metab = mat1, mu_b = mat1, ext_encounter = mat1, - diffusion = mat1, + ext_diffusion = mat1, ft_pred_kernel_e = ft_pred_kernel, ft_pred_kernel_p = ft_pred_kernel, pred_kernel = array(), @@ -717,6 +729,7 @@ emptyParams <- function(species_params, Mort = "mizerMort", ERepro = "mizerERepro", EGrowth = "mizerEGrowth", + Diffusion = "mizerDiffusion", ResourceMort = "mizerResourceMort", RDI = "mizerRDI", RDD = "BevertonHoltRDD"), @@ -726,7 +739,8 @@ emptyParams <- function(species_params, A = as.numeric(rep(NA, no_sp)), linecolour = linecolour, linetype = linetype, - ft_mask = ft_mask + ft_mask = ft_mask, + use_predation_diffusion = FALSE ) return(params) @@ -946,8 +960,8 @@ validParams <- function(params, info_level = 3) { if (!all(is.finite(params@ext_encounter))) { stop("ext_encounter must not contain non-finite values") } - if (!all(is.finite(params@diffusion))) { - stop("diffusion must not contain non-finite values") + if (!all(is.finite(params@ext_diffusion))) { + stop("ext_diffusion must not contain non-finite values") } if (!all(is.finite(params@selectivity))) { stop("selectivity must not contain non-finite values") diff --git a/R/diffusion.R b/R/diffusion.R new file mode 100644 index 000000000..3b803e0e3 --- /dev/null +++ b/R/diffusion.R @@ -0,0 +1,145 @@ +#' Get diffusion rate from predation +#' +#' @description +#' Calculates the diffusion rate \eqn{D_i(w)} (grams^2/year) for each species. +#' This is the rate at which the abundance density is diffused along the +#' size axis due to the variability in prey sizes. This is the diffusion +#' term from the jump-growth equation. +#' +#' @param params A MizerParams object +#' @param n A matrix of species abundances (species x size). Defaults to the +#' initial abundances in `params`. +#' @param n_pp A vector of the resource abundance by size. Defaults to the +#' initial resource abundances in `params`. +#' @param n_other A list of abundances for other dynamical components. +#' @param t The time for which to do the calculation. +#' @param ... Unused +#' +#' @return An array of dimensions species x size holding the diffusion rates. +#' @export +#' @family rate functions +#' @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 +getDiffusion <- function(params, n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + t = 0, + ...) { + UseMethod("getDiffusion") +} +#' @export +getDiffusion.MizerParams <- function(params, n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + t = 0, + ...) { + params <- validParams(params) + f <- get(params@rates_funcs$Diffusion) + d <- f(params, n = n, n_pp = n_pp, n_other = n_other, t = t, + feeding_level = getFeedingLevel(params, n = n, n_pp = n_pp, + n_other = n_other, t = t), ...) + ArraySpeciesBySize(d, value_name = "Diffusion rate", + units = "g^2/year", params = params) +} + + +#' Calculate diffusion rate +#' +#' @description +#' Calculates the diffusion rate \eqn{D_i(w)} (grams^2/year) for each species. +#' This diffusion rate has two components: +#' 1. Thethe 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()] +#' +#' You would not usually call this function directly but instead use +#' [getDiffusion()], which then calls this function unless an alternative +#' diffusion rate function has been registered, see [setRateFunction()]. +#' +#' @param params A MizerParams object +#' @param n A matrix of species abundances (species x size). +#' @param n_pp A vector of the resource abundance by size +#' @param n_other A list of abundances for other dynamical components +#' @param t The time for which to do the calculation (Not used by standard +#' mizer rate functions but useful for extensions.) +#' @param feeding_level An array (species x size) with the feeding level. +#' If not provided, it is calculated from the given abundances. +#' @param ... Unused +#' +#' @return A two dimensional array (species x size) holding the diffusion rate. +#' @export +mizerDiffusion <- function(params, n, n_pp, n_other, t, feeding_level, ...) { + + if (missing(feeding_level)) { + feeding_level <- getFeedingLevel(params, n = n, n_pp = n_pp, + n_other = n_other, t = t) + } + + # idx_sp are the indices into w_full that correspond to consumer sizes w + idx_sp <- (length(params@w_full) - length(params@w) + 1):length(params@w_full) + + if (isTRUE(params@use_predation_diffusion)) { + # Calculate the diffusion integral + # I_d(w) = sum_prey theta_i * N_prey(w_p) * w_p^2 * dw_p + # This is the same convolution as in mizerEncounter but with w_p^2 * dw_p + # weighting instead of w_p * dw_p. + prey_sq <- outer(params@species_params$interaction_resource, n_pp) + prey_sq[, idx_sp] <- prey_sq[, idx_sp] + params@interaction %*% n + prey_sq <- sweep(prey_sq, 2, params@w_full^2 * params@dw_full, "*") + + # Convolve with the predation kernel via FFT. + # mvfft() transforms each column, so we transpose to get row-wise FFTs, + # following the same pattern as mizerEncounter. + integral_d <- Re(base::t(mvfft(base::t(params@ft_pred_kernel_e) * + mvfft(base::t(prey_sq)), + inverse = TRUE))) / length(params@w_full) + # Keep only the consumer sizes + integral_d <- integral_d[, idx_sp, drop = FALSE] + # Remove numerical noise + integral_d[integral_d < 0] <- 0 + + # D(w) = (1 - f(w)) * gamma(w) * alpha^2 * I_d(w) + alpha <- params@species_params$alpha + D <- (1 - feeding_level) * params@search_vol * alpha^2 * integral_d + dimnames(D) <- dimnames(params@metab) + } else { + D <- matrix(0, nrow = nrow(feeding_level), ncol = ncol(feeding_level), + dimnames = dimnames(params@metab)) + } + + # Add any externally specified diffusion + D <- D + params@ext_diffusion + + return(D) +} + + +#' Get or set the use_predation_diffusion flag +#' +#' Controls whether predation-induced diffusion is included when calculating +#' rates with [mizerDiffusion()]. When `FALSE` (the default), the +#' predation-driven diffusion term is omitted, preserving the behaviour of +#' previous mizer versions. Set to `TRUE` to enable the diffusion term from +#' the jump-growth equation. +#' +#' @param params A MizerParams object. +#' @return `use_predation_diffusion()`: A single logical value. +#' @export +#' @family functions for setting parameters +use_predation_diffusion <- function(params) { + params@use_predation_diffusion +} + +#' @rdname use_predation_diffusion +#' @param value A single logical value (`TRUE` or `FALSE`). +#' @return `use_predation_diffusion<-`: A MizerParams object with the +#' `use_predation_diffusion` flag updated. +#' @export +`use_predation_diffusion<-` <- function(params, value) { + assert_that(is.flag(value)) + params@use_predation_diffusion <- value + params@time_modified <- lubridate::now() + params +} diff --git a/R/extension.R b/R/extension.R index f0def207a..d0d18dd43 100644 --- a/R/extension.R +++ b/R/extension.R @@ -47,6 +47,9 @@ EGrowth = f(params, n = n, n_pp = n_pp, n_other = n_other, t = t, e = rates$e, e_repro = rates$e_repro), + Diffusion = + f(params, n = n, n_pp = n_pp, n_other = n_other, t = t, + feeding_level = rates$feeding_level), PredRate = f(params, n = n, n_pp = n_pp, n_other = n_other, t = t, feeding_level = rates$feeding_level), @@ -89,7 +92,7 @@ } required <- c("encounter", "feeding_level", "pred_rate", "pred_mort", "f_mort", "mort", "resource_mort", "e", "e_repro", - "e_growth", "rdi", "rdd") + "e_growth", "diffusion", "rdi", "rdd") missing <- setdiff(required, names(result)) if (length(missing) > 0L) { stop("The list returned by '", fun, @@ -158,6 +161,7 @@ #' * `e` from [mizerEReproAndGrowth()] #' * `e_repro` from [mizerERepro()] #' * `e_growth` from [mizerEGrowth()] +#' * `diffusion` from [mizerDiffusion()] #' * `rdi` from [mizerRDI()] #' * `rdd` from [BevertonHoltRDD()] #' diff --git a/R/getRequiredRDD.R b/R/getRequiredRDD.R index 61c0392fc..19d7aeef3 100644 --- a/R/getRequiredRDD.R +++ b/R/getRequiredRDD.R @@ -19,7 +19,8 @@ getRequiredRDD.MizerParams <- function(params) { coefs <- get_transport_coefs(params, n = params@initial_n, g = getEGrowth(params), mu = getMort(params), dt, - recruitment_flux = numeric(no_sp)) + recruitment_flux = numeric(no_sp), + d = params@ext_diffusion) reproduction <- params@species_params$erepro # vector of correct length names(reproduction) <- params@species_params$species diff --git a/R/helpers.R b/R/helpers.R index ae8907be9..ad75cfdd6 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -106,7 +106,8 @@ get_steady_state_n <- function(params, g, mu, N0) { # Use get_transport_coefs to compute the coefficients with dt = 1 # and no recruitment flux (since we handle the boundary manually) coefs <- get_transport_coefs(params, n, g, mu, dt = 1, - recruitment_flux = rep(0, no_sp)) + recruitment_flux = rep(0, no_sp), + d = params@ext_diffusion) a <- coefs$a # For steady state, the diagonal term \tilde{B} is B - 1 diff --git a/R/manipulate_species.R b/R/manipulate_species.R index 59c7e8460..1c586be96 100644 --- a/R/manipulate_species.R +++ b/R/manipulate_species.R @@ -250,7 +250,7 @@ addSpecies.MizerParams <- function(params, species_params, gear_params = data.fr p@sc[old_w] <- params@sc p@mu_b[old_sp, old_w] <- params@mu_b p@ext_encounter[old_sp, old_w] <- params@ext_encounter - p@diffusion[old_sp, old_w] <- params@diffusion + p@ext_diffusion[old_sp, old_w] <- params@ext_diffusion p@intake_max[old_sp, old_w] <- params@intake_max p@search_vol[old_sp, old_w] <- params@search_vol p@metab[old_sp, old_w] <- params@metab @@ -381,7 +381,7 @@ removeSpecies.MizerParams <- function(params, species, ...) { p@ft_mask <- params@ft_mask[keep, , drop = FALSE] p@mu_b <- params@mu_b[keep, , drop = FALSE] p@ext_encounter <- params@ext_encounter[keep, , drop = FALSE] - p@diffusion <- params@diffusion[keep, , drop = FALSE] + p@ext_diffusion <- params@ext_diffusion[keep, , drop = FALSE] p@species_params <- p@species_params[keep, , drop = FALSE] p@given_species_params <- p@given_species_params[keep, , drop = FALSE] p@interaction <- params@interaction[keep, keep, drop = FALSE] @@ -490,7 +490,7 @@ renameSpecies.MizerParams <- function(params, replace, ...) { } dimnames(params@mu_b)$sp <- species dimnames(params@ext_encounter)$sp <- species - dimnames(params@diffusion)$sp <- species + dimnames(params@ext_diffusion)$sp <- species dimnames(params@interaction)$predator <- species dimnames(params@interaction)$prey <- species dimnames(params@selectivity)$sp <- species diff --git a/R/newMultispeciesParams.R b/R/newMultispeciesParams.R index 96c85388c..e87d986f1 100644 --- a/R/newMultispeciesParams.R +++ b/R/newMultispeciesParams.R @@ -70,6 +70,7 @@ #' @inheritSection setMetabolicRate Setting metabolic rate #' @inheritSection setExtMort Setting external mortality rate #' @inheritSection setExtEncounter Setting external encounter rate +#' @inheritSection setExtDiffusion Setting external diffusion rate #' @inheritSection setReproduction Setting reproduction #' @inheritSection setFishing Setting fishing #' @inheritSection setResource Setting resource dynamics @@ -327,6 +328,7 @@ newMultispeciesParams <- function( #' @inheritSection setMetabolicRate Setting metabolic rate #' @inheritSection setExtMort Setting external mortality rate #' @inheritSection setExtEncounter Setting external encounter rate +#' @inheritSection setExtDiffusion Setting external diffusion rate #' @inheritSection setReproduction Setting reproduction #' @inheritSection setFishing Setting fishing #' @export @@ -359,6 +361,7 @@ setParams.MizerParams <- function(object, interaction = NULL, params <- setMetabolicRate(params, ...) params <- setExtMort(params, ...) params <- setExtEncounter(params, ...) + params <- setExtDiffusion(params, ...) # setSearchVolume() should be called only after # setMaxIntakeRate() and setPredKernel() params <- setSearchVolume(params, ...) diff --git a/R/project.R b/R/project.R index 603e7fb48..275dc9862 100644 --- a/R/project.R +++ b/R/project.R @@ -470,9 +470,6 @@ project_simple.MizerParams <- c <- matrix(0, nrow = no_sp, ncol = no_w) S <- matrix(0, nrow = no_sp, ncol = no_w) - # Do we have diffusion? - has_diffusion <- any(params@diffusion > 0) - # Loop over time steps ---- for (i_time in 1:steps) { r <- rates_fns$Rates( @@ -509,7 +506,7 @@ project_simple.MizerParams <- ) # * Update species ---- - if (has_diffusion) { + if (any(r$diffusion > 0)) { n <- project_n(params, r, n, dt, a, b, c, S, idx, w_min_idx_array_ref, no_sp, no_w) } else { diff --git a/R/project_methods.R b/R/project_methods.R index ae489d82a..57a510615 100644 --- a/R/project_methods.R +++ b/R/project_methods.R @@ -68,6 +68,10 @@ mizerRates <- function(params, n, n_pp, n_other, r$e_growth <- rates_fns$EGrowth( params, n = n, n_pp = n_pp, n_other = n_other, e_repro = r$e_repro, e = r$e, t = t, ...) + # Calculate the diffusion rate D_i(w) + r$diffusion <- rates_fns$Diffusion( + params, n = n, n_pp = n_pp, n_other = n_other, + feeding_level = r$feeding_level, t = t, ...) ## Mortality ---- # Calculate the predation rate diff --git a/R/project_n.R b/R/project_n.R index 40f1d79fd..b6c28763f 100644 --- a/R/project_n.R +++ b/R/project_n.R @@ -29,7 +29,7 @@ project_n <- function(params, r, n, dt, a, b, c, S, idx, w_min_idx_array_ref, no_sp, no_w) { coefs <- get_transport_coefs(params, n, r$e_growth, r$mort, dt, - recruitment_flux = r$rdd) + recruitment_flux = r$rdd, d = r$diffusion) a <- coefs$a b <- coefs$b c <- coefs$c diff --git a/R/rate_functions.R b/R/rate_functions.R index 289ffd799..cc6882414 100644 --- a/R/rate_functions.R +++ b/R/rate_functions.R @@ -1048,7 +1048,7 @@ getFlux.MizerParams <- function(params, n = initialN(params), no_w <- length(params@w) g <- getEGrowth(params, n = n, n_pp = n_pp, n_other = n_other, t = t) - d <- params@diffusion + d <- params@ext_diffusion dw <- params@dw flux <- matrix(0, nrow = no_sp, ncol = no_w, diff --git a/R/setDiffusion.R b/R/setDiffusion.R deleted file mode 100644 index 8d92a5f71..000000000 --- a/R/setDiffusion.R +++ /dev/null @@ -1,67 +0,0 @@ -#' Set diffusion rate -#' -#' @section Setting diffusion rate: -#' The diffusion rate governs the redistribution of abundance density over body size -#' beyond deterministic growth dynamics. -#' -#' The `diffusion` argument allows you to specify a diffusion rate that depends on -#' species and body size. -#' -#' @param params MizerParams -#' @param diffusion Optional. An array (species x size) holding the diffusion rate. If -#' not supplied, the diffusion rate is left unchanged. Initially it is set to 0. -#' @param reset Unused. Included for interface consistency with other setter -#' functions. -#' @param ... Unused -#' -#' @return `setDiffusion()`: A MizerParams object with updated diffusion rate. -#' @export -#' @family functions for setting parameters -setDiffusion <- function(params, diffusion = NULL, reset = FALSE, ...) { - UseMethod("setDiffusion") -} -#' @rdname setDiffusion -#' @export -setDiffusion.MizerParams <- function(params, diffusion = NULL, reset = FALSE, ...) { - - if (is.null(diffusion)) { - diffusion <- params@diffusion - } - - assert_that(is.array(diffusion), - identical(dim(diffusion), dim(params@diffusion))) - params@diffusion[] <- diffusion - - # Keep old comment if new comment is NULL - if (!is.null(comment(diffusion))) { - comment(params@diffusion) <- comment(diffusion) - } - - params@time_modified <- lubridate::now() - return(params) -} - -#' @rdname setDiffusion -#' @return `diffusion()`: An `ArraySpeciesBySize` object (species x size) with -#' the diffusion rate. -#' @export -diffusion <- function(params) { - UseMethod("diffusion") -} -#' @export -diffusion.MizerParams <- function(params) { - ArraySpeciesBySize(params@diffusion, - value_name = "Diffusion rate", - params = params) -} - -#' @rdname setDiffusion -#' @param value diffusion -#' @export -`diffusion<-` <- function(params, value) { - UseMethod("diffusion<-") -} -#' @export -`diffusion<-.MizerParams` <- function(params, value) { - setDiffusion(params, diffusion = value) -} diff --git a/R/setExtDiffusion.R b/R/setExtDiffusion.R new file mode 100644 index 000000000..8bd941ed2 --- /dev/null +++ b/R/setExtDiffusion.R @@ -0,0 +1,122 @@ +#' Set external diffusion rate +#' +#' @section Setting external diffusion rate: +#' The external diffusion rate allows you to impose additional diffusion of +#' abundance density over body size beyond the predation-driven diffusion +#' computed by [mizerDiffusion()]. +#' +#' The `ext_diffusion` argument allows you to specify a diffusion rate that +#' depends on species and body size. +#' +#' If the `ext_diffusion` argument is not supplied, then the external diffusion +#' rate is calculated as a power law: +#' \deqn{D_{ext.i}(w) = D_{ext.i}\, w^{n_i+1}.}{ +#' D_{ext.i}(w) = D_{ext.i} * w^(n_i+1).} +#' The coefficient \eqn{D_{ext.i}} is taken from the `D_ext` column of the +#' species parameter data frame, which defaults to 0. The exponent +#' \eqn{n_i + 1} uses the `n` column of the species parameter data frame. +#' +#' If the `ext_diffusion` slot has a comment and `reset = FALSE`, then a +#' recalculation from the species parameters is suppressed and a message is +#' issued if the recalculated values would differ from the stored ones. +#' +#' @param params MizerParams +#' @param ext_diffusion Optional. An array (species x size) holding the +#' external diffusion rate. If not supplied, a default is calculated from the +#' `D_ext` and `n` species parameters as described in the section "Setting +#' external diffusion rate". +#' @param reset `r lifecycle::badge("experimental")` +#' If set to TRUE, then the external diffusion rate will be reset to the +#' value calculated from the species parameters, even if it was previously +#' overwritten with a custom value. If set to FALSE (default) then a +#' recalculation from the species parameters will take place only if no +#' custom value has been set. +#' @param ... Unused +#' +#' @return `setExtDiffusion()`: A MizerParams object with updated external +#' diffusion rate. +#' @export +#' @family functions for setting parameters +setExtDiffusion <- function(params, ext_diffusion = NULL, reset = FALSE, ...) { + UseMethod("setExtDiffusion") +} +#' @rdname setExtDiffusion +#' @export +setExtDiffusion.MizerParams <- function(params, ext_diffusion = NULL, + reset = FALSE, ...) { + assert_that(is.flag(reset)) + + if (reset) { + if (!is.null(ext_diffusion)) { + warning("Because you set `reset = TRUE`, the value you provided ", + "for `ext_diffusion` will be ignored and a value will be ", + "calculated from the species parameters.") + ext_diffusion <- NULL + } + comment(params@ext_diffusion) <- NULL + } + + # If ext_diffusion array is supplied, check it, store it and return + if (!is.null(ext_diffusion)) { + if (is.null(comment(ext_diffusion))) { + if (is.null(comment(params@ext_diffusion))) { + comment(ext_diffusion) <- "set manually" + } else { + comment(ext_diffusion) <- comment(params@ext_diffusion) + } + } + assert_that(is.array(ext_diffusion), + identical(dim(ext_diffusion), dim(params@ext_diffusion))) + params@ext_diffusion[] <- ext_diffusion + comment(params@ext_diffusion) <- comment(ext_diffusion) + + params@time_modified <- lubridate::now() + return(params) + } + + # Else recalculate from species params + params <- set_species_param_default(params, "D_ext", 0) + + ext_diffusion <- sweep(outer(params@species_params[["n"]], + params@w, function(x, y) y^(x + 1)), + 1, params@species_params[["D_ext"]], "*") + + # Prevent overwriting slot if it has been commented + if (!is.null(comment(params@ext_diffusion))) { + if (different(ext_diffusion, params@ext_diffusion)) { + message("The external diffusion rate has been commented and ", + "therefore will not be recalculated from the species ", + "parameters.") + } + return(params) + } + params@ext_diffusion[] <- ext_diffusion + + params@time_modified <- lubridate::now() + return(params) +} + +#' @rdname setExtDiffusion +#' @return `ext_diffusion()`: An `ArraySpeciesBySize` object (species x size) +#' with the external diffusion rate. +#' @export +ext_diffusion <- function(params) { + UseMethod("ext_diffusion") +} +#' @export +ext_diffusion.MizerParams <- function(params) { + ArraySpeciesBySize(params@ext_diffusion, + value_name = "Diffusion rate", + params = params) +} + +#' @rdname setExtDiffusion +#' @param value ext_diffusion +#' @export +`ext_diffusion<-` <- function(params, value) { + UseMethod("ext_diffusion<-") +} +#' @export +`ext_diffusion<-.MizerParams` <- function(params, value) { + setExtDiffusion(params, ext_diffusion = value) +} diff --git a/R/species_params.R b/R/species_params.R index 383ed60c5..a3d993f72 100644 --- a/R/species_params.R +++ b/R/species_params.R @@ -36,6 +36,8 @@ #' [setExtMort()]. #' * `E_ext` and `n` are used to set the external encounter rate, see #' [setExtEncounter()]. +#' * `D_ext` and `n` are used to set the external diffusion rate, see +#' [setExtDiffusion()]. #' * `w_mat`, `w_mat25`, `w_repro_max` and `m` are used to set the allocation to #' reproduction, see [setReproduction()]. #' * `pred_kernel_type` specifies the shape of the predation kernel. The default diff --git a/R/transport.R b/R/transport.R index 6d1da21bf..cfab6b6fa 100644 --- a/R/transport.R +++ b/R/transport.R @@ -12,14 +12,11 @@ #' #' @return A list with the coefficients A, B, C and S. #' @noRd -get_transport_coefs <- function(params, n, g, mu, dt, recruitment_flux) { - +get_transport_coefs <- function(params, n, g, mu, dt, recruitment_flux, d) { + no_sp <- nrow(params@species_params) no_w <- length(params@w) - # Diffusion coefficient D_i(w) - d <- params@diffusion # species x size - # Pre-calculate some common terms # dw_j dw <- params@dw diff --git a/R/upgrade.R b/R/upgrade.R index 2727453c7..835ccb2d3 100644 --- a/R/upgrade.R +++ b/R/upgrade.R @@ -17,7 +17,7 @@ needs_upgrading <- function(object) { stop("The object you supplied is neither a MizerParams nor a MizerSim object.") } !.hasSlot(params, "mizer_version") || - params@mizer_version < "2.5.4.9111" + params@mizer_version < "2.5.4.9123" } #' Upgrade MizerParams object from earlier mizer versions @@ -214,6 +214,9 @@ upgradeParams <- function(params) { pnew@initial_n <- params@initial_n pnew@initial_n_pp <- params@initial_n_pp } + if (.hasSlot(params, "diffusion")) { + pnew@ext_diffusion[] <- slot(params, "diffusion") + } if (.hasSlot(params, "initial_n_other")) { pnew@initial_n_other <- params@initial_n_other } @@ -309,11 +312,15 @@ upgradeParams <- function(params) { params <- pnew } - if (!.hasSlot(params, "diffusion")) { + if (!.hasSlot(params, "ext_diffusion")) { mat1 <- array(0, dim = c(nrow(params@species_params), length(params@w)), dimnames = list(sp = params@species_params$species, w = signif(params@w, 3))) - params@diffusion <- mat1 + if (.hasSlot(params, "diffusion")) { + params@ext_diffusion <- slot(params, "diffusion") + } else { + params@ext_diffusion <- mat1 + } } # Before version 2.4 ---- @@ -352,6 +359,17 @@ upgradeParams <- function(params) { params@species_params$is_background <- is.na(params@A) } + # Add Diffusion rate function if missing (added in 2.5.4.9122) + if (is.null(params@rates_funcs[["Diffusion"]])) { + params@rates_funcs[["Diffusion"]] <- "mizerDiffusion" + } + + # Add use_predation_diffusion slot if missing (added in 2.5.4.9xxx) + # Default to FALSE to preserve behaviour of previous mizer versions. + if (!.hasSlot(params, "use_predation_diffusion")) { + params@use_predation_diffusion <- FALSE + } + params@mizer_version <- packageVersion("mizer") params <- validParams(params, info_level = 0) params@time_modified <- lubridate::now() diff --git a/R/validSpeciesParams.R b/R/validSpeciesParams.R index 6dcad6fc9..0c72b62f1 100644 --- a/R/validSpeciesParams.R +++ b/R/validSpeciesParams.R @@ -52,6 +52,7 @@ #' * `z_ext` is set to `0` #' * `d` is set to `n - 1` #' * `E_ext` is set to `0` +#' * `D_ext` is set to `0` #' #' Note that the species parameters returned by these functions are not #' guaranteed to produce a viable model. More checks of the parameters are @@ -77,6 +78,7 @@ validSpeciesParams <- function(species_params) { sp <- set_species_param_default(sp, "z_ext", 0) sp <- set_species_param_default(sp, "d", sp$n - 1) sp <- set_species_param_default(sp, "E_ext", 0) + sp <- set_species_param_default(sp, "D_ext", 0) sp <- set_species_param_default(sp, "is_background", FALSE) return(sp) } diff --git a/data/NS_params.rda b/data/NS_params.rda index 368592d80..0b576f409 100644 Binary files a/data/NS_params.rda and b/data/NS_params.rda differ diff --git a/data/NS_sim.rda b/data/NS_sim.rda index 67f75d9ad..bd49c786b 100644 Binary files a/data/NS_sim.rda and b/data/NS_sim.rda differ diff --git a/man/MizerParams-class.Rd b/man/MizerParams-class.Rd index 5b45e0ee6..0bf576d6a 100644 --- a/man/MizerParams-class.Rd +++ b/man/MizerParams-class.Rd @@ -79,9 +79,10 @@ for each species at size. Changed with \code{\link[=setMetabolicRate]{setMetabol \item{\code{ext_encounter}}{An array (species x size) that holds the external encounter rate \eqn{E_{ext.i}(w)}. Changed with \code{\link[=setExtEncounter]{setExtEncounter()}}.} -\item{\code{diffusion}}{An array (species x size) that holds the rate at which the +\item{\code{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 \code{\link[=diffusion]{diffusion()}}.} +deterministic growth dynamics. Changed with \code{\link[=ext_diffusion]{ext_diffusion()}}.} \item{\code{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 @@ -179,6 +180,11 @@ Used to give consistent line types in plots.} \item{\code{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.} + +\item{\code{use_predation_diffusion}}{A logical flag controlling whether predation-induced +diffusion is included when calculating rates with \code{\link[=mizerDiffusion]{mizerDiffusion()}}. +Defaults to \code{FALSE} to preserve the behaviour of previous mizer versions. +Set to \code{TRUE} to enable the diffusion term from the jump-growth equation.} }} \seealso{ diff --git a/man/completeSpeciesParams.Rd b/man/completeSpeciesParams.Rd index f34e7f611..93fe1d301 100644 --- a/man/completeSpeciesParams.Rd +++ b/man/completeSpeciesParams.Rd @@ -72,6 +72,7 @@ parameters are missing or NA: \item \code{z_ext} is set to \code{0} \item \code{d} is set to \code{n - 1} \item \code{E_ext} is set to \code{0} +\item \code{D_ext} is set to \code{0} } Note that the species parameters returned by these functions are not diff --git a/man/gear_params.Rd b/man/gear_params.Rd index 9b1676b6a..ddc868ec0 100644 --- a/man/gear_params.Rd +++ b/man/gear_params.Rd @@ -79,7 +79,7 @@ gear_params(params)["Cod, gear1", "catchability"] <- 0.8 \code{\link[=validGearParams]{validGearParams()}} Other functions for setting parameters: -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/getDiffusion.Rd b/man/getDiffusion.Rd new file mode 100644 index 000000000..89d18a046 --- /dev/null +++ b/man/getDiffusion.Rd @@ -0,0 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/diffusion.R +\name{getDiffusion} +\alias{getDiffusion} +\title{Get diffusion rate from predation} +\usage{ +getDiffusion( + params, + n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + t = 0, + ... +) +} +\arguments{ +\item{params}{A MizerParams object} + +\item{n}{A matrix of species abundances (species x size). Defaults to the +initial abundances in \code{params}.} + +\item{n_pp}{A vector of the resource abundance by size. Defaults to the +initial resource abundances in \code{params}.} + +\item{n_other}{A list of abundances for other dynamical components.} + +\item{t}{The time for which to do the calculation.} + +\item{...}{Unused} +} +\value{ +An array of dimensions species x size holding the diffusion rates. +} +\description{ +Calculates the diffusion rate \eqn{D_i(w)} (grams^2/year) for each species. +This is the rate at which the abundance density is diffused along the +size axis due to the variability in prey sizes. This is the diffusion +term from the jump-growth equation. +} +\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 +} +\seealso{ +Other rate functions: +\code{\link{getEGrowth}()}, +\code{\link{getERepro}()}, +\code{\link{getEReproAndGrowth}()}, +\code{\link{getEncounter}()}, +\code{\link{getFMort}()}, +\code{\link{getFMortGear}()}, +\code{\link{getFeedingLevel}()}, +\code{\link{getFlux}()}, +\code{\link{getMort}()}, +\code{\link{getPredMort}()}, +\code{\link{getPredRate}()}, +\code{\link{getRDD}()}, +\code{\link{getRDI}()}, +\code{\link{getRates}()}, +\code{\link{getResourceMort}()} +} +\concept{rate functions} diff --git a/man/getEGrowth.Rd b/man/getEGrowth.Rd index 7fe270841..17b833024 100644 --- a/man/getEGrowth.Rd +++ b/man/getEGrowth.Rd @@ -73,6 +73,7 @@ growth["Sprat", "2"] \code{\link[=getERepro]{getERepro()}}, \code{\link[=getEReproAndGrowth]{getEReproAndGrowth()}} Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, \code{\link{getEncounter}()}, diff --git a/man/getERepro.Rd b/man/getERepro.Rd index 81d746fc4..ff83ea2ac 100644 --- a/man/getERepro.Rd +++ b/man/getERepro.Rd @@ -70,6 +70,7 @@ erepro["Sprat", "2"] } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getEReproAndGrowth}()}, \code{\link{getEncounter}()}, diff --git a/man/getEReproAndGrowth.Rd b/man/getEReproAndGrowth.Rd index 4dc88fdd7..537b6d2e9 100644 --- a/man/getEReproAndGrowth.Rd +++ b/man/getEReproAndGrowth.Rd @@ -70,6 +70,7 @@ calculated with \code{\link[=getEGrowth]{getEGrowth()}} and the part that is inv reproduction is calculated with \code{\link[=getERepro]{getERepro()}}. Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEncounter}()}, diff --git a/man/getESpawning.Rd b/man/getESpawning.Rd index 4aaa4a756..58a7d2348 100644 --- a/man/getESpawning.Rd +++ b/man/getESpawning.Rd @@ -70,6 +70,7 @@ erepro["Sprat", "2"] } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getEReproAndGrowth}()}, \code{\link{getEncounter}()}, diff --git a/man/getEncounter.Rd b/man/getEncounter.Rd index bc5fe719f..e79dee760 100644 --- a/man/getEncounter.Rd +++ b/man/getEncounter.Rd @@ -97,6 +97,7 @@ str(encounter) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getFMort.Rd b/man/getFMort.Rd index 032405446..296a4fb24 100644 --- a/man/getFMort.Rd +++ b/man/getFMort.Rd @@ -93,6 +93,7 @@ F <- getFMort(sim, time_range = c(10, 20)) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getFMortGear.Rd b/man/getFMortGear.Rd index e4b684164..5f129d29e 100644 --- a/man/getFMortGear.Rd +++ b/man/getFMortGear.Rd @@ -70,6 +70,7 @@ F <- getFMortGear(sim, time_range = c(10, 20)) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getFeedingLevel.Rd b/man/getFeedingLevel.Rd index 963e3e2e7..ff2b76b21 100644 --- a/man/getFeedingLevel.Rd +++ b/man/getFeedingLevel.Rd @@ -84,6 +84,7 @@ fl <- getFeedingLevel(sim, time_range = c(15, 20)) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getFlux.Rd b/man/getFlux.Rd index d673b31c8..b86486d0f 100644 --- a/man/getFlux.Rd +++ b/man/getFlux.Rd @@ -58,6 +58,7 @@ flux["Sprat", "2"] \code{\link[=getEGrowth]{getEGrowth()}}, \code{\link[=getRDD]{getRDD()}} Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getM2.Rd b/man/getM2.Rd index 9d2aa2037..539c66091 100644 --- a/man/getM2.Rd +++ b/man/getM2.Rd @@ -70,6 +70,7 @@ M2 <- getPredMort(sim, time_range = c(15, 20)) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getM2Background.Rd b/man/getM2Background.Rd index 95753e903..d5db2fda4 100644 --- a/man/getM2Background.Rd +++ b/man/getM2Background.Rd @@ -61,6 +61,7 @@ getResourceMort(params, n = N(sim)[15, , ], n_pp = NResource(sim)[15, ]) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getMort.Rd b/man/getMort.Rd index 93c2d471f..bf35fa92d 100644 --- a/man/getMort.Rd +++ b/man/getMort.Rd @@ -77,6 +77,7 @@ mort["Sprat", "2"] \code{\link[=getPredMort]{getPredMort()}}, \code{\link[=getFMort]{getFMort()}} Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getPredMort.Rd b/man/getPredMort.Rd index 6f19f6c46..c711a5d60 100644 --- a/man/getPredMort.Rd +++ b/man/getPredMort.Rd @@ -73,6 +73,7 @@ M2 <- getPredMort(sim, time_range = c(15, 20)) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getPredRate.Rd b/man/getPredRate.Rd index 07e5e4ee4..2ca1fd4ea 100644 --- a/man/getPredRate.Rd +++ b/man/getPredRate.Rd @@ -71,6 +71,7 @@ pred_rate <- getPredRate(params, n = N(sim)[15, , ], } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getRDD.Rd b/man/getRDD.Rd index a1210822a..a94193258 100644 --- a/man/getRDD.Rd +++ b/man/getRDD.Rd @@ -59,6 +59,7 @@ getRDD(params, n = N(sim)[15, , ], n_pp = NResource(sim)[15, ], t = 15) \code{\link[=getRDI]{getRDI()}} Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getRDI.Rd b/man/getRDI.Rd index f645b9d0a..f2e833c84 100644 --- a/man/getRDI.Rd +++ b/man/getRDI.Rd @@ -76,6 +76,7 @@ getRDI(params, n = N(sim)[15, , ], n_pp = NResource(sim)[15, ], t = 15) \code{\link[=getRDD]{getRDD()}} Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getRates.Rd b/man/getRates.Rd index 8c393a66e..ecb583ca7 100644 --- a/man/getRates.Rd +++ b/man/getRates.Rd @@ -64,6 +64,7 @@ identical(rates$encounter, getEncounter(NS_params)) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getResourceMort.Rd b/man/getResourceMort.Rd index 23e8d78a8..b311e56fb 100644 --- a/man/getResourceMort.Rd +++ b/man/getResourceMort.Rd @@ -61,6 +61,7 @@ getResourceMort(params, n = N(sim)[15, , ], n_pp = NResource(sim)[15, ]) } \seealso{ Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/getZ.Rd b/man/getZ.Rd index 4fc5c9220..282d80e2d 100644 --- a/man/getZ.Rd +++ b/man/getZ.Rd @@ -76,6 +76,7 @@ mort["Sprat", "2"] \code{\link[=getPredMort]{getPredMort()}}, \code{\link[=getFMort]{getFMort()}} Other rate functions: +\code{\link{getDiffusion}()}, \code{\link{getEGrowth}()}, \code{\link{getERepro}()}, \code{\link{getEReproAndGrowth}()}, diff --git a/man/mizerDiffusion.Rd b/man/mizerDiffusion.Rd new file mode 100644 index 000000000..471b6e6d1 --- /dev/null +++ b/man/mizerDiffusion.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/diffusion.R +\name{mizerDiffusion} +\alias{mizerDiffusion} +\title{Calculate diffusion rate} +\usage{ +mizerDiffusion(params, n, n_pp, n_other, t, feeding_level, ...) +} +\arguments{ +\item{params}{A MizerParams object} + +\item{n}{A matrix of species abundances (species x size).} + +\item{n_pp}{A vector of the resource abundance by size} + +\item{n_other}{A list of abundances for other dynamical components} + +\item{t}{The time for which to do the calculation (Not used by standard +mizer rate functions but useful for extensions.)} + +\item{feeding_level}{An array (species x size) with the feeding level. +If not provided, it is calculated from the given abundances.} + +\item{...}{Unused} +} +\value{ +A two dimensional array (species x size) holding the diffusion rate. +} +\description{ +Calculates the diffusion rate \eqn{D_i(w)} (grams^2/year) for each species. +This diffusion rate has two components: +\enumerate{ +\item Thethe diffusion due due to the variability in prey sizes. This is the +diffusion term from the jump-growth equation. +\item Any externally specified diffusion, which is added via \code{\link[=setExtDiffusion]{setExtDiffusion()}} +} + +You would not usually call this function directly but instead use +\code{\link[=getDiffusion]{getDiffusion()}}, which then calls this function unless an alternative +diffusion rate function has been registered, see \code{\link[=setRateFunction]{setRateFunction()}}. +} diff --git a/man/newMultispeciesParams.Rd b/man/newMultispeciesParams.Rd index a99377315..8e9c8292d 100644 --- a/man/newMultispeciesParams.Rd +++ b/man/newMultispeciesParams.Rd @@ -477,6 +477,28 @@ recalculation from the species parameters is suppressed and a message is issued if the recalculated values would differ from the stored ones. } +\section{Setting external diffusion rate}{ + +The external diffusion rate allows you to impose additional diffusion of +abundance density over body size beyond the predation-driven diffusion +computed by \code{\link[=mizerDiffusion]{mizerDiffusion()}}. + +The \code{ext_diffusion} argument allows you to specify a diffusion rate that +depends on species and body size. + +If the \code{ext_diffusion} argument is not supplied, then the external diffusion +rate is calculated as a power law: +\deqn{D_{ext.i}(w) = D_{ext.i}\, w^{n_i+1}.}{ + D_{ext.i}(w) = D_{ext.i} * w^(n_i+1).} +The coefficient \eqn{D_{ext.i}} is taken from the \code{D_ext} column of the +species parameter data frame, which defaults to 0. The exponent +\eqn{n_i + 1} uses the \code{n} column of the species parameter data frame. + +If the \code{ext_diffusion} slot has a comment and \code{reset = FALSE}, then a +recalculation from the species parameters is suppressed and a message is +issued if the recalculated values would differ from the stored ones. +} + \section{Setting reproduction}{ diff --git a/man/setDiffusion.Rd b/man/setDiffusion.Rd deleted file mode 100644 index 1ad15123f..000000000 --- a/man/setDiffusion.Rd +++ /dev/null @@ -1,65 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/setDiffusion.R -\name{setDiffusion} -\alias{setDiffusion} -\alias{setDiffusion.MizerParams} -\alias{diffusion} -\alias{diffusion<-} -\title{Set diffusion rate} -\usage{ -setDiffusion(params, diffusion = NULL, reset = FALSE, ...) - -\method{setDiffusion}{MizerParams}(params, diffusion = NULL, reset = FALSE, ...) - -diffusion(params) - -diffusion(params) <- value -} -\arguments{ -\item{params}{MizerParams} - -\item{diffusion}{Optional. An array (species x size) holding the diffusion rate. If -not supplied, the diffusion rate is left unchanged. Initially it is set to 0.} - -\item{reset}{Unused. Included for interface consistency with other setter -functions.} - -\item{...}{Unused} - -\item{value}{diffusion} -} -\value{ -\code{setDiffusion()}: A MizerParams object with updated diffusion rate. - -\code{diffusion()}: An \code{ArraySpeciesBySize} object (species x size) with -the diffusion rate. -} -\description{ -Set diffusion rate -} -\section{Setting diffusion rate}{ - -The diffusion rate governs the redistribution of abundance density over body size -beyond deterministic growth dynamics. - -The \code{diffusion} argument allows you to specify a diffusion rate that depends on -species and body size. -} - -\seealso{ -Other functions for setting parameters: -\code{\link{gear_params}()}, -\code{\link{setExtEncounter}()}, -\code{\link{setExtMort}()}, -\code{\link{setFishing}()}, -\code{\link{setInitialValues}()}, -\code{\link{setInteraction}()}, -\code{\link{setMaxIntakeRate}()}, -\code{\link{setMetabolicRate}()}, -\code{\link{setParams}()}, -\code{\link{setPredKernel}()}, -\code{\link{setReproduction}()}, -\code{\link{setSearchVolume}()}, -\code{\link{species_params}()} -} -\concept{functions for setting parameters} diff --git a/man/setExtDiffusion.Rd b/man/setExtDiffusion.Rd new file mode 100644 index 000000000..c52a40bc2 --- /dev/null +++ b/man/setExtDiffusion.Rd @@ -0,0 +1,85 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/setExtDiffusion.R +\name{setExtDiffusion} +\alias{setExtDiffusion} +\alias{setExtDiffusion.MizerParams} +\alias{ext_diffusion} +\alias{ext_diffusion<-} +\title{Set external diffusion rate} +\usage{ +setExtDiffusion(params, ext_diffusion = NULL, reset = FALSE, ...) + +\method{setExtDiffusion}{MizerParams}(params, ext_diffusion = NULL, reset = FALSE, ...) + +ext_diffusion(params) + +ext_diffusion(params) <- value +} +\arguments{ +\item{params}{MizerParams} + +\item{ext_diffusion}{Optional. An array (species x size) holding the +external diffusion rate. If not supplied, a default is calculated from the +\code{D_ext} and \code{n} species parameters as described in the section "Setting +external diffusion rate".} + +\item{reset}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +If set to TRUE, then the external diffusion rate will be reset to the +value calculated from the species parameters, even if it was previously +overwritten with a custom value. If set to FALSE (default) then a +recalculation from the species parameters will take place only if no +custom value has been set.} + +\item{...}{Unused} + +\item{value}{ext_diffusion} +} +\value{ +\code{setExtDiffusion()}: A MizerParams object with updated external +diffusion rate. + +\code{ext_diffusion()}: An \code{ArraySpeciesBySize} object (species x size) +with the external diffusion rate. +} +\description{ +Set external diffusion rate +} +\section{Setting external diffusion rate}{ + +The external diffusion rate allows you to impose additional diffusion of +abundance density over body size beyond the predation-driven diffusion +computed by \code{\link[=mizerDiffusion]{mizerDiffusion()}}. + +The \code{ext_diffusion} argument allows you to specify a diffusion rate that +depends on species and body size. + +If the \code{ext_diffusion} argument is not supplied, then the external diffusion +rate is calculated as a power law: +\deqn{D_{ext.i}(w) = D_{ext.i}\, w^{n_i+1}.}{ + D_{ext.i}(w) = D_{ext.i} * w^(n_i+1).} +The coefficient \eqn{D_{ext.i}} is taken from the \code{D_ext} column of the +species parameter data frame, which defaults to 0. The exponent +\eqn{n_i + 1} uses the \code{n} column of the species parameter data frame. + +If the \code{ext_diffusion} slot has a comment and \code{reset = FALSE}, then a +recalculation from the species parameters is suppressed and a message is +issued if the recalculated values would differ from the stored ones. +} + +\seealso{ +Other functions for setting parameters: +\code{\link{gear_params}()}, +\code{\link{setExtEncounter}()}, +\code{\link{setExtMort}()}, +\code{\link{setFishing}()}, +\code{\link{setInitialValues}()}, +\code{\link{setInteraction}()}, +\code{\link{setMaxIntakeRate}()}, +\code{\link{setMetabolicRate}()}, +\code{\link{setParams}()}, +\code{\link{setPredKernel}()}, +\code{\link{setReproduction}()}, +\code{\link{setSearchVolume}()}, +\code{\link{species_params}()} +} +\concept{functions for setting parameters} diff --git a/man/setExtEncounter.Rd b/man/setExtEncounter.Rd index 24ea573b6..f4c8f3740 100644 --- a/man/setExtEncounter.Rd +++ b/man/setExtEncounter.Rd @@ -83,7 +83,7 @@ ext_encounter(params) <- allo_encounter \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, \code{\link{setInitialValues}()}, diff --git a/man/setExtMort.Rd b/man/setExtMort.Rd index aee3df291..3000a02f7 100644 --- a/man/setExtMort.Rd +++ b/man/setExtMort.Rd @@ -103,7 +103,7 @@ ext_mort(params) <- allo_mort \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setFishing}()}, \code{\link{setInitialValues}()}, diff --git a/man/setFishing.Rd b/man/setFishing.Rd index 7ab9b797b..60d843643 100644 --- a/man/setFishing.Rd +++ b/man/setFishing.Rd @@ -179,7 +179,7 @@ str(getInitialEffort(NS_params)) Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setInitialValues}()}, diff --git a/man/setInitialValues.Rd b/man/setInitialValues.Rd index 1ee0d013d..b766a3ceb 100644 --- a/man/setInitialValues.Rd +++ b/man/setInitialValues.Rd @@ -63,7 +63,7 @@ params <- setInitialValues(params, sim) \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/setInteraction.Rd b/man/setInteraction.Rd index b0efcd1a4..42daadb1c 100644 --- a/man/setInteraction.Rd +++ b/man/setInteraction.Rd @@ -82,7 +82,7 @@ getInteraction(params) \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/setMaxIntakeRate.Rd b/man/setMaxIntakeRate.Rd index a22734435..d30c73db0 100644 --- a/man/setMaxIntakeRate.Rd +++ b/man/setMaxIntakeRate.Rd @@ -69,7 +69,7 @@ issued if the recalculated values would differ from the stored ones. \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/setMetabolicRate.Rd b/man/setMetabolicRate.Rd index c4903d985..89a357a5d 100644 --- a/man/setMetabolicRate.Rd +++ b/man/setMetabolicRate.Rd @@ -75,7 +75,7 @@ recalculated values would differ from the stored ones. \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/setParams.Rd b/man/setParams.Rd index ea3072490..a4417db48 100644 --- a/man/setParams.Rd +++ b/man/setParams.Rd @@ -367,6 +367,28 @@ recalculation from the species parameters is suppressed and a message is issued if the recalculated values would differ from the stored ones. } +\section{Setting external diffusion rate}{ + +The external diffusion rate allows you to impose additional diffusion of +abundance density over body size beyond the predation-driven diffusion +computed by \code{\link[=mizerDiffusion]{mizerDiffusion()}}. + +The \code{ext_diffusion} argument allows you to specify a diffusion rate that +depends on species and body size. + +If the \code{ext_diffusion} argument is not supplied, then the external diffusion +rate is calculated as a power law: +\deqn{D_{ext.i}(w) = D_{ext.i}\, w^{n_i+1}.}{ + D_{ext.i}(w) = D_{ext.i} * w^(n_i+1).} +The coefficient \eqn{D_{ext.i}} is taken from the \code{D_ext} column of the +species parameter data frame, which defaults to 0. The exponent +\eqn{n_i + 1} uses the \code{n} column of the species parameter data frame. + +If the \code{ext_diffusion} slot has a comment and \code{reset = FALSE}, then a +recalculation from the species parameters is suppressed and a message is +issued if the recalculated values would differ from the stored ones. +} + \section{Setting reproduction}{ @@ -562,7 +584,7 @@ an effort that varies through time. \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/setPredKernel.Rd b/man/setPredKernel.Rd index 53a5b6ffb..55b2bf6d0 100644 --- a/man/setPredKernel.Rd +++ b/man/setPredKernel.Rd @@ -124,7 +124,7 @@ params<- setPredKernel(params, pred_kernel = pred_kernel) \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/setRateFunction.Rd b/man/setRateFunction.Rd index d1d3a439d..fa32346bb 100644 --- a/man/setRateFunction.Rd +++ b/man/setRateFunction.Rd @@ -58,6 +58,7 @@ following components: \item \code{e} from \code{\link[=mizerEReproAndGrowth]{mizerEReproAndGrowth()}} \item \code{e_repro} from \code{\link[=mizerERepro]{mizerERepro()}} \item \code{e_growth} from \code{\link[=mizerEGrowth]{mizerEGrowth()}} +\item \code{diffusion} from \code{\link[=mizerDiffusion]{mizerDiffusion()}} \item \code{rdi} from \code{\link[=mizerRDI]{mizerRDI()}} \item \code{rdd} from \code{\link[=BevertonHoltRDD]{BevertonHoltRDD()}} } diff --git a/man/setReproduction.Rd b/man/setReproduction.Rd index cd6929e35..1054f4d34 100644 --- a/man/setReproduction.Rd +++ b/man/setReproduction.Rd @@ -197,7 +197,7 @@ ggplot(dff) + geom_line(aes(x = Size, y = Proportion, colour = Type)) \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/setSearchVolume.Rd b/man/setSearchVolume.Rd index f7bda4839..1d4dd089c 100644 --- a/man/setSearchVolume.Rd +++ b/man/setSearchVolume.Rd @@ -77,7 +77,7 @@ issued if the recalculated values would differ from the stored ones. \seealso{ Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/species_params.Rd b/man/species_params.Rd index ab73edf3f..4dbbb6a46 100644 --- a/man/species_params.Rd +++ b/man/species_params.Rd @@ -80,6 +80,8 @@ see \code{\link[=setMetabolicRate]{setMetabolicRate()}}. \code{\link[=setExtMort]{setExtMort()}}. \item \code{E_ext} and \code{n} are used to set the external encounter rate, see \code{\link[=setExtEncounter]{setExtEncounter()}}. +\item \code{D_ext} and \code{n} are used to set the external diffusion rate, see +\code{\link[=setExtDiffusion]{setExtDiffusion()}}. \item \code{w_mat}, \code{w_mat25}, \code{w_repro_max} and \code{m} are used to set the allocation to reproduction, see \code{\link[=setReproduction]{setReproduction()}}. \item \code{pred_kernel_type} specifies the shape of the predation kernel. The default @@ -185,7 +187,7 @@ MizerParams object, in case your own code makes use of them. Other functions for setting parameters: \code{\link{gear_params}()}, -\code{\link{setDiffusion}()}, +\code{\link{setExtDiffusion}()}, \code{\link{setExtEncounter}()}, \code{\link{setExtMort}()}, \code{\link{setFishing}()}, diff --git a/man/use_predation_diffusion.Rd b/man/use_predation_diffusion.Rd new file mode 100644 index 000000000..83788c2de --- /dev/null +++ b/man/use_predation_diffusion.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/diffusion.R +\name{use_predation_diffusion} +\alias{use_predation_diffusion} +\alias{use_predation_diffusion<-} +\title{Get or set the use_predation_diffusion flag} +\usage{ +use_predation_diffusion(params) + +use_predation_diffusion(params) <- value +} +\arguments{ +\item{params}{A MizerParams object.} + +\item{value}{A single logical value (\code{TRUE} or \code{FALSE}).} +} +\value{ +\code{use_predation_diffusion()}: A single logical value. + +\verb{use_predation_diffusion<-}: A MizerParams object with the +\code{use_predation_diffusion} flag updated. +} +\description{ +Controls whether predation-induced diffusion is included when calculating +rates with \code{\link[=mizerDiffusion]{mizerDiffusion()}}. When \code{FALSE} (the default), the +predation-driven diffusion term is omitted, preserving the behaviour of +previous mizer versions. Set to \code{TRUE} to enable the diffusion term from +the jump-growth equation. +} +\seealso{ +Other functions for setting parameters: +\code{\link[=gear_params]{gear_params()}}, +\code{\link[=setExtDiffusion]{setExtDiffusion()}}, +\code{\link[=setExtEncounter]{setExtEncounter()}}, +\code{\link[=setExtMort]{setExtMort()}}, +\code{\link[=setFishing]{setFishing()}}, +\code{\link[=setInitialValues]{setInitialValues()}}, +\code{\link[=setInteraction]{setInteraction()}}, +\code{\link[=setMaxIntakeRate]{setMaxIntakeRate()}}, +\code{\link[=setMetabolicRate]{setMetabolicRate()}}, +\code{\link[=setParams]{setParams()}}, +\code{\link[=setPredKernel]{setPredKernel()}}, +\code{\link[=setReproduction]{setReproduction()}}, +\code{\link[=setSearchVolume]{setSearchVolume()}}, +\code{\link[=species_params]{species_params()}} +} +\concept{functions for setting parameters} diff --git a/man/validSpeciesParams.Rd b/man/validSpeciesParams.Rd index fe74dca4e..8fd3fde41 100644 --- a/man/validSpeciesParams.Rd +++ b/man/validSpeciesParams.Rd @@ -78,6 +78,7 @@ parameters are missing or NA: \item \code{z_ext} is set to \code{0} \item \code{d} is set to \code{n - 1} \item \code{E_ext} is set to \code{0} +\item \code{D_ext} is set to \code{0} } Note that the species parameters returned by these functions are not diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index fdae48870..f7ab80384 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -73,6 +73,7 @@ articles: - cheatsheet - a_multispecies_model_of_the_north_sea - plotting + - predation_diffusion - cohort_dynamics_and_diffusion - extensions - title: Publications @@ -190,7 +191,7 @@ reference: - setInteraction - setMaxIntakeRate - setMetabolicRate - - setDiffusion + - setExtDiffusion - setExtMort - setExtEncounter - setReproduction diff --git a/tests/testthat/_snaps/project.md b/tests/testthat/_snaps/project.md index e66985dfb..6ddac60a3 100644 --- a/tests/testthat/_snaps/project.md +++ b/tests/testthat/_snaps/project.md @@ -118,3 +118,125 @@ 33400 39900 0.000000e+00 0.000000e+00 +# Simulation gives same numerical results with predation diffusion + + Code + sim_d@n[2, 3, ] + Output + 0.001 0.00119 0.00142 0.0017 0.00203 0.00242 + 1.501688e+14 1.325437e+14 1.168046e+14 1.027527e+14 9.021193e+13 7.902661e+13 + 0.00289 0.00345 0.00411 0.00491 0.00586 0.00699 + 6.905886e+13 6.018633e+13 5.230021e+13 4.530334e+13 3.910866e+13 3.363781e+13 + 0.00834 0.00995 0.0119 0.0142 0.0169 0.0202 + 2.881988e+13 2.459048e+13 2.089080e+13 1.766696e+13 1.486949e+13 1.245282e+13 + 0.0241 0.0288 0.0343 0.0409 0.0489 0.0583 + 1.037501e+13 8.597468e+12 7.084756e+12 5.804444e+12 4.726978e+12 3.825561e+12 + 0.0696 0.083 0.0991 0.118 0.141 0.168 + 3.076036e+12 2.456761e+12 1.948477e+12 1.534158e+12 1.198858e+12 9.295455e+11 + 0.201 0.24 0.286 0.342 0.408 0.486 + 7.149409e+11 5.453458e+11 4.124826e+11 3.093384e+11 2.300186e+11 1.696108e+11 + 0.58 0.693 0.827 0.987 1.18 1.4 + 1.240594e+11 9.005054e+10 6.490737e+10 4.649423e+10 3.312988e+10 2.350924e+10 + 1.68 2 2.39 2.85 3.4 4.06 + 1.663364e+10 1.174984e+10 8.297590e+09 5.865751e+09 4.156309e+09 2.955478e+09 + 4.84 5.78 6.9 8.23 9.82 11.7 + 2.110584e+09 1.513925e+09 1.091356e+09 7.910411e+08 5.767238e+08 4.233583e+08 + 14 16.7 19.9 23.8 28.4 33.8 + 3.175404e+08 2.476159e+08 2.094759e+08 1.921884e+08 1.612193e+08 1.142373e+08 + 40.4 48.2 57.5 68.7 81.9 97.8 + 7.495069e+07 4.889503e+07 3.155271e+07 1.963895e+07 1.155456e+07 6.221359e+06 + 117 139 166 198 237 282 + 2.883095e+06 1.248335e+06 5.118298e+05 1.963936e+05 6.930678e+04 2.217540e+04 + 337 402 480 573 683 816 + 6.362020e+03 1.622872e+03 3.656946e+02 7.241314e+01 1.254511e+01 1.894644e+00 + 973 1160 1390 1650 1970 2360 + 2.487720e-01 2.835263e-02 2.805023e-03 2.413702e-04 1.814194e-05 1.199238e-06 + 2810 3350 4000 4780 5700 6800 + 7.042058e-08 3.717669e-09 1.789030e-10 7.969878e-12 3.341405e-13 1.340437e-14 + 8120 9690 11600 13800 16500 19600 + 5.226928e-16 2.007798e-17 7.671555e-19 2.931387e-20 1.121033e-21 4.274254e-23 + 23400 28000 33400 39900 + 1.611671e-24 5.938407e-26 2.105831e-27 7.050712e-29 + +--- + + Code + sim_d@n_pp[2, ] + Output + 2.12e-13 2.53e-13 3.02e-13 3.61e-13 4.3e-13 5.14e-13 + 1.087309e+38 7.457141e+37 5.114366e+37 3.507609e+37 2.405640e+37 1.649871e+37 + 6.13e-13 7.31e-13 8.73e-13 1.04e-12 1.24e-12 1.48e-12 + 1.131539e+37 7.760486e+36 5.322410e+36 3.650293e+36 2.503497e+36 1.716985e+36 + 1.77e-12 2.11e-12 2.52e-12 3.01e-12 3.59e-12 4.28e-12 + 1.177568e+36 8.076169e+35 5.538917e+35 3.798781e+35 2.605336e+35 1.786829e+35 + 5.11e-12 6.1e-12 7.28e-12 8.69e-12 1.04e-11 1.24e-11 + 1.225470e+35 8.404695e+34 5.764231e+34 3.953309e+34 2.711316e+34 1.859515e+34 + 1.48e-11 1.76e-11 2.1e-11 2.51e-11 2.99e-11 3.57e-11 + 1.275320e+34 8.746584e+33 5.998710e+33 4.114123e+33 2.821608e+33 1.935157e+33 + 4.26e-11 5.09e-11 6.07e-11 7.25e-11 8.65e-11 1.03e-10 + 1.327197e+33 9.102380e+32 6.242727e+32 4.281478e+32 2.936386e+32 2.013875e+32 + 1.23e-10 1.47e-10 1.75e-10 2.09e-10 2.5e-10 2.98e-10 + 1.381185e+32 9.472647e+31 6.496669e+31 4.455640e+31 3.055832e+31 2.095795e+31 + 3.56e-10 4.25e-10 5.07e-10 6.05e-10 7.21e-10 8.61e-10 + 1.437368e+31 9.857968e+30 6.760934e+30 4.636881e+30 3.180132e+30 2.181044e+30 + 1.03e-09 1.23e-09 1.46e-09 1.75e-09 2.08e-09 2.49e-09 + 1.495834e+30 1.025894e+30 7.035930e+29 4.825479e+29 3.309476e+29 2.269749e+29 + 2.97e-09 3.54e-09 4.23e-09 5.04e-09 6.02e-09 7.18e-09 + 1.556670e+29 1.067615e+29 7.322054e+28 5.021701e+28 3.444042e+28 2.362032e+28 + 8.57e-09 1.02e-08 1.22e-08 1.46e-08 1.74e-08 2.07e-08 + 1.619954e+28 1.111013e+28 7.619645e+27 5.225764e+27 3.583968e+27 2.457974e+27 + 2.48e-08 2.95e-08 3.53e-08 4.21e-08 5.02e-08 5.99e-08 + 1.685735e+27 1.156112e+27 7.928823e+26 5.437703e+26 3.729234e+26 2.557530e+26 + 7.15e-08 8.53e-08 1.02e-07 1.22e-07 1.45e-07 1.73e-07 + 1.753953e+26 1.202848e+26 8.248939e+25 5.656901e+25 3.879278e+25 2.660194e+25 + 2.07e-07 2.46e-07 2.94e-07 3.51e-07 4.19e-07 5e-07 + 1.824163e+25 1.250833e+25 8.576651e+24 5.880510e+24 4.031691e+24 2.763947e+24 + 5.96e-07 7.12e-07 8.49e-07 1.01e-06 1.21e-06 1.44e-06 + 1.894684e+24 1.298682e+24 8.900623e+23 6.099333e+23 4.179076e+23 2.862896e+23 + 1.72e-06 2.06e-06 2.45e-06 2.93e-06 3.49e-06 4.17e-06 + 1.960869e+23 1.342766e+23 9.192860e+22 6.292028e+22 4.305379e+22 2.945134e+22 + 4.98e-06 5.94e-06 7.09e-06 8.46e-06 1.01e-05 1.2e-05 + 2.014026e+22 1.376851e+22 9.409523e+21 6.428462e+21 4.390436e+21 2.997608e+21 + 1.44e-05 1.71e-05 2.05e-05 2.44e-05 2.91e-05 3.48e-05 + 2.046053e+21 1.396190e+21 9.525155e+20 6.497027e+20 4.430876e+20 3.021442e+20 + 4.15e-05 4.95e-05 5.91e-05 7.05e-05 8.42e-05 1e-04 + 2.060184e+20 1.404689e+20 9.577464e+19 6.530224e+19 4.452685e+19 3.036272e+19 + 0.00012 0.000143 0.000171 0.000204 0.000243 0.00029 + 2.070589e+19 1.412193e+19 9.632943e+18 6.572197e+18 4.485111e+18 3.061742e+18 + 0.000346 0.000413 0.000493 0.000588 0.000702 0.000838 + 2.090786e+18 1.428211e+18 9.758692e+17 6.668941e+17 4.557383e+17 3.113712e+17 + 0.001 0.00119 0.00142 0.0017 0.00203 0.00242 + 2.126433e+17 1.451242e+17 9.896050e+16 6.741486e+16 4.587534e+16 3.118252e+16 + 0.00289 0.00345 0.00411 0.00491 0.00586 0.00699 + 2.117124e+16 1.435793e+16 9.726646e+15 6.582347e+15 4.450080e+15 3.005725e+15 + 0.00834 0.00995 0.0119 0.0142 0.0169 0.0202 + 2.028387e+15 1.367727e+15 9.215574e+14 6.205082e+14 4.175446e+14 2.808125e+14 + 0.0241 0.0288 0.0343 0.0409 0.0489 0.0583 + 1.887624e+14 1.268316e+14 8.518808e+13 5.719992e+13 3.839755e+13 2.577090e+13 + 0.0696 0.083 0.0991 0.118 0.141 0.168 + 1.729424e+13 1.160511e+13 7.787633e+12 5.226453e+12 3.508291e+12 2.355690e+12 + 0.201 0.24 0.286 0.342 0.408 0.486 + 1.582443e+12 1.063616e+12 7.154081e+11 4.816224e+11 3.245777e+11 2.190123e+11 + 0.58 0.693 0.827 0.987 1.18 1.4 + 1.479910e+11 1.001607e+11 6.790902e+10 4.613081e+10 3.140086e+10 2.141986e+10 + 1.68 2 2.39 2.85 3.4 4.06 + 1.464313e+10 1.003196e+10 6.887105e+09 4.737140e+09 3.263756e+09 2.251621e+09 + 4.84 5.78 6.9 8.23 9.82 11.7 + 1.554755e+09 1.073971e+09 7.416993e+08 5.117838e+08 3.525962e+08 0.000000e+00 + 14 16.7 19.9 23.8 28.4 33.8 + 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 40.4 48.2 57.5 68.7 81.9 97.8 + 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 117 139 166 198 237 282 + 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 337 402 480 573 683 816 + 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 973 1160 1390 1650 1970 2360 + 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 2810 3350 4000 4780 5700 6800 + 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 8120 9690 11600 13800 16500 19600 + 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 23400 28000 33400 39900 + 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + diff --git a/tests/testthat/_snaps/project_methods.md b/tests/testthat/_snaps/project_methods.md index 9b88f77ee..54717ecc6 100644 --- a/tests/testthat/_snaps/project_methods.md +++ b/tests/testthat/_snaps/project_methods.md @@ -472,3 +472,54 @@ "value": [0, 0.14003569, 0, 0.16166917, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.15516856, 0, 0.18171766, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.17074229, 0, 0.20425055, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18592246, 0, 0.22957537, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19931068, 0, 0.25803764, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20878045, 0, 0.29002552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.21147359, 0, 0.32597508, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20420514, 0, 0.36637613, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18447594, 0, 0.4117788, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.15185233, 0, 0.46280092, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.10878655, 0, 0.52013618, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.05993608, 0, 0.58456336, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0102479, 0, 0.65695649, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.73829629, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.82968277, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.93234916, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.04767728, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.17721437, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.32269134, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.48604242, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.6694259, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.87524576, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.10617333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.36516804, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.65549556, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.9807409, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.34481256, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.75193245, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.2066034, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.71354288, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.27756672, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9034, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6.5953844, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7.35704041, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8.19043179, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9.09526805, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10.06767486, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11.09856898, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12.171609, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13.26077642, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14.32779915, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.31987889, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.16851041, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.7904892, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 17.09228898, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.97855166, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.36423489, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.18817278, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13.42424084, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11.08601057, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8.22269342, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.90710338, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.21936681, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.88840108, 0, 0, 5.33227098, 0, 0, 0, 0, 0, 0, 0, 0, 10.41794867, 0, 0, 15.88911147, 0, 0, 0, 0, 0, 0, 0, 0, 20.56993074, 0, 0, 30.46998309, 0, 0, 0, 0, 0, 0, 0, 0, 26.28703851, 0, 0, 50.6072171, 0, 0, 0, 0, 0, 0, 0, 0, 24.79342486, 0, 0, 71.34445202, 0, 0, 0, 0, 2.72897895, 0, 0, 0, 26.75122552, 0, 0.92710799, 98.11889641, 0, 0, 1.45013185, 0, 31.81585224, 2.92502179, 0, 0, 29.75220119, 0, 9.41209694, 129.30297448, 0, 0, 10.14835947, 0, 75.24893094, 32.17985068, 0, 0, 30.8220051, 0, 19.47218051, 159.92317178, 0, 0, 22.03912144, 0, 133.90652196, 72.527281, 0, 0, 29.37826288, 0, 30.9453227, 189.3292094, 0, 4.27131921, 37.27028959, 16.09850483, 208.74236874, 126.93900325, 0, 0, 24.78026773, 0, 44.8451552, 205.21505321, 9.79868942, 13.99116607, 57.71317442, 51.49405374, 303.1305361, 200.78963382, 0, 0, 15.87705746, 0, 58.41914878, 191.99385549, 24.15583441, 25.62525209, 80.93854984, 100.70956403, 412.70482579, 289.92944402, 0, 0, 2.14698647, 0, 71.39439996, 185.56713478, 36.1751081, 38.98360259, 100.61503761, 163.10469882, 537.02072415, 393.92626312, 0, 0, 0, 0, 81.93863831, 195.37605628, 53.6114257, 55.1763245, 109.76253421, 248.32556956, 675.30724861, 512.80373286, 0, 0, 0, 0, 86.53720178, 209.28248901, 73.88605796, 70.07380214, 120.24424356, 328.28423769, 820.71610967, 638.1068284, 0, 0, 0, 0, 85.01445398, 223.1921164, 95.10610295, 83.16040954, 136.79530406, 372.19368039, 975.40727828, 772.76990541, 0, 0, 0, 0, 75.85271, 235.76150483, 112.88157355, 93.3260379, 154.45379903, 390.30289943, 1139.62035544, 914.68757371, 0, 0, 0, 0, 57.73642202, 246.15605099, 125.4533447, 99.05671479, 171.29544645, 430.34488554, 1313.65843626, 1065.78738938, 0, 0, 0, 0, 29.97659847, 253.87354977, 135.74868181, 100.83747354, 188.29339933, 487.85371961, 1505.28396301, 1229.57345808, 0, 0, 0, 0, 0, 257.55805055, 139.0388638, 96.09168458, 203.84178094, 543.11684759, 1712.89537302, 1405.2201359, 0, 0, 0, 0, 0, 256.08971982, 135.04844245, 84.49591331, 218.1496525, 597.20946007, 1940.66661839, 1597.22532912, 0, 0, 0, 0, 0, 247.81643868, 122.2920949, 64.53590282, 230.72898118, 646.1698563, 2191.51805938, 1808.14311058, 0, 0, 0, 0, 0, 230.81400745, 100.60672949, 35.14130001, 241.54269781, 692.98935327, 2471.42447934, 2041.55033844, 0, 0, 0, 0, 0, 202.50079329, 67.36222277, 0, 249.70745159, 735.92318086, 2783.87574136, 2296.87321407, 0, 0, 0, 0, 0, 159.74313136, 19.89923682, 0, 254.21601293, 772.00541879, 3132.55286257, 2552.25512079, 0, 0, 0, 0, 0, 98.69243802, 0, 0, 253.82903127, 799.55095652, 3517.93678532, 2691.66310154, 0, 0, 0, 0, 0, 14.62119233, 0, 0, 247.06564323, 815.68403187, 3917.29631252, 2593.03717017, 0, 0, 0, 0, 0, 0, 0, 0, 231.97594545, 816.51348629, 4183.98492225, 2573.67886296, 0, 0, 0, 0, 0, 0, 0, 0, 206.15267973, 797.18604087, 4013.39408382, 2729.80529622, 0, 0, 0, 0, 0, 0, 0, 0, 166.63310148, 751.89385796, 3793.26873821, 2951.952997, 0, 0, 0, 0, 0, 0, 0, 0, 109.73628415, 672.93914153, 3911.65586668, 3196.35492018, 0, 0, 0, 0, 0, 0, 0, 0, 30.96562576, 551.00679017, 4178.98238601, 3451.10567661, 0, 0, 0, 0, 0, 0, 0, 0, 0, 374.60393446, 4484.92069023, 3710.98064955, 0, 0, 0, 0, 0, 0, 0, 0, 0, 129.5109971, 4798.01192983, 3971.03476089, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5105.50393101, 4225.08908879, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5396.69233327, 4465.17096837, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5658.89001281, 4680.9053887, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5876.06412374, 4859.00417877, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6027.87034069, 4982.55926504, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6088.63269474, 5030.22749581, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6026.14024946, 4975.24789547, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5800.20408305, 4784.22682034, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5360.90275283, 4415.69319245, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4646.47393952, 3818.34814797, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3580.76989806, 2928.95030802, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2070.19173473, 1669.77453057, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.02041075e-10, 0] } +# getDiffusion snapshot with use_predation_diffusion + + { + "type": "double", + "attributes": { + "dim": { + "type": "integer", + "attributes": {}, + "value": [12, 100] + }, + "dimnames": { + "type": "list", + "attributes": { + "names": { + "type": "character", + "attributes": {}, + "value": ["sp", "w"] + } + }, + "value": [ + { + "type": "character", + "attributes": {}, + "value": ["Sprat", "Sandeel", "N.pout", "Herring", "Dab", "Whiting", "Sole", "Gurnard", "Plaice", "Haddock", "Cod", "Saithe"] + }, + { + "type": "character", + "attributes": {}, + "value": ["0.001", "0.00119", "0.00142", "0.0017", "0.00203", "0.00242", "0.00289", "0.00345", "0.00411", "0.00491", "0.00586", "0.00699", "0.00834", "0.00995", "0.0119", "0.0142", "0.0169", "0.0202", "0.0241", "0.0288", "0.0343", "0.0409", "0.0489", "0.0583", "0.0696", "0.083", "0.0991", "0.118", "0.141", "0.168", "0.201", "0.24", "0.286", "0.342", "0.408", "0.486", "0.58", "0.693", "0.827", "0.987", "1.18", "1.4", "1.68", "2", "2.39", "2.85", "3.4", "4.06", "4.84", "5.78", "6.9", "8.23", "9.82", "11.7", "14", "16.7", "19.9", "23.8", "28.4", "33.8", "40.4", "48.2", "57.5", "68.7", "81.9", "97.8", "117", "139", "166", "198", "237", "282", "337", "402", "480", "573", "683", "816", "973", "1160", "1390", "1650", "1970", "2360", "2810", "3350", "4000", "4780", "5700", "6800", "8120", "9690", "11600", "13800", "16500", "19600", "23400", "28000", "33400", "39900"] + } + ] + }, + "class": { + "type": "character", + "attributes": {}, + "value": ["ArraySpeciesBySize", "matrix", "array"] + }, + "value_name": { + "type": "character", + "attributes": {}, + "value": ["Diffusion rate"] + }, + "units": { + "type": "character", + "attributes": {}, + "value": ["g^2/year"] + } + }, + "value": [0, 3671.22474205, 2.45346456e-07, 3664.73433203, 6.3017814e-06, 1.81690835e-06, 0.00005985, 2.34988764e-07, 2.32894613e-07, 0.11730095, 3.92111199e-06, 3.82926603e-06, 0, 4087.09299053, 5.45960151e-07, 4128.81082917, 4.70420119e-06, 3.64640934e-06, 0.00003305, 2.40096661e-07, 5.72795095e-07, 0.0744323, 6.22107344e-06, 4.64998237e-06, 8.00249722e-08, 4523.67632403, 0, 4651.57458712, 2.94695605e-06, 2.49029514e-06, 0.00001845, 0, 4.78042441e-07, 0.04689794, 5.08549131e-06, 4.001913e-06, 6.73376057e-08, 4963.78495517, 0, 5240.43629762, 3.55152463e-06, 3.8625661e-06, 0.00001063, 0, 8.6259499e-07, 0.02934275, 7.55069456e-06, 4.27436742e-06, 8.70590159e-08, 5377.76200374, 0, 5903.7404642, 2.93392882e-06, 3.91094987e-06, 6.76535926e-06, 0, 5.96754674e-07, 0.01823021, 8.53674018e-06, 5.17149151e-06, 2.52613105e-07, 5719.77688089, 0, 6650.8825712, 3.6028721e-06, 6.74645733e-06, 4.38500211e-06, 1.49367613e-08, 8.08901445e-07, 0.01124806, 0.00001409, 4.64437033e-06, 3.74629995e-07, 5927.68262824, 0, 7492.44079181, 3.5079453e-06, 4.50880253e-06, 4.60040445e-06, 4.03568541e-07, 8.21614685e-07, 0.00689174, 0.00001084, 2.28145627e-06, 7.54956748e-07, 5932.00395857, 0, 8440.32397966, 4.61149444e-06, 7.98066575e-06, 3.55172207e-06, 1.82030727e-07, 2.64926416e-07, 0.00419565, 0.00001561, 4.57695344e-06, 9.64121202e-07, 5678.58706095, 0, 9507.93787306, 3.20349963e-06, 5.34332037e-06, 4.02037509e-06, 0, 0, 0.00253627, 9.18925363e-06, 1.75136161e-06, 1.47834931e-06, 5159.62725582, 0, 10710.37163654, 4.89762018e-06, 8.63757276e-06, 4.00475762e-06, 0, 0, 0.0015247, 0.00001806, 1.03256671e-06, 1.75208049e-06, 4432.36752084, 0, 12064.60706413, 2.04208074e-06, 5.79024534e-06, 4.73372174e-06, 0, 0, 0.0009079, 0.00001242, 2.4954189e-06, 2.24189267e-06, 3604.22831066, 0, 13589.7529666, 3.81949853e-06, 0.00001144, 3.64231209e-06, 0, 0, 0.00054202, 0.00002136, 7.38183236e-06, 1.64810407e-06, 2790.08710492, 0, 15307.30744853, 2.86986497e-06, 7.85762918e-06, 4.66948185e-06, 0, 0, 0.00031842, 0.00001523, 7.34216482e-06, 1.94081414e-06, 2072.95575498, 0, 17241.45093212, 1.01418675e-06, 0.00001233, 4.22668331e-06, 0, 0, 0.00019022, 0.00002556, 7.93489313e-06, 2.14766225e-06, 1490.67780553, 0, 19419.37287627, 0, 6.90526454e-06, 5.59628017e-06, 0, 0, 0.00011275, 0.00001365, 0, 2.12849193e-06, 1045.08179528, 0, 21871.63512954, 0, 0.00001068, 4.00328984e-06, 0, 0, 0.00006976, 0.00003565, 9.2457521e-06, 2.37345367e-06, 718.30695237, 0, 24632.57468396, 0, 0, 6.63801794e-06, 0, 0, 0.00003995, 0.0000146, 8.68357539e-07, 2.80048124e-06, 485.94184309, 0, 27740.74817348, 0, 0, 4.21889155e-06, 0, 0, 0.00002606, 0.0000281, 0, 5.10561138e-06, 324.43982837, 0, 31239.41965574, 0, 0, 8.95146941e-06, 0, 0, 0.00002009, 0.00001626, 5.90160745e-06, 5.70938772e-06, 214.14846879, 0, 35177.09183824, 0, 0, 6.51078925e-06, 0, 0, 0.00001499, 0.00006215, 7.14422726e-06, 9.82275586e-06, 139.8959837, 0, 39608.07868348, 0, 2.96462876e-06, 0.00001294, 0, 0, 0, 0.00006985, 0.00001634, 0.00001096, 90.51116625, 0, 44593.11385704, 0, 0, 4.59762833e-06, 0, 0, 0.00001674, 0.00010369, 0.00001451, 0.00001391, 58.02128853, 0, 50199.98419121, 0, 0, 0.00001585, 0, 0, 7.63416487e-06, 0.00006482, 0.00001192, 8.12187502e-06, 36.86130085, 0, 56504.16940918, 0, 0, 8.03615838e-06, 0, 0, 0.00001226, 0.0000866, 9.04777118e-06, 6.25380941e-06, 23.21223163, 0, 63589.45762927, 7.9687971e-06, 0, 0.0000284, 0, 0, 0.00001676, 0.00006276, 0.00001147, 4.11332155e-06, 14.48989373, 0, 71548.48899875, 0.00001584, 0, 7.72640572e-06, 0, 0, 0.00003231, 0.00008055, 1.88818879e-06, 9.02025405e-06, 8.96681014, 6.2022775e-06, 80483.15491891, 9.28823956e-06, 0, 0.00003197, 0, 0, 0.00001733, 0.00006827, 0.00002693, 3.73591046e-06, 5.50103804, 0.00001665, 90504.74459564, 0.0000295, 0, 0.00001831, 0, 0, 0.00004714, 0.00004466, 0, 0, 3.34578913, 0.00003126, 101733.67988411, 9.06517965e-06, 0, 0.00003913, 0, 0, 0.00002937, 0, 0.0000452, 0, 2.01745984, 0.00001523, 114298.60827335, 0.00002079, 0, 0.00001559, 0, 0, 0.00004548, 0, 0, 0.0000157, 1.20602636, 0.00005229, 128334.5257819, 0.00001008, 0, 0.00003894, 0, 0, 0.00006507, 0.00003129, 0, 0, 0.71475877, 0.00006221, 143979.46953367, 0.00004822, 0, 0.00003505, 0, 0, 0.00010544, 0, 0, 0, 0.4199596, 0.00005764, 161369.14780663, 0.00003109, 0, 0.00005952, 0, 0, 0.00006926, 0, 0, 0, 0.24461762, 0.00009209, 180628.66235738, 0.00005122, 0, 0.00005555, 0, 0, 0.00013822, 0, 0, 0.00003016, 0.14127466, 0.00032966, 201860.23590235, 0.00001771, 0.00007507, 0.0000659, 0, 0.00001827, 0.0000699, 0, 0, 0, 0.08083119, 0.000422, 225125.62406274, 0.0000522, 0.00016325, 0.00003918, 0, 0.00004052, 0.00017595, 0, 0, 0, 0.04596142, 0.00075568, 250421.75628688, 0, 0.00059968, 0.00009026, 0, 0.00007921, 0.00003779, 0, 0, 0, 0.0258909, 0.00137404, 277648.28360901, 0.00016112, 0.00148384, 0.00002209, 0, 0.00014717, 0.00023742, 0, 0, 0, 0.01435088, 0.00320091, 306566.41433386, 0.00035002, 0.00368743, 0.00016468, 0, 0.00029521, 0.0002279, 0.00034184, 0.00035215, 0, 0.007905, 0.00555587, 336750.1178498, 0.00083138, 0.00731905, 0.00018107, 0, 0.00049909, 0.00063003, 0.00113749, 0.00079701, 0, 0.00406612, 0.01013644, 367534.02230089, 0.00149183, 0.01280565, 0.00038933, 0.00010773, 0.00093653, 0.00089506, 0.00266584, 0.00180335, 0, 0.00240787, 0.01798858, 397967.48773848, 0.00300654, 0.02468916, 0.00062007, 0.00037338, 0.00170129, 0.00167188, 0.00555776, 0.00294593, 0, 0.00134883, 0.03867103, 426791.11252668, 0.00593903, 0.04995222, 0.00133358, 0.00117719, 0.00309686, 0.0029207, 0.01159994, 0.00646669, 0.0001436, 0.0007878, 0.06812319, 452458.41036068, 0.01151879, 0.09187486, 0.00230714, 0.00239454, 0.00560347, 0.00558227, 0.021934, 0.01169905, 0, 0.00027728, 0.13599879, 473227.26887974, 0.02053845, 0.17052994, 0.00415396, 0.00472303, 0.00992154, 0.01068754, 0.04121365, 0.02294491, 0.00009527, 0.00037748, 0.26544919, 487336.8616595, 0.04281941, 0.34086412, 0.00811892, 0.00904485, 0.02015852, 0.02005763, 0.08152433, 0.04482501, 0.00021948, 0, 0.57897909, 493261.44517459, 0.08770988, 0.69240817, 0.01688271, 0.01902149, 0.04045817, 0.04100188, 0.16623097, 0.09208804, 0, 0, 1.03471072, 489995.36235587, 0.15836825, 1.24701447, 0.03008282, 0.03573023, 0.0729175, 0.07557133, 0.31363692, 0.17684684, 0, 0.00004438, 1.92856643, 477289.51652621, 0.29068241, 2.23869064, 0.05496319, 0.06644219, 0.13394181, 0.13493929, 0.58552276, 0.3434318, 0, 0, 3.49103071, 455755.76937232, 0.53097936, 4.12884556, 0.10067457, 0.12282225, 0.24570715, 0.25602879, 1.09793123, 0.62847415, 0.00024151, 0, 7.66907204, 426790.00219807, 1.1315821, 8.99507136, 0.21205825, 0.27446317, 0.51699534, 0.55031814, 2.27988133, 1.2553926, 0, 0, 15.37893301, 392334.56640771, 2.4412716, 18.52252024, 0.45170641, 0.56557108, 1.09566627, 1.18221476, 4.74734073, 2.51706656, 0.00035251, 0, 25.61530693, 354555.76576772, 4.50443036, 33.30195516, 0.87294224, 1.01632089, 2.0403471, 2.0496999, 8.83056261, 4.58802833, 0.00103361, 0, 50.13012825, 315526.67248304, 9.18401679, 62.69989599, 1.84486361, 2.05837273, 4.10623997, 4.11805932, 17.64480178, 8.9983734, 0.00057405, 0, 89.91346667, 277010.83797187, 18.14978419, 112.25805877, 3.65987484, 4.03356713, 8.00164902, 7.88189803, 34.07812578, 17.2003973, 0.00051431, 0, 176.61947847, 240340.77416104, 35.29689367, 202.32807716, 6.95994182, 7.97947577, 15.29077529, 16.32090971, 66.01886605, 34.81329291, 0.00078991, 0, 303.53648794, 206409.83879402, 67.65272811, 347.87684901, 13.85739496, 15.27589999, 29.21076186, 32.43009725, 125.24311516, 68.17076568, 0.00006499, 0, 523.74569908, 175715.96661418, 133.21577481, 594.71883092, 27.42618102, 31.85559082, 57.92198696, 66.47047224, 245.0406616, 132.52607516, 0.0004258, 0, 780.76782094, 148460.26518867, 222.58259004, 880.19475219, 46.57036493, 55.50164946, 98.37992274, 114.99053014, 428.40544385, 232.54364127, 0, 0, 1218.7778229, 124598.6103856, 410.87075644, 1366.48511125, 90.59040882, 107.41987923, 184.49835406, 220.72793522, 766.97122857, 423.71677485, 0, 0, 1880.26235905, 103959.46134658, 736.89201312, 2115.1128539, 172.80982272, 209.53499348, 332.78812664, 423.77416244, 1346.92040033, 801.54443235, 0, 0, 2692.90494654, 86286.37885392, 1218.86641421, 2986.36728888, 309.22245653, 369.17339505, 558.54852675, 776.61315154, 2241.37357747, 1354.68417715, 0, 0, 3814.86882869, 71274.69147303, 1930.51972132, 4149.20670838, 513.12518489, 627.82126646, 883.4026017, 1384.2854157, 3570.25930248, 2194.14216107, 0, 0, 5560.85843885, 58605.63637927, 3181.18422011, 5963.56590512, 949.00144816, 1169.39279678, 1443.1791784, 2485.0583473, 5719.64296537, 3687.69688833, 0, 0, 7644.65408973, 47990.01879263, 4989.66623972, 8241.14934149, 1735.94497636, 1983.09891179, 2279.69821829, 4271.38242682, 8822.54343895, 5882.91753209, 0, 0, 10390.91719066, 39140.90784419, 7893.98447507, 11305.41248519, 2953.93621901, 3208.16846361, 3587.39958144, 6942.15182328, 13445.03538225, 9163.86753866, 0, 0, 15216.64888419, 31832.88724429, 12511.42801271, 16083.2585271, 5197.15949582, 5459.3355251, 5699.3796258, 11929.86972182, 20579.52988031, 14389.72237903, 0, 0, 20346.06829247, 25890.4351607, 18335.47880524, 21973.68237922, 8620.65143089, 8442.37840934, 8394.04987516, 18415.70946506, 29918.27717237, 20973.69151732, 0, 0, 27624.62824038, 21207.59632586, 26554.45984821, 29845.53113159, 13934.48851952, 12606.20934775, 12189.2479224, 27646.21837947, 42852.62687927, 30401.56856628, 0, 0, 37104.3925841, 17812.76513281, 37287.00455432, 40414.94350493, 20682.93453015, 18427.8403352, 17030.98899472, 40369.48498355, 60180.69406035, 42603.74918128, 0, 0, 49403.43531625, 15681.9421436, 48644.63311158, 53289.89608753, 28180.15956348, 25694.95396319, 22493.07270811, 56406.05237058, 80581.716569, 58776.18441071, 0, 0.00654518, 70095.03303922, 16229.8357019, 70370.5447554, 75697.50761684, 42566.66064957, 38421.5804856, 32309.72763779, 82878.12054798, 114855.21690036, 83498.44176756, 0.00255168, 0.01894924, 92268.79390415, 19524.61864111, 95015.05939341, 100345.69681928, 60178.72850206, 53404.12618431, 43592.03283278, 112265.13658718, 157478.0923, 112314.91795527, 0.00366588, 0.05205452, 124435.59515119, 27612.99289701, 125981.71039233, 132307.98684876, 82582.59649423, 73895.46173614, 57911.7683488, 153771.42917618, 211409.76612131, 150058.11875327, 0.01067804, 0.11107304, 158104.23882742, 38052.50118934, 155033.7212158, 162531.99371489, 101430.04954781, 91723.07369489, 71567.70923364, 193605.42921344, 265247.87327095, 194398.6307087, 0.01099226, 0.23340087, 208205.46965899, 61432.80657841, 205456.01078883, 218020.84690934, 135954.08317421, 123262.97819235, 94638.52490455, 256757.58642868, 346494.69945549, 256070.0951917, 0.01327006, 0.46376461, 281717.47083414, 102484.90707679, 283836.69943699, 297234.75215091, 184753.19498061, 168590.37904045, 129233.89714778, 358704.86492389, 462494.66707971, 342204.70047955, 0.01016524, 0.92081997, 376836.12589862, 169131.30895087, 388930.22289783, 403395.4415051, 249280.33108828, 221488.06710519, 176871.51982665, 483088.44126256, 623888.53111355, 461852.27234754, 0.02005746, 1.80347544, 511666.97979077, 270398.83430127, 517209.75122705, 551146.41762208, 328484.72561869, 302889.14341247, 234081.83532892, 651894.32170713, 827853.7232406, 604250.19082904, 0.01731938, 3.59473485, 703493.07706975, 424524.81754771, 695050.43296251, 745016.04560688, 452524.69323542, 408551.46689662, 318166.00971512, 889779.44163768, 1108485.15540879, 839561.56158748, 0.01986635, 7.11172972, 939888.55092016, 651600.53027346, 955105.96373362, 1008015.3856989, 621296.87852454, 544087.44058052, 433113.78655125, 1203719.22134031, 1536493.96177724, 1137879.37804099, 0.02720567, 14.0663821, 1208821.17163974, 906782.14673104, 1257990.63307435, 1313985.68727009, 811472.40729395, 695834.99362118, 577726.52276587, 1562470.58740019, 2008426.14441814, 1478164.05261052, 0.04969386, 27.97449453, 1675644.40551546, 1325684.73919883, 1735922.00141615, 1786885.18068007, 1127982.95628237, 950186.1305416, 799238.75607471, 2131648.56124591, 2731804.2062668, 1998022.30645471, 0.07672454, 55.22003021, 2244725.17444066, 1821091.56447156, 2247220.29885956, 2423006.68066188, 1469425.09998245, 1262291.56896417, 1037881.25959464, 2827333.37998341, 3611661.04299921, 2633796.7535235, 0.12119996, 107.66091168, 2906709.34678089, 2449510.65561756, 2952098.67158503, 3199234.77903282, 1907953.53825024, 1685971.61349771, 1369617.88599269, 3650223.00512727, 4738658.07358305, 3381846.39955074, 0.2349558, 214.74172046, 3974746.32495928, 3512273.81621836, 4278783.61755226, 4391293.30001398, 2802414.69657314, 2302689.97423803, 1958525.32823069, 5019715.24530596, 6530387.68862364, 4581979.83285836, 0.41392746, 428.89317989, 5736883.77325683, 5023655.66034893, 5790941.47123096, 6099592.8269074, 3855120.32624461, 3259626.50651127, 2659253.19301739, 7203888.65527265, 9034272.06431425, 6464132.7832732, 0.81025046, 838.85861403, 7431788.7967373, 6919459.35260311, 7934095.06311153, 8159155.61600027, 5357286.1023839, 4307311.96236156, 3656578.91111336, 9390150.97228902, 12302189.01800773, 8683949.69994465, 1.51998026, 1646.0712153, 10235039.96798785, 9457072.38121275, 10627689.07002106, 11045604.78390486, 7063000.51778886, 5897459.50481315, 4912955.36982989, 13172242.63967287, 16761941.85855336, 11944154.56868635, 2.97884065, 3223.22493017, 13432389.32390384, 12430015.61128897, 14194451.66712877, 14484677.64596281, 9513184.20783649, 7770971.01891774, 6612963.16177152, 17242136.19282264, 22429768.61259142, 15975406.58774081, 5.73783119, 6411.21667866, 17759621.11442146, 16381351.73897, 18514803.49466322, 19119508.09818694, 12659641.91004088, 10336773.50567668, 8624611.95822755, 22792838.62335711, 29600130.32258716, 21599981.359332, 11.10553222, 12735.94440331, 23827242.21440469, 21661330.28951951, 24378054.48567673, 25331815.39913674, 16787242.7419357, 13776767.94854832, 11360566.54890397, 30264749.21425445, 39206490.74343798, 29573782.38531723, 21.39713242, 24851.90197865, 32741470.3348749, 29527970.09962671, 32776381.43948571, 34344992.38288016, 22068204.48384048, 18613944.54665855, 15190151.24789516, 41645600.46035212, 52815720.16928751, 38988455.27739758, 41.44100286, 48669.46929824, 43456633.72053333, 38552022.47062275, 43423052.4868163, 45380303.20794479, 28985748.40474013, 24697802.65893549, 20193491.96976492, 55276116.82126442, 70059315.7063974, 52503111.49599912, 80.05893239, 94641.0460601, 58017980.12369667, 51888598.17176614, 58398417.79651085, 61623903.03347099, 38375467.31928742, 33705650.28174309, 26874396.47710833, 74172799.17692499, 94309953.4355568, 70617074.47317019, 155.1456068, 180603.02006985, 77162318.75067097, 68306748.51547684, 77527496.46716143, 81421177.57439694, 50196725.16073133, 45907489.72006904, 35344456.32881003, 98516162.45434283, 127436577.09141569, 91817513.13714676, 300.65750712, 343443.9697776, 105330981.8667376, 92170492.34141624, 100500982.26708592, 108326796.2438293, 64462551.25747406, 60193514.50381986, 45753710.30612082, 132005010.91982752, 170227306.82412043, 123917140.36270706, 582.344316, 649030.80485455, 144190153.04240641, 123981618.641054, 133887798.77104129, 145262883.43120223, 86646515.1363081, 80716267.06574172, 60949159.57822782, 180308885.58043122, 226140117.59066874, 174330727.95899472, 1124.28430541, 1197923.64017148, 191634512.71291131, 168467983.07627687, 182690471.74343029, 195311968.47787145, 116679019.68680321, 106935204.10362458, 83283166.47738604, 241938857.52164772, 303069593.6895442, 232051844.62218204, 2153.45289729, 2163428.8659767, 252938124.4144589, 227560832.8655715, 249570444.07573673, 268123496.27124065, 156339478.02820808, 144460757.1358318, 112429266.13690938, 321778087.95028317, 405688791.61676925, 305241113.83719373] + } + diff --git a/tests/testthat/_snaps/project_n.md b/tests/testthat/_snaps/project_n.md new file mode 100644 index 000000000..09106d93b --- /dev/null +++ b/tests/testthat/_snaps/project_n.md @@ -0,0 +1,36 @@ +# project_n with nonzero diffusion produces a different result + + { + "type": "double", + "attributes": { + "dim": { + "type": "integer", + "attributes": {}, + "value": [1, 100] + }, + "dimnames": { + "type": "list", + "attributes": { + "names": { + "type": "character", + "attributes": {}, + "value": ["sp", "w"] + } + }, + "value": [ + { + "type": "character", + "attributes": {}, + "value": ["Saithe"] + }, + { + "type": "character", + "attributes": {}, + "value": ["0.001", "0.00119", "0.00142", "0.0017", "0.00203", "0.00242", "0.00289", "0.00345", "0.00411", "0.00491", "0.00585", "0.00699", "0.00834", "0.00995", "0.0119", "0.0142", "0.0169", "0.0202", "0.0241", "0.0287", "0.0343", "0.0409", "0.0488", "0.0582", "0.0695", "0.0829", "0.099", "0.118", "0.141", "0.168", "0.201", "0.239", "0.286", "0.341", "0.407", "0.486", "0.579", "0.691", "0.825", "0.985", "1.17", "1.4", "1.67", "2", "2.38", "2.84", "3.39", "4.05", "4.83", "5.76", "6.88", "8.21", "9.8", "11.7", "13.9", "16.6", "19.9", "23.7", "28.3", "33.8", "40.3", "48.1", "57.4", "68.4", "81.7", "97.5", "116", "139", "166", "198", "236", "281", "336", "401", "478", "571", "681", "812", "970", "1160", "1380", "1650", "1970", "2350", "2800", "3340", "3990", "4760", "5680", "6770", "8080", "9650", "11500", "13700", "16400", "19600", "23300", "27900", "33200", "39700"] + } + ] + } + }, + "value": [34.47062159, 17.31184266, 11.60313231, 8.75859551, 7.06137563, 5.93950144, 5.14816892, 4.56532548, 4.1235451, 3.78279304, 3.51806219, 3.31320879, 3.1576526, 3.04451246, 2.96951082, 2.93031401, 2.92613204, 2.95747964, 3.02604103, 3.13460291, 3.28703176, 3.48827782, 3.74439034, 4.06252899, 4.45095437, 4.91897777, 5.47684609, 6.13553348, 6.9064076, 7.8007363, 8.82900327, 10.00001063, 11.31976733, 12.7901972, 14.40775271, 16.16208846, 18.03502371, 20.00008819, 22.02297029, 24.06313423, 26.07670953, 28.02047499, 29.85640415, 31.55591752, 33.10285441, 34.49436466, 35.73944779, 36.85557195, 37.86439026, 38.78775415, 39.6449248, 40.45128873, 41.21832112, 41.95425714, 42.66496011, 43.35467778, 44.02658427, 44.68312877, 45.32625079, 45.95751525, 46.57820148, 47.1893654, 47.79188515, 48.38649549, 48.97381457, 49.5543649, 50.1285902, 50.69686887, 51.25952498, 51.81683731, 52.36904667, 52.91636207, 53.45896576, 53.99701745, 54.53065785, 55.06001158, 55.58518967, 56.10629162, 56.62340711, 57.13661752, 57.64599711, 58.1516141, 58.65353151, 59.15180791, 59.64649806, 60.13765342, 60.62532261, 61.10955177, 61.59038488, 62.06786404, 62.5420297, 63.01292085, 63.48057517, 63.94502921, 64.40631845, 64.86447744, 65.31953991, 65.77153878, 66.22050627, 66.66647389] + } + diff --git a/tests/testthat/helper.R b/tests/testthat/helper.R index 650bcd483..c3b756801 100644 --- a/tests/testthat/helper.R +++ b/tests/testthat/helper.R @@ -26,7 +26,7 @@ example_params <- function() { # Give diffusion to one species n <- params@species_params$n[1] d <- 0.1 * params@w^(n + 1) - diffusion(params)[1, ] <- d + ext_diffusion(params)[1, ] <- d params } diff --git a/tests/testthat/test-ArraySpeciesBySize.R b/tests/testthat/test-ArraySpeciesBySize.R index 2c5e979bd..4090e08f2 100644 --- a/tests/testthat/test-ArraySpeciesBySize.R +++ b/tests/testthat/test-ArraySpeciesBySize.R @@ -38,7 +38,7 @@ test_that("Rate functions return ArraySpeciesBySize", { expect_true(is.ArraySpeciesBySize(getExtEncounter(params))) expect_true(is.ArraySpeciesBySize(getMaturityProportion(params))) expect_true(is.ArraySpeciesBySize(getReproductionProportion(params))) - expect_true(is.ArraySpeciesBySize(diffusion(params))) + expect_true(is.ArraySpeciesBySize(ext_diffusion(params))) }) test_that("All rate functions have consistent dimnames", { diff --git a/tests/testthat/test-analytic_results.R b/tests/testthat/test-analytic_results.R index 3f313f5eb..8024446be 100644 --- a/tests/testthat/test-analytic_results.R +++ b/tests/testthat/test-analytic_results.R @@ -219,3 +219,65 @@ test_that("getFeedingLevel approximates analytic result", { # # so we only test the others # # expect_lt(max(relative_error[1:(no_w - 1)]), 0.02) # }) + + +# mizerDiffusion ---- +# The setup (p, n0, n_pp, sp, sigma, beta, gamma, q, n, lm2) is inherited from +# the top of this file. We only need a few extra variables: +alpha_sp <- p@species_params$alpha[sp] +lambda <- p@resource_params$lambda +lm3 <- lambda - 3 # = lambda - 3, used in D(w) analytic formula +Dx <- p@w[2] / p@w[1] - 1 +dx <- log(p@w[2] / p@w[1]) + +test_that("mizerDiffusion gives power law in w^(n+1) when feeding on resource only", { + p_d <- p + p_d@use_predation_diffusion <- TRUE + d <- getDiffusion(p_d, n0, n_pp)[sp, ] + # D(w) is proportional to w^(n+1), so D(w) / w^(n+1) should be constant + # (eq-dw2 in vignette predation_diffusion.qmd) + d_coeff <- d / p_d@w^(n + 1) + expect_equal(d_coeff, rep(d_coeff[1], length(d_coeff)), + tolerance = 1e-3, ignore_attr = TRUE) +}) + +test_that("mizerDiffusion matches analytic formula (eq-dw2)", { + p_d <- p + p_d@use_predation_diffusion <- TRUE + d <- getDiffusion(p_d, n0, n_pp)[sp, ] + f0 <- getFeedingLevel(p_d, n0, n_pp)[sp, 1] + # Coefficient from the Gaussian integral (eq-dw2), with Dx/dx correction + # for the log-spaced Riemann sum used in the FFT convolution + diffusion_analytic <- (1 - f0) * p_d@resource_params$kappa * + exp(lm3^2 * sigma^2 / 2) * + beta^lm3 * sqrt(2 * pi) * sigma * + alpha_sp^2 * gamma * Dx / dx + d_coeff <- d / p_d@w^(n + 1) + # The diffusion integral uses w_p^2 weighting, making the Riemann sum + # slightly less accurate than for the encounter rate; 2% tolerance is adequate. + expect_equal(d_coeff[1], diffusion_analytic, tolerance = 2e-2, ignore_attr = TRUE) +}) + +test_that("mizerDiffusion satisfies D(w) = A(w) * w * alpha/beta * exp((5-2*lambda)*sigma^2/2) (eq-dw3)", { + # This tests the analytic relationship between diffusion and assimilation + # rates derived in the vignette predation_diffusion.qmd (eq-dw3). + p_d <- p + p_d@use_predation_diffusion <- TRUE + d <- getDiffusion(p_d, n0, n_pp)[sp, ] + E <- getEncounter(p_d, n0, n_pp)[sp, ] + f <- getFeedingLevel(p_d, n0, n_pp)[sp, ] + # Assimilation rate A(w) = alpha * (1 - f(w)) * E(w) + A <- alpha_sp * (1 - f) * E + ratio <- d / (A * p_d@w) + expected_ratio <- alpha_sp / beta * exp((5 - 2 * lambda) * sigma^2 / 2) + expect_equal(ratio, rep(expected_ratio, length(ratio)), + tolerance = 1e-2, ignore_attr = TRUE) +}) + +test_that("mizerDiffusion is zero when use_predation_diffusion is FALSE", { + # Default is FALSE, so only ext_diffusion (zero by default) contributes + p_d <- p + expect_false(p_d@use_predation_diffusion) + d <- getDiffusion(p_d, n0, n_pp)[sp, ] + expect_equal(unname(d), rep(0, length(d))) +}) diff --git a/tests/testthat/test-analytic_transport.R b/tests/testthat/test-analytic_transport.R index 5da021a4c..3a3a625f7 100644 --- a/tests/testthat/test-analytic_transport.R +++ b/tests/testthat/test-analytic_transport.R @@ -98,7 +98,7 @@ test_that("Exact steady state is maintained", { params <- setRateFunction(params, "EGrowth", "start_growth") params <- setRateFunction(params, "Mort", "start_mort") params <- setRateFunction(params, "RDD", "constant_rdd") - params@diffusion[1, ] <- K * params@w^(p + 1) + params@ext_diffusion[1, ] <- K * params@w^(p + 1) params <- setResource(params, resource_dynamics = "resource_constant") initialNResource(params) <- 0 @@ -133,7 +133,7 @@ test_that("Exact time-dependent solution is followed", { params <- setRateFunction(params, "EGrowth", "start_growth") params <- setRateFunction(params, "Mort", "start_mort") params <- setRateFunction(params, "RDD", "time_dep_rdd") - params@diffusion[1, ] <- K * params@w^(p + 1) + params@ext_diffusion[1, ] <- K * params@w^(p + 1) params <- setResource(params, resource_dynamics = "resource_constant") initialNResource(params) <- 0 diff --git a/tests/testthat/test-diffusion.R b/tests/testthat/test-diffusion.R new file mode 100644 index 000000000..4d8ea7cc2 --- /dev/null +++ b/tests/testthat/test-diffusion.R @@ -0,0 +1,98 @@ +# getDiffusion ---- +test_that("getDiffusion returns ArraySpeciesBySize with correct dimensions", { + params <- newSingleSpeciesParams() + d <- getDiffusion(params) + expect_true(is.ArraySpeciesBySize(d)) + expect_identical(dim(d), dim(params@initial_n)) +}) + +test_that("getDiffusion includes ext_diffusion", { + params <- newSingleSpeciesParams() + d_base <- getDiffusion(params) + # Adding a constant to ext_diffusion should shift getDiffusion by the same amount + params@ext_diffusion[] <- 1 + d_with_ext <- getDiffusion(params) + expect_equal(d_with_ext, d_base + 1, ignore_attr = TRUE) +}) + +test_that("getDiffusion dispatches via rates_funcs", { + params <- newSingleSpeciesParams() + e <- globalenv() + e$constant_diffusion <- function(params, n, n_pp, n_other, t, feeding_level, ...) { + array(42, dim = dim(params@initial_n), dimnames = dimnames(params@initial_n)) + } + params@rates_funcs$Diffusion <- "constant_diffusion" + d <- getDiffusion(params) + expect_true(all(d == 42)) +}) + +test_that("r$diffusion is included in getRates output", { + params <- newSingleSpeciesParams() + r <- getRates(params) + expect_true("diffusion" %in% names(r)) + expect_identical(dim(r$diffusion), dim(params@initial_n)) +}) + +test_that("r$diffusion matches getDiffusion", { + params <- newSingleSpeciesParams() + r <- getRates(params) + expect_equal(r$diffusion, getDiffusion(params), ignore_attr = TRUE) +}) + +# mizerDiffusion behaviour ---- + +test_that("mizerDiffusion accepts pre-computed feeding_level and gives same result", { + params <- newSingleSpeciesParams() + params@use_predation_diffusion <- TRUE + n <- initialN(params) + n_pp <- initialNResource(params) + fl <- getFeedingLevel(params, n = n, n_pp = n_pp) + d_auto <- getDiffusion(params, n, n_pp) + d_given <- mizerDiffusion(params, n = n, n_pp = n_pp, + n_other = initialNOther(params), + t = 0, feeding_level = fl) + expect_equal(d_auto, d_given, ignore_attr = TRUE) +}) + +test_that("mizerDiffusion is zero when feeding level is 1 everywhere", { + # D(w) = (1 - f(w)) * ... so f = 1 gives D = ext_diffusion + params <- newSingleSpeciesParams() + params@use_predation_diffusion <- TRUE + n <- initialN(params) + n_pp <- initialNResource(params) + fl_one <- matrix(1, nrow = nrow(n), ncol = ncol(n), dimnames = dimnames(n)) + d <- mizerDiffusion(params, n = n, n_pp = n_pp, + n_other = initialNOther(params), + t = 0, feeding_level = fl_one) + expect_equal(d, params@ext_diffusion, ignore_attr = TRUE) +}) + +test_that("mizerDiffusion scales as alpha^2", { + params <- newSingleSpeciesParams() + params@use_predation_diffusion <- TRUE + n <- initialN(params) + n_pp <- initialNResource(params) + d1 <- getDiffusion(params, n, n_pp) + # Double alpha → diffusion should quadruple (alpha enters as alpha^2) + params2 <- params + params2@species_params$alpha <- params@species_params$alpha * 2 + # Keep search_vol unchanged to isolate the alpha^2 factor + params2@search_vol[] <- params@search_vol + d2 <- getDiffusion(params2, n, n_pp) + expect_equal(d2, d1 * 4, tolerance = 1e-12, ignore_attr = TRUE) +}) + +test_that("mizerDiffusion increases when fish prey are present", { + # With use_predation_diffusion = TRUE, adding fish as prey increases D. + # This exercises the params@interaction %*% n term. + params <- newMultispeciesParams(NS_species_params_gears, inter, + no_w = 30, info_level = 0) + params@use_predation_diffusion <- TRUE + n_zero <- initialN(params) + n_zero[] <- 0 + n_pp <- initialNResource(params) + d_no_fish <- getDiffusion(params, n_zero, n_pp) + d_with_fish <- getDiffusion(params, initialN(params), n_pp) + # At least one species should have larger diffusion with fish present + expect_true(any(d_with_fish > d_no_fish)) +}) diff --git a/tests/testthat/test-getFlux.R b/tests/testthat/test-getFlux.R index 946208a95..2bd31379c 100644 --- a/tests/testthat/test-getFlux.R +++ b/tests/testthat/test-getFlux.R @@ -9,7 +9,7 @@ test_that("getFlux works correctly", { t <- 0 g <- getEGrowth(params, n = n, t = t) - d <- params@diffusion + d <- params@ext_diffusion dw <- params@dw rdd <- getRDD(params, n = n, t = t) diff --git a/tests/testthat/test-getRequiredRDD.R b/tests/testthat/test-getRequiredRDD.R index ddf68ff6f..3326b04f9 100644 --- a/tests/testthat/test-getRequiredRDD.R +++ b/tests/testthat/test-getRequiredRDD.R @@ -24,9 +24,9 @@ test_that("getRequiredRDD handles diffusion", { params <- newSingleSpeciesParams() # Add diffusion - diffusion <- params@diffusion - diffusion[] <- 1e9 * params@w - params <- setDiffusion(params, diffusion = diffusion) + ext_diff <- params@ext_diffusion + ext_diff[] <- 1e9 * params@w + params <- setExtDiffusion(params, ext_diffusion = ext_diff) # Update initial_n to be the steady state solution with this diffusion # We need to recalculate it using get_steady_state_n diff --git a/tests/testthat/test-get_steady_state_n.R b/tests/testthat/test-get_steady_state_n.R index e3714e9a5..4d6bbf27b 100644 --- a/tests/testthat/test-get_steady_state_n.R +++ b/tests/testthat/test-get_steady_state_n.R @@ -9,7 +9,7 @@ test_that("get_steady_state_n works with no diffusion", { N0_vec <- rep(100, no_sp) # Zero diffusion - params@diffusion[] <- 0 + params@ext_diffusion[] <- 0 n_calc <- mizer:::get_steady_state_n(params, growth, mort, N0_vec) @@ -53,7 +53,7 @@ test_that("get_steady_state_n works with diffusion", { } # Compare with no diffusion - params@diffusion[] <- 0 + params@ext_diffusion[] <- 0 n_nodiff <- mizer:::get_steady_state_n(params, growth, mort, N0_vec) expect_false(isTRUE(all.equal(n_calc, n_nodiff))) }) diff --git a/tests/testthat/test-project.R b/tests/testthat/test-project.R index 0d705f67d..c6a4358ff 100644 --- a/tests/testthat/test-project.R +++ b/tests/testthat/test-project.R @@ -388,8 +388,38 @@ test_that("t_max less than effort array duration uses effort times", { test_that("project does not change the params object", { params <- NS_params - params@diffusion[] <- 1 + params@ext_diffusion[] <- 1 old_params <- unserialize(serialize(params, NULL)) sim <- project(params, t_max = 1) expect_identical(params, old_params) -}) \ No newline at end of file +}) + + +# predation diffusion ---- +test_that("Simulation gives same numerical results with predation diffusion", { + params_d <- NS_params + params_d@use_predation_diffusion <- TRUE + sim_d <- project(params_d, t_max = 1) + expect_snapshot(sim_d@n[2, 3, ]) + expect_snapshot(sim_d@n_pp[2, ]) +}) + +test_that("Final result the same when called with sim or params, predation diffusion on", { + params_d <- NS_params + params_d@use_predation_diffusion <- TRUE + sim_d <- project(params_d, t_max = 1) + params_d@initial_n[] <- sim_d@n[2, , ] + params_d@initial_n_pp[] <- sim_d@n_pp[2, ] + params_d@initial_n_other <- sim_d@n_other[2, ] + sim1 <- project(params_d, t_max = 1) + sim2 <- project(sim_d, t_max = 1) + expect_identical(sim1@n[2, 3, ], sim2@n[3, 3, ]) +}) + +test_that("Predation diffusion changes simulation trajectory", { + params_d <- NS_params + params_d@use_predation_diffusion <- TRUE + sim_d <- project(params_d, t_max = 1) + sim_base <- project(NS_params, t_max = 1) + expect_false(identical(sim_d@n, sim_base@n)) +}) diff --git a/tests/testthat/test-project_methods.R b/tests/testthat/test-project_methods.R index 033f37a8d..daa154d5f 100644 --- a/tests/testthat/test-project_methods.R +++ b/tests/testthat/test-project_methods.R @@ -42,10 +42,10 @@ params2@initial_effort <- params2@initial_effort / 2 # getRates ---- test_that("getRates works", { r <- getRates(params) - expect_identical(names(r), - c("encounter", "feeding_level", "e", "e_repro", - "e_growth", "pred_rate", "pred_mort", "f_mort", - "mort", "rdi", "rdd", "resource_mort")) + expect_identical(names(r), + c("encounter", "feeding_level", "e", "e_repro", + "e_growth", "diffusion", "pred_rate", "pred_mort", + "f_mort", "mort", "rdi", "rdd", "resource_mort")) # test that the optional parameters take the correct defaults expect_identical(r, getRates(params, n = params@initial_n, @@ -642,8 +642,8 @@ test_that("mizerRates returns the standard rate list from registered functions", rates_fns = rates_fns) expect_named(r, c("encounter", "feeding_level", "e", "e_repro", "e_growth", - "pred_rate", "pred_mort", "f_mort", "mort", "rdi", - "rdd", "resource_mort")) + "diffusion", "pred_rate", "pred_mort", "f_mort", "mort", + "rdi", "rdd", "resource_mort")) expect_identical(r$encounter, rates_fns$Encounter(params, n = params@initial_n, @@ -710,6 +710,7 @@ test_that("project function returns objects of correct dimension when community expect_equal(dim(getEGrowth(params, n, n_pp)), c(1, no_w)) expect_length(getRDI(params, n, n_pp), 1) expect_length(getRDD(params, n, n_pp), 1) + expect_equal(dim(getDiffusion(params, n, n_pp)), c(1, no_w)) # MizerSim functions # time x species x size @@ -723,6 +724,33 @@ test_that("project function returns objects of correct dimension when community expect_equal(dim(getFMortGear(sim)), c(t_max + 1, 1, 1, no_w)) # time x species x size - note drop = TRUE expect_equal(dim(getFMort(sim)), c(t_max + 1, no_w)) - # time x species x size + # time x species x size expect_equal(dim(getFMort(sim, drop = FALSE)), c(t_max + 1, 1, no_w)) }) + + +# getDiffusion with predation diffusion ---- +test_that("getDiffusion snapshot with use_predation_diffusion", { + params_d <- params + params_d@use_predation_diffusion <- TRUE + d <- getDiffusion(params_d, n, n_full) + expect_snapshot_value(drop_params(d), style = 'json2', tolerance = 1e-5) +}) + +test_that("predation diffusion is non-negative and adds to baseline", { + params_d <- params + params_d@use_predation_diffusion <- TRUE + d_pred <- getDiffusion(params_d, n, n_full) + d_base <- getDiffusion(params, n, n_full) + # predation diffusion is non-negative, so total must not fall below baseline + expect_true(all(d_pred >= d_base - .Machine$double.eps)) + # with non-zero abundances it must be strictly larger somewhere + expect_true(any(d_pred > d_base)) +}) + +test_that("getRates diffusion is nonzero with use_predation_diffusion", { + params_d <- params + params_d@use_predation_diffusion <- TRUE + r_d <- getRates(params_d) + expect_true(any(r_d$diffusion > 0)) +}) diff --git a/tests/testthat/test-project_n.R b/tests/testthat/test-project_n.R index ca33255a2..d0ece3b95 100644 --- a/tests/testthat/test-project_n.R +++ b/tests/testthat/test-project_n.R @@ -12,13 +12,16 @@ test_that("project_n follows the documented one-step update", { dimnames = dimnames(params@initial_n)) mort <- matrix(seq(0.05, 0.05 * no_w, by = 0.05), nrow = no_sp, dimnames = dimnames(params@initial_n)) - r <- list(e_growth = e_growth, mort = mort, rdd = 3) + diffusion <- matrix(0, nrow = no_sp, ncol = no_w, + dimnames = dimnames(params@initial_n)) + r <- list(e_growth = e_growth, mort = mort, rdd = 3, diffusion = diffusion) dt <- 0.1 a <- matrix(0, nrow = no_sp, ncol = no_w) b <- matrix(0, nrow = no_sp, ncol = no_w) + c <- matrix(0, nrow = no_sp, ncol = no_w) S <- matrix(0, nrow = no_sp, ncol = no_w) - result <- project_n(params, r, n, dt, a, b, S, idx, w_min_idx_array_ref, + result <- project_n(params, r, n, dt, a, b, c, S, idx, w_min_idx_array_ref, no_sp, no_w) a_expected <- a @@ -41,3 +44,55 @@ test_that("project_n follows the documented one-step update", { expect_equal(result, expected) }) + +test_that("project_n with nonzero diffusion produces a different result", { + params <- newMultispeciesParams(NS_species_params_gears[12, , drop = FALSE], + info_level = 0) + no_sp <- nrow(params@species_params) + no_w <- length(params@w) + idx <- 2:no_w + w_min_idx_array_ref <- (params@w_min_idx - 1) * no_sp + seq_len(no_sp) + + n <- matrix(seq_len(no_w), nrow = no_sp, + dimnames = dimnames(params@initial_n)) + e_growth <- matrix(seq(0.2, 0.2 * no_w, by = 0.2), nrow = no_sp, + dimnames = dimnames(params@initial_n)) + mort <- matrix(seq(0.05, 0.05 * no_w, by = 0.05), nrow = no_sp, + dimnames = dimnames(params@initial_n)) + dt <- 0.1 + a <- b <- c <- S <- matrix(0, nrow = no_sp, ncol = no_w) + + r_zero <- list(e_growth = e_growth, mort = mort, rdd = 3, + diffusion = matrix(0, nrow = no_sp, ncol = no_w, + dimnames = dimnames(params@initial_n))) + r_diff <- list(e_growth = e_growth, mort = mort, rdd = 3, + diffusion = matrix(0.01 * seq_len(no_w), nrow = no_sp, ncol = no_w, + dimnames = dimnames(params@initial_n))) + + result_zero <- project_n(params, r_zero, n, dt, a, b, c, S, idx, + w_min_idx_array_ref, no_sp, no_w) + result_diff <- project_n(params, r_diff, n, dt, a, b, c, S, idx, + w_min_idx_array_ref, no_sp, no_w) + + expect_false(identical(result_zero, result_diff)) + expect_snapshot_value(result_diff, style = 'json2', tolerance = 1e-5) +}) + +test_that("project_n with predation diffusion stays non-negative and finite", { + params_d <- NS_params + params_d@use_predation_diffusion <- TRUE + r <- getRates(params_d) + # confirm predation diffusion is active + expect_true(any(r$diffusion > 0)) + + no_sp <- nrow(params_d@species_params) + no_w <- length(params_d@w) + idx <- 2:no_w + w_min_idx_array_ref <- (params_d@w_min_idx - 1) * no_sp + seq_len(no_sp) + a <- b <- c <- S <- matrix(0, nrow = no_sp, ncol = no_w) + + n_new <- project_n(params_d, r, params_d@initial_n, dt = 0.1, + a, b, c, S, idx, w_min_idx_array_ref, no_sp, no_w) + expect_true(all(is.finite(n_new))) + expect_true(all(n_new >= 0)) +}) diff --git a/tests/testthat/test-setDiffusion.R b/tests/testthat/test-setDiffusion.R deleted file mode 100644 index 970f8ee17..000000000 --- a/tests/testthat/test-setDiffusion.R +++ /dev/null @@ -1,30 +0,0 @@ -test_that("setDiffusion sets and returns the diffusion array", { - params <- NS_params - new <- params@diffusion + 1 - - updated <- setDiffusion(params, diffusion = new) - - expect_true(is.ArraySpeciesBySize(diffusion(updated))) - expect_equal(diffusion(updated), new, ignore_attr = TRUE) - expect_equal(updated@diffusion, new, ignore_attr = TRUE) -}) - -test_that("setDiffusion preserves comments from the supplied array", { - params <- NS_params - new <- params@diffusion - comment(new) <- "custom" - - updated <- setDiffusion(params, diffusion = new) - - expect_identical(comment(updated@diffusion), "custom") -}) - -test_that("setDiffusion validates dimensions and accessors delegate correctly", { - params <- NS_params - - expect_error(setDiffusion(params, diffusion = array(0, dim = c(1, 1)))) - - new <- params@diffusion + 2 - diffusion(params) <- new - expect_identical(params@diffusion, new) -}) diff --git a/tests/testthat/test-setExtDiffusion.R b/tests/testthat/test-setExtDiffusion.R new file mode 100644 index 000000000..1e105bf9e --- /dev/null +++ b/tests/testthat/test-setExtDiffusion.R @@ -0,0 +1,101 @@ +test_that("setExtDiffusion sets and returns the ext_diffusion array", { + params <- newSingleSpeciesParams() + new <- ext_diffusion(params) + 1 + + updated <- setExtDiffusion(params, ext_diffusion = new) + + expect_true(is.ArraySpeciesBySize(ext_diffusion(updated))) + expect_equal(ext_diffusion(updated), new, ignore_attr = TRUE) + expect_equal(updated@ext_diffusion, new, ignore_attr = TRUE) +}) + +test_that("setExtDiffusion preserves comments from the supplied array", { + params <- newSingleSpeciesParams() + new <- ext_diffusion(params) + comment(new) <- "custom" + + updated <- setExtDiffusion(params, ext_diffusion = new) + + expect_identical(comment(updated@ext_diffusion), "custom") +}) + +test_that("setExtDiffusion works", { + params <- NS_params + + # Without ext_diffusion argument, recalculates from D_ext (default 0) and n + p2 <- setExtDiffusion(params) + zero_diffusion <- p2@ext_diffusion + zero_diffusion[] <- 0 + expect_identical(zero_diffusion, p2@ext_diffusion) + + # supplying ext_diffusion + p2 <- setExtDiffusion(params, 3 * params@ext_diffusion + 1) + expect_equal(p2@ext_diffusion, 3 * params@ext_diffusion + 1, + ignore_attr = TRUE) + + # only ext_diffusion changed + p2@ext_diffusion <- params@ext_diffusion + expect_unchanged(p2, params) + + # has updated time_modified + expect_false(identical(params@time_modified, p2@time_modified)) +}) + +test_that("setExtDiffusion uses D_ext species param", { + params <- NS_params + species_params(params)$D_ext <- 0.1 + p2 <- setExtDiffusion(params, reset = TRUE) + expected <- sweep(outer(species_params(params)[["n"]], + w(params), function(x, y) y^(x + 1)), + 1, species_params(params)[["D_ext"]], "*") + expect_equal(p2@ext_diffusion, expected, ignore_attr = TRUE) +}) + +test_that("reset works on ext_diffusion", { + params <- NS_params + # Set a custom ext_diffusion with a comment + custom <- params@ext_diffusion + custom[] <- 1 + comment(custom) <- "custom" + params <- setExtDiffusion(params, ext_diffusion = custom) + expect_identical(comment(params@ext_diffusion), "custom") + + # reset = TRUE ignores comment and recalculates from species params (D_ext=0) + p2 <- setExtDiffusion(params, reset = TRUE) + expect_null(comment(p2@ext_diffusion)) + expect_equal(p2@ext_diffusion, params@ext_diffusion * 0, ignore_attr = TRUE) +}) + +test_that("Comment works on ext_diffusion", { + params <- NS_params + ext_diffusion <- params@ext_diffusion + # comment is stored + comment(ext_diffusion) <- "test" + params <- setExtDiffusion(params, ext_diffusion = ext_diffusion) + expect_identical(comment(params@ext_diffusion), "test") + + # if no comment, previous comment is kept + comment(ext_diffusion) <- NULL + params <- setExtDiffusion(params, ext_diffusion = ext_diffusion) + expect_identical(comment(params@ext_diffusion), "test") +}) + +test_that("ext_diffusion() returns ArraySpeciesBySize", { + expect_true(is.ArraySpeciesBySize(ext_diffusion(NS_params))) + expect_equal(ext_diffusion(NS_params), NS_params@ext_diffusion, + ignore_attr = TRUE) +}) + +test_that("setExtDiffusion validates dimensions", { + expect_error(setExtDiffusion(NS_params, array(0, dim = c(1, 1)))) +}) + +test_that("Can get and set ext_diffusion slot", { + params <- NS_params + ed <- ext_diffusion(params) + new <- ed + 1 + comment(new) <- "test" + ext_diffusion(params) <- new + expect_equal(ext_diffusion(params), new, ignore_attr = TRUE) + expect_identical(comment(params@ext_diffusion), "test") +}) diff --git a/tests/testthat/test-steady.R b/tests/testthat/test-steady.R index eb40e1675..a129d4517 100644 --- a/tests/testthat/test-steady.R +++ b/tests/testthat/test-steady.R @@ -167,6 +167,16 @@ test_that("valid_species_arg works", { NS_params@species_params$species)) }) +test_that("projectToSteady() converges with use_predation_diffusion", { + params_d <- NS_params + params_d@use_predation_diffusion <- TRUE + initialN(params_d)[1, ] <- initialN(params_d)[1, ] * 3 + expect_message( + projectToSteady(params_d, tol = 10), + "Convergence was achieved" + ) +}) + test_that("valid_gears_arg works", { all_gears <- unique(NS_params@gear_params$gear) expect_identical(valid_gears_arg(NS_params), all_gears) diff --git a/tests/testthat/test-steadySingleSpecies.R b/tests/testthat/test-steadySingleSpecies.R index 982554c10..182aafb1d 100644 --- a/tests/testthat/test-steadySingleSpecies.R +++ b/tests/testthat/test-steadySingleSpecies.R @@ -45,7 +45,7 @@ test_that("steadySingleSpecies produces steady state with diffusion", { species <- params@species_params$species[1] n <- params@species_params[species, "n"] d <- 0.1 * params@w^(n + 1) - diffusion(params)[species, ] <- d + ext_diffusion(params)[species, ] <- d # Increase minimum size to test boundary condition params@w_min_idx[species] <- 10 diff --git a/tests/testthat/test-transport.R b/tests/testthat/test-transport.R index 79a650780..5038df5d8 100644 --- a/tests/testthat/test-transport.R +++ b/tests/testthat/test-transport.R @@ -15,7 +15,8 @@ test_that("get_transport_coefs works correctly", { getEGrowth <- mizer:::getEGrowth getMort <- mizer:::getMort - coefs <- get_transport_coefs(params, n, getEGrowth(params), getMort(params), dt, recruitment_flux) + d <- params@ext_diffusion + coefs <- get_transport_coefs(params, n, getEGrowth(params), getMort(params), dt, recruitment_flux, d) # Check dimensions expect_equal(dim(coefs$a), dim(n)) @@ -58,8 +59,8 @@ test_that("get_transport_coefs works correctly", { dw <- params@dw # Species 2 - expected_b_2 <- 1 + dt * mu[2, j_start_2] + - (dt / dw[j_start_2]) * (g[2, j_start_2] + 0.5 * params@diffusion[2, j_start_2] / dw[j_start_2]) + expected_b_2 <- 1 + dt * mu[2, j_start_2] + + (dt / dw[j_start_2]) * (g[2, j_start_2] + 0.5 * d[2, j_start_2] / dw[j_start_2]) expect_equal(coefs$b[2, j_start_2], expected_b_2, ignore_attr = TRUE) }) diff --git a/vignettes/_quarto.yml b/vignettes/_quarto.yml new file mode 100644 index 000000000..1f499b51c --- /dev/null +++ b/vignettes/_quarto.yml @@ -0,0 +1,2 @@ +project: + render: ['*.qmd'] diff --git a/vignettes/analytic_test.Rmd b/vignettes/analytic_test.Rmd index c4e40a5e5..b424dabcb 100644 --- a/vignettes/analytic_test.Rmd +++ b/vignettes/analytic_test.Rmd @@ -144,7 +144,7 @@ params <- setRateFunction(params, "Mort", "start_mort") params <- setRateFunction(params, "RDD", "constant_rdd") # Set diffusion -params@diffusion[1, ] <- K * w(params)^(p+1) +ext_diffusion(params)[1, ] <- K * w(params)^(p + 1) # We also need to switch off resource dynamics and other things to avoid interference params <- setResource(params, resource_dynamics = "resource_constant") @@ -443,7 +443,7 @@ run_numerical_diffusion_test <- function(no_w) { # NO physical diffusion # We set the diffusion matrix to zero - params@diffusion[] <- 0 + ext_diffusion(params)[] <- 0 params <- setResource(params, resource_dynamics = "resource_constant") initialNResource(params) <- 0 diff --git a/vignettes/cohort_dynamics_and_diffusion.Rmd b/vignettes/cohort_dynamics_and_diffusion.Rmd index 7defb2b73..0e3f8fd15 100644 --- a/vignettes/cohort_dynamics_and_diffusion.Rmd +++ b/vignettes/cohort_dynamics_and_diffusion.Rmd @@ -151,7 +151,7 @@ moves to the right (towards larger sizes) as the fish grow. # Adding diffusion Now let's add diffusion to the model. Diffusion is set as an array with -dimensions species × size via `setDiffusion()`. We'll set a constant +dimensions species × size via `setExtDiffusion()`. We'll set a constant diffusion rate across all sizes. We define a helper function that runs the simulation for a given @@ -160,9 +160,9 @@ diffusion coefficient: ```{r} run_with_diffusion <- function(params_base, diff_coeff, diff_exp, t_max = 5) { p <- params_base - d <- p@diffusion + d <- p@ext_diffusion d[] <- diff_coeff * w ^ diff_exp - p <- setDiffusion(p, diffusion = d) + p <- setExtDiffusion(p, ext_diffusion = d) initialN(p)[] <- 0 sim <- project(p, t_max = t_max, dt = 0.05, t_save = 0.1, progress_bar = FALSE) @@ -308,6 +308,6 @@ This vignette demonstrated: 2. How to implement pulsed annual reproduction using a custom RDD function. 3. How cohorts of fish grow through the size spectrum over time. -4. How diffusion, set via `setDiffusion()`, controls the spreading of +4. How diffusion, set via `setExtDiffusion()`, controls the spreading of cohorts — representing individual variability in growth rates. diff --git a/vignettes/diffusion_references.bib b/vignettes/diffusion_references.bib new file mode 100644 index 000000000..c310969e2 --- /dev/null +++ b/vignettes/diffusion_references.bib @@ -0,0 +1,61 @@ + +@article{datta2016, + title = {The effects of seasonal processes on size spectrum dynamics}, + author = {Datta, Samik and Blanchard, Julia L.}, + year = {2016}, + month = {04}, + date = {2016-04}, + journal = {Canadian Journal of Fisheries and Aquatic Sciences}, + pages = {598--610}, + volume = {73}, + number = {4}, + doi = {10.1139/cjfas-2015-0468}, + url = {https://cdnsciencepub.com/doi/full/10.1139/cjfas-2015-0468}, + note = {26 citations (Crossref) [2024-11-14] +00000 +Publisher: NRC Research Press} +} + +@article{datta2010, + title = {A Jump-Growth Model for Predator-Prey Dynamics: Derivation and Application to Marine Ecosystems}, + author = {Datta, Samik and Delius, Gustav W. and Law, Richard}, + year = {2010}, + month = {01}, + date = {2010-01}, + journal = {Bulletin of Mathematical Biology}, + pages = {1361--1382}, + volume = {72}, + number = {6}, + doi = {10.1007/s11538-009-9496-5}, + url = {http://arxiv.org/abs/0812.4968}, + note = {00036} +} + +@article{capitan2010, + title = {Scale-invariant model of marine population dynamics}, + author = {Capitan, Jose A and Delius, Gustav W.}, + year = {2010}, + month = {06}, + date = {2010-06-01}, + journal = {Physical Review E}, + pages = {061901}, + volume = {81}, + number = {6}, + doi = {10.1103/PhysRevE.81.061901}, + url = {http://arxiv.org/abs/1003.4674}, + note = {00018} +} + +@article{datta2010a, + title = {A stability analysis of the power-law steady state of marine size spectra}, + author = {Datta, Samik and Delius, Gustav W. and Law, Richard and Plank, Michael J.}, + year = {2010}, + month = {12}, + date = {2010-12-14}, + journal = {Journal of Mathematical Biology}, + pages = {779--799}, + volume = {63}, + doi = {10.1007/s00285-010-0387-z}, + url = {http://www.springerlink.com/content/70g2073851314707/}, + note = {00032} +} diff --git a/vignettes/mizer.Rmd b/vignettes/mizer.Rmd index 30505e98e..628cd05ed 100644 --- a/vignettes/mizer.Rmd +++ b/vignettes/mizer.Rmd @@ -79,7 +79,7 @@ installing a development version of mizer, click on the triangle below to reveal further details:
-Mizer is compatible with R versions 3.1 and later. You can install R on your +Mizer is compatible with R versions 3.5 and later. You can install R on your computer by following the instructions at for your particular platform. diff --git a/vignettes/predation_diffusion.qmd b/vignettes/predation_diffusion.qmd new file mode 100644 index 000000000..2eeb76a74 --- /dev/null +++ b/vignettes/predation_diffusion.qmd @@ -0,0 +1,186 @@ +--- +title: "Incorporating predation diffusion" +bibliography: diffusion_references.bib +toc: true +crossref: + eq-prefix: Eq. +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r} +#| message: false +#| warning: false +library(mizerEcopath) +library(ggplot2) +library(tidyr) +library(dplyr) +``` + +# The jump-growth equation {#sec-jump-growth} + +As was observed in [@datta2010], the variability in prey size leads to a diffusion term in the PDE for the abundance density $N(w,t)$. This term was neglected in mizer in the past but can now be included. It becomes important when matching mizer abundances to observed abundances, because it will account for fish that are larger than average. The PDE including the diffusion term is +$$ +\frac{\partial N}{\partial t} = \frac12 \frac{\partial^2}{\partial w^2}(d N) - \frac{\partial}{\partial w}(g N)-\mu N +$$ {#eq-PDE} +where $g$ is the growth rate, $\mu$ is the death rate and $d$ is a new diffusion rate. + +The diffusion rate has an expression that is similar to that of the rate $A(w)$ at which prey biomass is assimilated and thus available for metabolism, growth and reproduction. Recall that this rate is given by +$$ +\begin{split} +A(w)=&(1-f(w))\,\gamma(w)\int N_c(w_p)\,\phi(w/w_p)\,\alpha\, w_p\,dw_p +\end{split} +$$ {#eq-gw} +where $N_c(w)$ is the abundance density of prey, $\phi(w/w_p)$ is the predation kernel, $\gamma(w)$ is the search volume, $f(w)$ is the feeding level and $\alpha$ is the assimilation efficiency. + +The factor $\alpha\,w_p$ in the integral in @eq-gw is the increase in assimilated biomass resulting from the ingestion of an individual prey of weight $w_p$. In the expression for the diffusion rate $D(w)$ this factor is squared: +$$ +D(w) = (1-f(w))\,\gamma(w) \int N_c(w_p)\phi(w/w_p)(\alpha\,w_p)^2\,dw_p. +$$ {#eq-dw} + +This diffusion rate captures only the randomness in the size of encountered prey. There will be many other sources of stochasticity that we do not model explicitly in mizer. They should be added as external diffusion with `ext_diffusion()`. + +The increase $\alpha\,w_p$ in predator mass is typically only a small proportion of the predator mass because the preferred prey are typically much smaller than the predator and also $\alpha$ is smaller than $1$. Because this factor is squared in the expression for $D(w)$ it might be expected that the term in the PDE involving $D(w)$ can be safely neglected. However it is worth testing this intuition. We will do this now by first determining the diffusion rate in a model with allometric encounter and mortality rates. We will then use that to determine its effect on the slope of the juvenile spectrum. Finally we will look at the numerical solution for the steady state. + +# Example calculation of diffusion rate {#sec-diffusion-example} + +Let us assume that the prey abundance is given by a power law: $N_c(w)=N_0w^{-\lambda}$ and that the predation kernel is +$$ +\phi(w/w_p) = \exp\left(-\frac{\log(w/w_p/\beta)^2}{2\sigma^2}\right). +$$ {#eq-phi} + +and hence the integral in the expression @eq-gw for the assimilation rate becomes +$$ +\begin{split} +I_A&:=\int N_c(w_p)\phi(w/w_p)\alpha\,w_p\,dw_p\\ +&=\alpha\int w_p^{1-\lambda}\exp\left(-\frac{\log(w/w_p/\beta)^2}{2\sigma^2}\right)dw_p. +\end{split} +$$ {#eq-Ig} +This integral can be evaluated most easily by changing integration variable to $x=\log(w_p/w_0)$for an arbitrary reference weight $w_0$ and then recognising the resulting integral +$$ +I_A=\alpha\, w_0^{2-\lambda}\int e^{(2-\lambda)x}\exp\left(-\frac{(x-\log(w/w_0)+\log(\beta))^2}{2\sigma^2}\right)dx +$$ {#eq-Ig2} +as a Gaussian integral. Using the general result that +$$ +\int e^{ax}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)dx = \sqrt{\frac{2\pi}{\sigma^2}}\exp\left(a\mu+\frac{a^2\sigma^2}{2}\right) +$$ {#eq-gaussian} +with $a = 2-\lambda$ and $\mu = \log(w/w_0)-\log(\beta)$ we find that +$$ +\begin{split} +I_A&= +\alpha\, w_0^{2-\lambda}\sqrt{\frac{2\pi}{\sigma^2}}\exp\left((2-\lambda)(\log(w/w_0)-\log(\beta))+\frac{(2-\lambda)^2\sigma^2}{2}\right)\\ +&=\alpha w^{2-\lambda}\sqrt{\frac{2\pi}{\sigma^2}}\beta^{\lambda - 2}\exp\left(\frac{(2-\lambda)^2\sigma^2}{2}\right) +\end{split} +$$ {#eq-Ig3} + +Assuming an allometric search volume $\gamma(w)=\gamma w^q$ with an exponent of $q = n - 2 + \lambda$ we obtain +$$ +A(w) = (1-f(w))\,\gamma \alpha \sqrt{\frac{2\pi}{\sigma^2}}\beta^{\lambda - 2}\exp\left(\frac{(2-\lambda)^2\sigma^2}{2}\right) w^{n}. +$$ {#eq-gw2} + +We can evaluate the diffusion rate $D(w)$ in the same manner. The extra factor of $w_p$ changes $a$ from $2-\lambda$ to $3-\lambda$ and we obtain +$$ +D(w) = (1-f(w))\gamma \alpha^2 \sqrt{\frac{2\pi}{\sigma^2}}\beta^{\lambda - 3}\exp\left(\frac{(3-\lambda)^2\sigma^2}{2}\right) w^{n+1}. +$$ {#eq-dw2} +Comparing this to the expression @eq-gw2 for $A(w)$ we find that +$$ +D(w) = A(w)w\frac{\alpha}{\beta}\exp\left(\frac{(5-2\lambda)\sigma^2}{2}\right). +$$ {#eq-dw3} + +To get a feel for the typical magnitude of the factor let's consider concrete values +$$ +f=0.6, \qquad\alpha = 0.8, \qquad\beta = 10, \qquad\sigma = 2, \qquad\lambda = 2. +$$ {#eq-params} +```{r} +#| include: false +feeding_coeff <- 0.6 +alpha <- 0.8 +beta <- 10 +sigma <- 2 +lambda <- 2 +factor <- alpha / beta * exp((5 - 2 * lambda) * sigma^2 / 2) +``` + +Then +$$ +D(w) \approx `r round(factor, digits = 2)` A(w)w. +$$ {#eq-dw4} + +```{r} +#| include: false +sigma <- 1 +beta <- 100 +factor_1 <- alpha / beta * exp((5 - 2 * lambda) * sigma^2 / 2) +``` + +However we see that the value of $D(w)/A(w)$ is strongly influenced by the location and width of the feeding kernel. It decreases with increasing $\beta$ and decreasing $\sigma$ If we choose $\beta = 100$ and $\sigma = 1$ then $D(w)\approx `r round(factor_1, digits = 2)` A(w)w$. This relative decrease in the diffusion rate with increasing $\beta$ and decreasing $\sigma$ may explains the result from [@datta2010a] about how the stability of the system depends on these parameters. + +# Effect on juvenile slope {#sec-diffusion-juvenile} + +In this section we will calculate the effect of the diffusion on the slope of the juvenile spectrum in the steady state. You can skip this section if you are not interested in the details of the calculation. The outcome is that the change in the juvenile slope is small. However, the change in the size spectrum of the adults is much larger, as we will see in @sec-diffusion-numerical. + +Without diffusion we find the juvenile spectrum by solving the steady-state equation +$$ +\frac{\partial}{\partial w}(g(w) N(w)) =-\mu(w) N(w). +$$ {#eq-ss} + +This has the solution +$$ +N(w) = \frac{g(w_0)}{g(w)}N(w_0)\exp\left(-\int_{w_0}^w\frac{\mu(w')}{g(w')}dw'\right). +$$ {#eq-Nw} + +With allometric growth and death rates $g(w)=g_0w^n$ and $\mu(w)=\mu_0w^{n-1}$ this gives +$$ +N(w)=\left(\frac{w}{w_0}\right)^{-n} N(w_0)\exp\left(-\frac{\mu_0}{g_0}\int_{w_0}^w\frac{1}{w'}dw'\right)=N(w_0)\left(\frac{w}{w_0}\right)^{-\mu_0/g_0-n}. +$$ {#eq-Nwss} + +Thus the juvenile steady state abundance density is given by a power law with exponent $-\mu_0/g_0-n$. + +In the presence of diffusion the steady state equation becomes the second-order ODE +$$ +\frac12 \frac{\partial^2}{\partial w^2}(D N) - \frac{\partial}{\partial w}(g N)-\mu N=0. +$$ {#eq-diffss} + +We have seen in the previous section that with $g(w)=g_0w^n$ the diffusion rate is also a power law with one extra factor of $w$: $D(w)=D_0w^{n+1}$. This makes the equation @eq-diffss scale invariant and hence we again expect the solution to be a power law, So we make the Ansatz $N(w)=N(w_0)(w/w_0)^a$ with the exponent $a$ to be determined. Substituting this Ansatz into @eq-diffss gives +$$ +\frac12 D_0N_0(n+1+a)(n+a)w^{n+a-1} -g_0N_0(n+a)w^{n+a-1}-\mu_0 N_0w^{n-1+a} = 0 +$$ {#eq-diffssansatz} +which requires that +$$ +\frac12 D_0 (n+a)^2+\left(\frac12 D_0 -g_0\right)(n+a)-\mu_0=0. +$$ {#eq-diffssansatz2} +This is a quadratic equation for $n+a$: +$$ +\frac12 D_0 (n+a)^2+\left(\frac12 D_0 -g_0\right)(n+a)-\mu_0=0. +$$ {#eq-diffssansatz3} +This has two solutions +$$ +(n+a_\pm) = \frac1{D_0}\left(g_0-D_0/2\pm\sqrt{(g_0-D_0/2)^2+2D_0\mu_0}\right) +$$ {#eq-a} +where $a_+$ is the solution with the + sign and $a_-$ is the solution with the - sign. + +We are only interested in the solution that goes to $a=-\mu_0/g_0-n$ when $D_0\to 0$. This means we are interested in the solution with the - sign. To check that indeed the solution $a_{-}$ satisfies $\lim_{a_-\to 0}= -\mu_0/g_0-n$ it is helpful to expand the square root term around $D_0=0$: +$$ +\sqrt{(g_0-D_0/2)^2+2D_0\mu_0}=g_0+\frac{-g_0+2\mu_0}{2g_0}D_0+\frac{g_0^2-(-g_0+2\mu_0)^2}{8g_0^3}D_0^2+\dots +$$ {#eq-sqrt} +So we find +$$ +\begin{split} +a &= -n+\frac1{D_0}\left(g_0-D_0/2-\sqrt{(g_0-D_0/2)^2+2D_0\mu_0}\right)\\ +&\approx -n-\frac{\mu_0}{g_0} - \frac18\left(1-\left(-1+2\frac{\mu_0}{g_0}\right)^2\right) D_0+\dots\\ +&=-n-\frac{\mu_0}{g_0} +\frac12\left(\frac{\mu_0}{g_0}\left(1-\frac{\mu_0}{g_0}\right)\right)\frac{D_0}{g_0}+\dots +\end{split} +$$ {#eq-a-2} + +We see that when $\mu_0=g_0$ then the first correction term to the juvenile slope vanishes. The largest increase in slope (i.e., the least negative slope) is achieved when $\mu_0/g_0=1/2$. In that case the slope without diffusion is $-1.25$ (assuming $n=0.75$). The correction term then is $D_0/(8g_0)$. We had already seen in the previous section that $D_0/g_0$ is typically very small, so the change in slope is also very small. + +Using the example value of $D_0/g_0\approx`r round(factor, digits=2)`$ and $n=0.75$ we get a slope correction of approximately $`r round(factor/8, digits=3)`$ from $-1.25$ to $`r round(factor/8-1.25, digits=3)`$. + +Using the value $D_0/g_0\approx`r round(factor_1, digits=2)`$ we get a smaller slope correction of approximately $`r round(factor_1/8, digits=3)`$ from $-1.25$ to $`r round(factor_1/8-1.25, digits=3)`$. + +When $\mu_0/g_0>1$ then the diffusion correction is negative, i.e., the juvenile slope becomes more negative due to diffusion, meaning fewer large fish.