R/calibrate.R

In mizer: Dynamic Multi-Species Size Spectrum Modelling

#' Calibrate the model scale to match total observed biomass
#' 
#' `r lifecycle::badge("experimental")`
#' Given a MizerParams object `params` for which biomass observations are
#' available for at least some species via the `biomass_observed` column in the
#' species_params data frame, this function returns an updated MizerParams
#' object which is rescaled with [scaleModel()] so that the total biomass in
#' the model agrees with the total observed biomass.
#' 
#' Biomass observations usually only include individuals above a certain size.
#' This size should be specified in a biomass_cutoff column of the species
#' parameter data frame. If this is missing, it is assumed that all sizes are
#' included in the observed biomass, i.e., it includes larval biomass.
#' 
#' After using this function the total biomass in the model will match the
#' total biomass, summed over all species. However the biomasses of the
#' individual species will not match observations yet, with some species
#' having biomasses that are too high and others too low. So after this
#' function you may want to use [matchBiomasses()]. This is described in the
#' blog post at \url{https://blog.mizer.sizespectrum.org/posts/2021-08-20-a-5-step-recipe-for-tuning-the-model-steady-state/}.
#'
#' @param params A MizerParams object
#' @param ... Additional arguments passed to the method.
#' @return A MizerParams object. If no non-missing observed biomass values are
#'   provided, the original object is returned unchanged.
#' @export
#' @examples 
#' params <- NS_params
#' species_params(params)$biomass_observed <- 
#'     c(0.8, 61, 12, 35, 1.6, 20, 10, 7.6, 135, 60, 30, 78)
#' species_params(params)$biomass_cutoff <- 10
#' params2 <- calibrateBiomass(params)
#' plotBiomassObservedVsModel(params2)
calibrateBiomass <- function(params, ...)
    UseMethod("calibrateBiomass")

#' @export
calibrateBiomass.MizerParams <- function(params, ...) {
    calibrate_to(params, to = "biomass")
}

#' Calibrate the model scale to match total observed number
#'
#' `r lifecycle::badge("experimental")`
#' Given a MizerParams object `params` for which number observations are
#' available for at least some species via the `number_observed` column in the
#' species_params data frame, this function returns an updated MizerParams
#' object which is rescaled with [scaleModel()] so that the total number in
#' the model agrees with the total observed number.
#'
#' Number observations usually only include individuals above a certain size.
#' This size should be specified in a number_cutoff column of the species
#' parameter data frame. If this is missing, it is assumed that all sizes are
#' included in the observed number, i.e., it includes larval number.
#'
#' After using this function the total number in the model will match the
#' total number, summed over all species. However the numbers of the
#' individual species will not match observations yet, with some species
#' having numbers that are too high and others too low. So after this
#' function you may want to use [matchNumbers()]. This is described in the
#' blog post at \url{https://blog.mizer.sizespectrum.org/posts/2021-08-20-a-5-step-recipe-for-tuning-the-model-steady-state/}.
#'
#' @param params A MizerParams object
#' @param ... Additional arguments passed to the method.
#' @return A MizerParams object. If no non-missing observed number values are
#'   provided, the original object is returned unchanged.
#' @export
#' @examples
#' params <- NS_params
#' species_params(params)$number_observed <-
#'     c(0.8, 61, 12, 35, 1.6, 20, 10, 7.6, 135, 60, 30, 78)
#' species_params(params)$number_cutoff <- 10
#' params2 <- calibrateNumber(params)
calibrateNumber <- function(params, ...)
    UseMethod("calibrateNumber")

#' @export
calibrateNumber.MizerParams <- function(params, ...) {
    calibrate_to(params, to = "number")
}

