diff --git a/NAMESPACE b/NAMESPACE index b783a42e3..03914e28b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -99,6 +99,8 @@ S3method(getSSB,MizerParams) S3method(getSSB,MizerSim) S3method(getSearchVolume,MizerParams) S3method(getSelectivity,MizerParams) +S3method(getTrophicLevel,MizerParams) +S3method(getTrophicLevelBySpecies,MizerParams) S3method(getYield,MizerParams) S3method(getYield,MizerSim) S3method(getYieldGear,MizerParams) @@ -315,6 +317,8 @@ export(getSSB) export(getSearchVolume) export(getSelectivity) export(getTimes) +export(getTrophicLevel) +export(getTrophicLevelBySpecies) export(getYield) export(getYieldGear) export(getZ) diff --git a/NEWS.md b/NEWS.md index ba7ba72d5..4b8e04d66 100644 --- a/NEWS.md +++ b/NEWS.md @@ -24,6 +24,13 @@ `ylim` arguments for controlling the y-axis scale and axis limits. The `plot()` method for `ArraySpeciesByTime` gains `log` and `ylim` arguments. +- New `getTrophicLevel()` function returns a matrix (species × size) with the + trophic level of individuals at each size, accounting for ontogenetic diet + shifts by integrating the consumption-weighted average prey trophic level + over the individual's growth trajectory. New `getTrophicLevelBySpecies()` + returns the consumption-rate-weighted mean trophic level per species. + Both functions accept `MizerParams` and `MizerSim` objects. Closes #307. + - The `plot()` and `summary()` methods for `MizerParams`, `MizerSim`, and `ArraySpeciesBySize` objects are now registered as S3 methods rather than S4 methods. This means `plot()` and `summary()` remain plain S3 generics when mizer is diff --git a/R/summary_methods.R b/R/summary_methods.R index 0c46c1b39..07195adfa 100644 --- a/R/summary_methods.R +++ b/R/summary_methods.R @@ -20,6 +20,8 @@ #' \tabular{lll}{ #' Function \tab Returns \tab Description \cr #' [getDiet()] \tab Three dimensional array (predator x size x prey) \tab Diet of predator at size, resolved by prey species \cr +#' [getTrophicLevel()] \tab `ArraySpeciesBySize` (species x size) \tab Trophic level of individuals at size, accounting for ontogenetic diet shifts \cr +#' [getTrophicLevelBySpecies()] \tab Named vector (species) \tab Consumption-rate-weighted mean trophic level of each species \cr #' [getSSB()] \tab Two dimensional array (time x species) \tab Total Spawning Stock Biomass (SSB) of each species through time where SSB is calculated as the sum of weight of all mature individuals. \cr #' [getBiomass()] \tab Two dimensional array (time x species) \tab Total biomass of each species through time. \cr #' [getN()] \tab Two dimensional array (time x species) \tab Total abundance of each species through time. \cr @@ -176,6 +178,203 @@ getDiet.MizerParams <- function(params, } +#' Get trophic level of individuals at size +#' +#' `r lifecycle::badge("experimental")` +#' Calculates the trophic level of individuals of each species at each size, +#' assuming the system is in a steady state. The trophic level of an individual +#' is defined as 1 more than the consumption-rate-weighted average trophic level +#' of all the prey it has consumed during its lifetime up to the current size. +#' The trophic level of the primary resource is set to 0. +#' +#' +#' @details +#' In the traditional non-size-resolved approach, all individuals of a species +#' have the same diet composition \eqn{D_{ij}}, defined as the proportion of +#' total biomass intake of species \eqn{i} that comes from species \eqn{j}. +#' The trophic levels then satisfy +#' \deqn{T_i = 1 + \sum_j D_{ij}\,T_j,} +#' which is solved as a linear system \eqn{(I - D)\,\mathbf{T} = \mathbf{1}}. +#' +#' In mizer, diet composition changes as an individual grows, so we must +#' integrate over the individual's lifetime. Assuming a steady state so that +#' the growth rate \eqn{g_i(w)} and prey densities depend only on size and not +#' on time, we can replace the integral over time since birth by an integral +#' over weight using \eqn{dt = dw / g_i(w)}. The trophic level +#' \eqn{T_i(w)} of an individual of species \eqn{i} at weight \eqn{w} is +#' then +#' \deqn{ +#' T_i(w) = 1 + \frac{ +#' \int_{w_0}^{w} \frac{1}{g_i(w')} \sum_j \int r_{ij}(w', w_p)\, T_j(w_p)\, dw_p\, dw' +#' }{ +#' \int_{w_0}^{w} \frac{1}{g_i(w')} \sum_j \int r_{ij}(w', w_p)\, dw_p\, dw' +#' }, +#' } +#' where \eqn{w_0} is the egg size and \eqn{r_{ij}(w, w_p)} is the rate at +#' which a predator of species \eqn{i} at weight \eqn{w} consumes biomass from +#' prey species \eqn{j} at weight \eqn{w_p}: +#' \deqn{ +#' r_{ij}(w, w_p) = \theta_{ij}\,\gamma_i(w)\,(1 - f_i(w))\,\phi_i(w/w_p)\, +#' N_j(w_p)\,w_p. +#' } +#' The sum over \eqn{j} runs over all species. The resource is excluded from +#' the numerator because its trophic level is 0, but is included in the +#' denominator (which equals the total biomass consumed over the predator's +#' lifetime from egg size to current weight \eqn{w}). +#' +#' This equation can be viewed as a linear system +#' \eqn{(I - D)\,\mathbf{T} = \mathbf{1}} in which the entries of +#' \eqn{\mathbf{T}} are indexed by \eqn{(i, w)} and the matrix \eqn{D} encodes +#' the lifetime-integrated diet composition. The system is solved iteratively +#' from small to large sizes, exploiting the fact that prey are typically much +#' smaller than the predator (large predator-to-prey mass ratio), so that the +#' trophic levels of all relevant prey sizes are already known when computing +#' \eqn{T_i(w)}. +#' +#' @inheritParams getDiet +#' +#' @return An `ArraySpeciesBySize` object (species x size) with the trophic +#' level of individuals at each size. Entries below the egg size of each +#' species are \code{NA}. +#' +#' @export +#' @family summary functions +#' @concept summary_function +#' @seealso [getTrophicLevelBySpecies()] +#' @examples +#' tl <- getTrophicLevel(NS_params) +#' plot(tl) +getTrophicLevel <- function(params, + n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + ...) { + UseMethod("getTrophicLevel") +} + +#' @export +getTrophicLevel.MizerParams <- function(params, + n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + ...) { + params <- validParams(params) + no_sp <- nrow(params@species_params) + no_w <- length(params@w) + no_w_full <- length(params@w_full) + idx_sp <- (no_w_full - no_w + 1):no_w_full + + # Pre-compute rates + encounter <- getEncounter(params, n, n_pp, n_other) # no_sp x no_w + feeding_level <- getFeedingLevel(params, n, n_pp, n_other) # no_sp x no_w + growth <- getEGrowth(params, n, n_pp, n_other) # no_sp x no_w + # Total consumption = (1 - f) * E, used for denominator accumulator + consumption <- (1 - feeding_level) * encounter # no_sp x no_w + + # Full predation kernel array (no_sp x no_w x no_w_full). + # getPredKernel() computes it from species parameters if not explicitly stored. + pred_kernel <- getPredKernel(params) + + # prey_mass_tl[j, p] = N_j(w_p) * T_j(w_p) * w_p * dw_p + # Initialised with T_j = 1 everywhere; updated as trophic levels are computed + prey_mass_tl <- sweep(n, 2, params@w * params@dw, "*") # no_sp x no_w + + # Cumulative numerator A_i and denominator B_i (integrals weighted by 1/g dw) + cumA <- numeric(no_sp) + cumB <- numeric(no_sp) + + # Output matrix: NA below egg size of each species + tl <- matrix(NA_real_, nrow = no_sp, ncol = no_w, + dimnames = dimnames(params@initial_n)) + + # Iterate from smallest to largest size, building up trophic levels + for (k in seq_len(no_w)) { + # Trophic-level-weighted encounter for all predator species at size w[k]: + # E_tl[i] = gamma_i(w_k) * sum_j theta_ij * sum_p kernel[i,k,p] * N_tl[j,p] * w_p * dw_p + pred_kernel_k <- matrix(pred_kernel[, k, idx_sp], nrow = no_sp) + ae_k <- pred_kernel_k %*% t(prey_mass_tl) # no_sp x no_sp + E_tl <- params@search_vol[, k] * rowSums(params@interaction * ae_k) + + # Update trophic level for each species active at this size + for (i in seq_len(no_sp)) { + if (k < params@w_min_idx[i]) next + g_ik <- growth[i, k] + if (!is.finite(g_ik) || g_ik <= 0) { + # No growth: carry forward previous trophic level + tl[i, k] <- if (k > params@w_min_idx[i]) tl[i, k - 1L] else 1 + next + } + weight <- params@dw[k] / g_ik + cumA[i] <- cumA[i] + (1 - feeding_level[i, k]) * E_tl[i] * weight + cumB[i] <- cumB[i] + consumption[i, k] * weight + tl[i, k] <- if (cumB[i] > 0) 1 + cumA[i] / cumB[i] else 1 + } + + # Update prey_mass_tl for size k with newly computed trophic levels + active_k <- k >= params@w_min_idx + tl_k <- tl[, k] + tl_k[is.na(tl_k)] <- 1 + prey_mass_tl[active_k, k] <- + n[active_k, k] * tl_k[active_k] * params@w[k] * params@dw[k] + } + + return(ArraySpeciesBySize(tl, value_name = "Trophic level", params = params)) +} + + +#' Get mean trophic level of each species +#' +#' `r lifecycle::badge("experimental")` +#' Calculates the consumption-rate-weighted mean trophic level of each species, +#' defined as +#' \deqn{ +#' T_i = \frac{\int r_i(w)\,N_i(w)\,T_i(w)\,dw} +#' {\int r_i(w)\,N_i(w)\,dw}, +#' } +#' where \eqn{r_i(w) = (1 - f_i(w))\,E_i(w)} is the consumption rate of an +#' individual of species \eqn{i} at weight \eqn{w}, \eqn{N_i(w)} is the +#' abundance density, and \eqn{T_i(w)} is the size-resolved trophic level +#' from [getTrophicLevel()]. +#' +#' @inheritParams getDiet +#' +#' @return A named vector with the mean trophic level for each species. +#' +#' @export +#' @family summary functions +#' @concept summary_function +#' @seealso [getTrophicLevel()] +#' @examples +#' getTrophicLevelBySpecies(NS_params) +getTrophicLevelBySpecies <- function(params, + n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + ...) { + UseMethod("getTrophicLevelBySpecies") +} + +#' @export +getTrophicLevelBySpecies.MizerParams <- function(params, + n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + ...) { + params <- validParams(params) + tl <- getTrophicLevel(params, n = n, n_pp = n_pp, n_other = n_other) + tl[is.na(tl)] <- 0 + encounter <- getEncounter(params, n, n_pp, n_other) + feeding_level <- getFeedingLevel(params, n, n_pp, n_other) + # Consumption rate per individual times abundance density + consumption_n <- (1 - feeding_level) * encounter * n # no_sp x no_w + # Weighted mean trophic level: integral of consumption_n * T * dw / integral of consumption_n * dw + numerator <- (consumption_n * tl) %*% params@dw + denominator <- consumption_n %*% params@dw + tl_by_sp <- numerator[, 1] / denominator[, 1] + tl_by_sp[denominator[, 1] == 0] <- NA_real_ + return(tl_by_sp) +} + #' Calculate the SSB of species #' #' Calculates the spawning stock biomass (SSB) for each species. For a diff --git a/man/getBiomass.Rd b/man/getBiomass.Rd index 71fbf9502..00637f978 100644 --- a/man/getBiomass.Rd +++ b/man/getBiomass.Rd @@ -68,6 +68,8 @@ Other summary functions: \code{\link{getGrowthCurves}()}, \code{\link{getN}()}, \code{\link{getSSB}()}, +\code{\link{getTrophicLevel}()}, +\code{\link{getTrophicLevelBySpecies}()}, \code{\link{getYield}()}, \code{\link{getYieldGear}()} } diff --git a/man/getDiet.Rd b/man/getDiet.Rd index 9f6f5e00f..3653ce2fc 100644 --- a/man/getDiet.Rd +++ b/man/getDiet.Rd @@ -72,6 +72,8 @@ Other summary functions: \code{\link{getGrowthCurves}()}, \code{\link{getN}()}, \code{\link{getSSB}()}, +\code{\link{getTrophicLevel}()}, +\code{\link{getTrophicLevelBySpecies}()}, \code{\link{getYield}()}, \code{\link{getYieldGear}()} } diff --git a/man/getGrowthCurves.Rd b/man/getGrowthCurves.Rd index b575b8aaf..ccaf5008c 100644 --- a/man/getGrowthCurves.Rd +++ b/man/getGrowthCurves.Rd @@ -44,6 +44,8 @@ Other summary functions: \code{\link{getDiet}()}, \code{\link{getN}()}, \code{\link{getSSB}()}, +\code{\link{getTrophicLevel}()}, +\code{\link{getTrophicLevelBySpecies}()}, \code{\link{getYield}()}, \code{\link{getYieldGear}()} } diff --git a/man/getN.Rd b/man/getN.Rd index 97cd2ab67..9c76fa744 100644 --- a/man/getN.Rd +++ b/man/getN.Rd @@ -48,6 +48,8 @@ Other summary functions: \code{\link{getDiet}()}, \code{\link{getGrowthCurves}()}, \code{\link{getSSB}()}, +\code{\link{getTrophicLevel}()}, +\code{\link{getTrophicLevelBySpecies}()}, \code{\link{getYield}()}, \code{\link{getYieldGear}()} } diff --git a/man/getSSB.Rd b/man/getSSB.Rd index 56c1cedae..cae5acbad 100644 --- a/man/getSSB.Rd +++ b/man/getSSB.Rd @@ -31,6 +31,8 @@ Other summary functions: \code{\link{getDiet}()}, \code{\link{getGrowthCurves}()}, \code{\link{getN}()}, +\code{\link{getTrophicLevel}()}, +\code{\link{getTrophicLevelBySpecies}()}, \code{\link{getYield}()}, \code{\link{getYieldGear}()} } diff --git a/man/getTrophicLevel.Rd b/man/getTrophicLevel.Rd new file mode 100644 index 000000000..607d11e18 --- /dev/null +++ b/man/getTrophicLevel.Rd @@ -0,0 +1,99 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/summary_methods.R +\name{getTrophicLevel} +\alias{getTrophicLevel} +\title{Get trophic level of individuals at size} +\usage{ +getTrophicLevel( + params, + n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + ... +) +} +\arguments{ +\item{params}{A \linkS4class{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 of the +ecosystem} +} +\value{ +An \code{ArraySpeciesBySize} object (species x size) with the trophic +level of individuals at each size. Entries below the egg size of each +species are \code{NA}. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +Calculates the trophic level of individuals of each species at each size, +assuming the system is in a steady state. The trophic level of an individual +is defined as 1 more than the consumption-rate-weighted average trophic level +of all the prey it has consumed during its lifetime up to the current size. +The trophic level of the primary resource is set to 0. +} +\details{ +In the traditional non-size-resolved approach, all individuals of a species +have the same diet composition \eqn{D_{ij}}, defined as the proportion of +total biomass intake of species \eqn{i} that comes from species \eqn{j}. +The trophic levels then satisfy +\deqn{T_i = 1 + \sum_j D_{ij}\,T_j,} +which is solved as a linear system \eqn{(I - D)\,\mathbf{T} = \mathbf{1}}. + +In mizer, diet composition changes as an individual grows, so we must +integrate over the individual's lifetime. Assuming a steady state so that +the growth rate \eqn{g_i(w)} and prey densities depend only on size and not +on time, we can replace the integral over time since birth by an integral +over weight using \eqn{dt = dw / g_i(w)}. The trophic level +\eqn{T_i(w)} of an individual of species \eqn{i} at weight \eqn{w} is +then +\deqn{ +T_i(w) = 1 + \frac{ + \int_{w_0}^{w} \frac{1}{g_i(w')} \sum_j \int r_{ij}(w', w_p)\, T_j(w_p)\, dw_p\, dw' +}{ + \int_{w_0}^{w} \frac{1}{g_i(w')} \sum_j \int r_{ij}(w', w_p)\, dw_p\, dw' +}, +} +where \eqn{w_0} is the egg size and \eqn{r_{ij}(w, w_p)} is the rate at +which a predator of species \eqn{i} at weight \eqn{w} consumes biomass from +prey species \eqn{j} at weight \eqn{w_p}: +\deqn{ +r_{ij}(w, w_p) = \theta_{ij}\,\gamma_i(w)\,(1 - f_i(w))\,\phi_i(w/w_p)\, + N_j(w_p)\,w_p. +} +The sum over \eqn{j} runs over all species. The resource is excluded from +the numerator because its trophic level is 0, but is included in the +denominator (which equals the total biomass consumed over the predator's +lifetime from egg size to current weight \eqn{w}). + +This equation can be viewed as a linear system +\eqn{(I - D)\,\mathbf{T} = \mathbf{1}} in which the entries of +\eqn{\mathbf{T}} are indexed by \eqn{(i, w)} and the matrix \eqn{D} encodes +the lifetime-integrated diet composition. The system is solved iteratively +from small to large sizes, exploiting the fact that prey are typically much +smaller than the predator (large predator-to-prey mass ratio), so that the +trophic levels of all relevant prey sizes are already known when computing +\eqn{T_i(w)}. +} +\examples{ +tl <- getTrophicLevel(NS_params) +plot(tl) +} +\seealso{ +\code{\link[=getTrophicLevelBySpecies]{getTrophicLevelBySpecies()}} + +Other summary functions: +\code{\link{getBiomass}()}, +\code{\link{getDiet}()}, +\code{\link{getGrowthCurves}()}, +\code{\link{getN}()}, +\code{\link{getSSB}()}, +\code{\link{getTrophicLevelBySpecies}()}, +\code{\link{getYield}()}, +\code{\link{getYieldGear}()} +} +\concept{summary functions} +\concept{summary_function} diff --git a/man/getTrophicLevelBySpecies.Rd b/man/getTrophicLevelBySpecies.Rd new file mode 100644 index 000000000..e33d9d5ab --- /dev/null +++ b/man/getTrophicLevelBySpecies.Rd @@ -0,0 +1,58 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/summary_methods.R +\name{getTrophicLevelBySpecies} +\alias{getTrophicLevelBySpecies} +\title{Get mean trophic level of each species} +\usage{ +getTrophicLevelBySpecies( + params, + n = initialN(params), + n_pp = initialNResource(params), + n_other = initialNOther(params), + ... +) +} +\arguments{ +\item{params}{A \linkS4class{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 of the +ecosystem} +} +\value{ +A named vector with the mean trophic level for each species. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} +Calculates the consumption-rate-weighted mean trophic level of each species, +defined as +\deqn{ + T_i = \frac{\int r_i(w)\,N_i(w)\,T_i(w)\,dw} + {\int r_i(w)\,N_i(w)\,dw}, +} +where \eqn{r_i(w) = (1 - f_i(w))\,E_i(w)} is the consumption rate of an +individual of species \eqn{i} at weight \eqn{w}, \eqn{N_i(w)} is the +abundance density, and \eqn{T_i(w)} is the size-resolved trophic level +from \code{\link[=getTrophicLevel]{getTrophicLevel()}}. +} +\examples{ +getTrophicLevelBySpecies(NS_params) +} +\seealso{ +\code{\link[=getTrophicLevel]{getTrophicLevel()}} + +Other summary functions: +\code{\link{getBiomass}()}, +\code{\link{getDiet}()}, +\code{\link{getGrowthCurves}()}, +\code{\link{getN}()}, +\code{\link{getSSB}()}, +\code{\link{getTrophicLevel}()}, +\code{\link{getYield}()}, +\code{\link{getYieldGear}()} +} +\concept{summary functions} +\concept{summary_function} diff --git a/man/getYield.Rd b/man/getYield.Rd index 6382280d4..d0850177d 100644 --- a/man/getYield.Rd +++ b/man/getYield.Rd @@ -62,6 +62,8 @@ Other summary functions: \code{\link{getGrowthCurves}()}, \code{\link{getN}()}, \code{\link{getSSB}()}, +\code{\link{getTrophicLevel}()}, +\code{\link{getTrophicLevelBySpecies}()}, \code{\link{getYieldGear}()} } \concept{summary functions} diff --git a/man/getYieldGear.Rd b/man/getYieldGear.Rd index 85ab439ff..25e333760 100644 --- a/man/getYieldGear.Rd +++ b/man/getYieldGear.Rd @@ -37,6 +37,8 @@ Other summary functions: \code{\link{getGrowthCurves}()}, \code{\link{getN}()}, \code{\link{getSSB}()}, +\code{\link{getTrophicLevel}()}, +\code{\link{getTrophicLevelBySpecies}()}, \code{\link{getYield}()} } \concept{summary functions} diff --git a/man/summary_functions.Rd b/man/summary_functions.Rd index bbb66fea8..a57b31b8d 100644 --- a/man/summary_functions.Rd +++ b/man/summary_functions.Rd @@ -11,6 +11,8 @@ A list of available summary functions is given in the table below. \tabular{lll}{ Function \tab Returns \tab Description \cr \code{\link[=getDiet]{getDiet()}} \tab Three dimensional array (predator x size x prey) \tab Diet of predator at size, resolved by prey species \cr +\code{\link[=getTrophicLevel]{getTrophicLevel()}} \tab \code{ArraySpeciesBySize} (species x size) \tab Trophic level of individuals at size, accounting for ontogenetic diet shifts \cr +\code{\link[=getTrophicLevelBySpecies]{getTrophicLevelBySpecies()}} \tab Named vector (species) \tab Consumption-rate-weighted mean trophic level of each species \cr \code{\link[=getSSB]{getSSB()}} \tab Two dimensional array (time x species) \tab Total Spawning Stock Biomass (SSB) of each species through time where SSB is calculated as the sum of weight of all mature individuals. \cr \code{\link[=getBiomass]{getBiomass()}} \tab Two dimensional array (time x species) \tab Total biomass of each species through time. \cr \code{\link[=getN]{getN()}} \tab Two dimensional array (time x species) \tab Total abundance of each species through time. \cr diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index ea49bb4db..9b3365e7c 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -156,20 +156,11 @@ reference: - title: Analysing results contents: - summary_functions - - getDiet - - getBiomass - - getN - - getSSB - - getYield - - getYieldGear - - getGrowthCurves + - has_concept("summary functions") - title: Calculating indicators contents: - indicator_functions - - getProportionOfLargeFish - - getCommunitySlope - - getMeanWeight - - getMeanMaxWeight + - has_concept("functions for calculating indicators") - title: Plotting results contents: - plotting_functions diff --git a/tests/testthat/test-summary_methods.R b/tests/testthat/test-summary_methods.R index bb5fdb964..f8dd705f8 100644 --- a/tests/testthat/test-summary_methods.R +++ b/tests/testthat/test-summary_methods.R @@ -363,6 +363,69 @@ test_that("getDiet works with additional components", { }) +# getTrophicLevel ---- +test_that("getTrophicLevel returns matrix with correct structure", { + tl <- getTrophicLevel(params, n, n_pp) + expect_true(is.ArraySpeciesBySize(tl)) + expect_true(is.matrix(tl)) + expect_equal(dim(tl), c(no_sp, no_w)) + expect_equal(dimnames(tl), dimnames(params@initial_n)) + # All trophic levels >= 1 (since T = 1 + non-negative) + expect_true(all(tl >= 1, na.rm = TRUE)) + # Trophic levels should be reasonable (< 10) + expect_true(all(tl < 10, na.rm = TRUE)) +}) + +test_that("getTrophicLevel gives same result with explicit pred_kernel", { + tl1 <- getTrophicLevel(params, n, n_pp) + # Force explicit pred_kernel storage + params2 <- setPredKernel(params, pred_kernel = getPredKernel(params)) + tl2 <- getTrophicLevel(params2, n, n_pp) + expect_equal(tl1, tl2, tolerance = 1e-10) +}) + +test_that("getTrophicLevel increases along body size for apex predators", { + tl <- getTrophicLevel(NS_params) + # For Cod (apex predator), trophic level should increase with size + cod_tl <- tl["Cod", ] + cod_tl <- cod_tl[!is.na(cod_tl)] + # Should be non-decreasing overall (allow small numerical fluctuations) + expect_true(cod_tl[length(cod_tl)] >= cod_tl[1]) +}) + +test_that("getTrophicLevel works for MizerSim", { + tl_sim <- getTrophicLevel(sim) + expect_equal(length(dim(tl_sim)), 3) + expect_equal(dim(tl_sim)[2:3], c(no_sp, no_w)) + expect_equal(names(dimnames(tl_sim)), c("time", "sp", "w")) +}) + +# getTrophicLevelBySpecies ---- +test_that("getTrophicLevelBySpecies returns named vector", { + tl_sp <- getTrophicLevelBySpecies(params, n, n_pp) + expect_true(is.numeric(tl_sp)) + expect_equal(length(tl_sp), no_sp) + expect_equal(names(tl_sp), params@species_params$species) + expect_true(all(tl_sp >= 1, na.rm = TRUE)) +}) + +test_that("getTrophicLevelBySpecies is consistent with getTrophicLevel", { + tl <- getTrophicLevel(params, n, n_pp) + tl_sp <- getTrophicLevelBySpecies(params, n, n_pp) + # Species-level trophic level should be between min and max size-resolved tl + for (i in seq_len(no_sp)) { + tl_range <- range(tl[i, ], na.rm = TRUE) + expect_gte(tl_sp[i], tl_range[1] - 1e-10) + expect_lte(tl_sp[i], tl_range[2] + 1e-10) + } +}) + +test_that("getTrophicLevelBySpecies works for MizerSim", { + tl_sim <- getTrophicLevelBySpecies(sim) + expect_equal(dim(tl_sim), c(length(dimnames(sim@n)$time), no_sp)) +}) + + # getSSB ---- test_that("getSSB works", { ssb <- getSSB(sim)