#' Calibrate the model scale to match a total observation
#'
#' Internal implementation shared by [calibrateBiomass()] and
#' [calibrateNumber()]. Rescales the model with [scaleModel()] so that the
#' total over all observed species of the modelled quantity agrees with the
#' total of the observations. Species with no observation are left out of both
#' totals.
#'
#' @param params A MizerParams object.
#' @param to The type of observation, either "biomass" or "number".
#' @return A MizerParams object. If no non-missing observations are provided,
#'   the original object is returned unchanged.
#' @concept helper
#' @keywords internal
calibrate_to <- function(params, to = c("biomass", "number")) {
    cols <- observation_columns(to)
    observed <- params@species_params[[cols$observed]]
    if (is.null(observed) || all(is.na(observed))) {
        return(params)
    }
    observed_sel <- !is.na(observed)
    model <- model_observation(params, cols$to)
    scaleModel(params, factor = sum(observed[observed_sel]) /
                   sum(model[observed_sel]))
}

#' Change scale of the model
#'
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' The abundances in mizer and some rates depend on the size of the area to
#' which they refer. So they could be given per square meter or per square
#' kilometre or for an entire study area or any other choice of yours. This
#' function allows you to change the scale of the model by automatically
#' changing the abundances and rates accordingly.
#'
#' @details
#' If you rescale the model by a factor \eqn{c} then this function makes the
#' following rescalings in the params object:
#' \itemize{
#' \item The initial abundances are rescaled by \eqn{c}.
#' \item The search volume is rescaled by \eqn{1/c}.
#' \item The resource carrying capacity is rescaled by \eqn{c}
#' \item The maximum reproduction rate \eqn{R_{max}} is rescaled by
#'   \eqn{c}.
#' }
#' The effect of this is that the dynamics of the rescaled model are identical
#' to those of the unscaled model, in the sense that it does not matter whether
#' one first calls [scaleModel()] and then runs a simulation with
#' [project()] or whether one first runs a simulation and then rescales the
#' resulting abundances.
#'
#' Note that if you use non-standard resource dynamics or other components then
#' you may need to rescale additional parameters that appear in those dynamics.
#' 
#' In practice you will need to use some observations to set the scale for your
#' model. If you have biomass observations you can use [calibrateBiomass()],
#' if you have observed numbers you can use [calibrateNumber()].
#'
#' @param params A MizerParams object
#' @param factor The factor by which the scale is multiplied
#' @param ... Additional arguments passed to the method.
#'
#' @return The rescaled MizerParams object
#' @export
scaleModel <- function(params, factor, ...)
    UseMethod("scaleModel")

#' @export
scaleModel.MizerParams <- function(params, factor, ...) {
    params <- validParams(params)
    assert_that(is.number(factor),
                factor > 0)
    
    # Resource carrying capacity
    params@cc_pp <- params@cc_pp * factor
    
    # Search volume
    params@search_vol <- params@search_vol / factor
    
    sp <- params@species_params
    # Rmax
    # r_max is a deprecated spelling of R_max. Get rid of it.
    if ("r_max" %in% names(sp)) {
        sp$R_max <- sp$r_max
        sp$r_max <- NULL
        signal_info("R_max",
                    "The 'r_max' column has been renamed to 'R_max'.",
                    level = 1, unhandled = "show")
    }
    if ("R_max" %in% names(sp)) {
        sp$R_max <- sp$R_max * factor
    }
    if ("gamma" %in% names(sp)) {
        sp$gamma <- sp$gamma / factor
    }
    # `recalculate = FALSE` records the scaled parameters, so that a later
    # recalculation does not undo the rescaling, without recalculating the
    # rates: the search volume has already been scaled by hand above.
    species_params(params, recalculate = FALSE) <- sp
    
    # Initial values
    initial_n_other <- params@initial_n_other
    for (res in names(initial_n_other)) {
        initial_n_other[[res]] <- initial_n_other[[res]] * factor
    }
    initialN(params) <- params@initial_n * factor
    initialNResource(params) <- params@initial_n_pp * factor
    initialNOther(params) <- initial_n_other
    # Resource abundance coefficient
    params@resource_params$kappa <- params@resource_params$kappa * factor
    
    # community
    params@sc <- params@sc * factor
    
    params@time_modified <- lubridate::now()
    return(params)
}