R/inventory.R

In paleolimbot/pb210: Lead-210 dating utilities

#' Calculate cumulative lead-210 activity
#'
#' The cumulative content of lead-210 from the bottom of the core is the basis
#' for the constant rate of supply model, and is a required input to
#' [pb210_crs()]. There are as many ways to calculate this as there are people
#' who interpret lead-210 activities. This function models the bottom (below
#' measured lead-210), middle (between measured values of lead-210), and top
#' (above measured lead-210) with separate models to accomodate the variety of
#' methods. By default, lead-210 activity is estimated for samples in which it
#' was not measured by an exponential fit of lead-210 activity vs. cumulative dry mass
#' (bottom), and repeating the maximum value (top). Inventory for all values is
#' estimated using a trapezoidal approximation.
#'
#' @inheritParams pb210_cic
#' @param model_bottom A fit object that will be used to model activities below
#'   the last positive finite lead-210 activity. Use a constant to specify a custom
#'   deep inventory of lead-210. Can also be a lambda-function of
#'   `cumulative_dry_mass`, and `excess`. By default, this uses
#'   [pb210_fit_loglinear()] and [finite_tail_prop()] to generate an exponential
#'   fit using the last half of the finite values.
#' @param object An inventory calculator generated by
#'   [pb210_inventory_calculator()].
#'
#' @return [pb210_inventory()] returns a vector with [errors::errors()] of
#'   cumulative lead-210 activities for each sample (in Bq) that can be passed
#'   as the `inventory` to [pb210_crs()]. [pb210_inventory_calculator()] returns
#'   a fit object that can be used to calculate inventory given a
#'   `cumulative_dry_mass` and an `excess`.
#' @export
#'
#' @examples
#' fake_mass <- 1:10
#' fake_pb210 <- exp(5 - fake_mass)
#' pb210_inventory(fake_mass, fake_pb210)
#'
#' # compare with known inventory from integrating
#' # exp(5 - fake_mass) to +Inf
#' exp(-1 * fake_mass  + 5) / -(-1)
#'
pb210_inventory <- function(cumulative_dry_mass, excess,
                            model_top = max_finite(excess),
                            model_bottom = pb210_fit_loglinear(
                                cumulative_dry_mass, excess,
                                subset = ~finite_tail(..1, ..2, n_tail = 3)
                              )
                            ) {
  fit <- pb210_inventory_calculator(
    model_top = model_top,
    model_bottom = model_bottom
  )

  predict.inventory_calculator(
    fit,
    cumulative_dry_mass = cumulative_dry_mass,
    excess = excess
  )
}

#' @rdname pb210_inventory
#' @export
pb210_inventory_calculator <- function(model_top = ~max_finite(..2),
                                       model_bottom = ~pb210_fit_loglinear(
                                           ..1, ..2,
                                           ~finite_tail(..1, ..2, n_tail = 3)
                                         )
                                       ) {
  structure(
    list(
      model_top = pb210_as_fit(model_top),
      model_bottom = pb210_as_fit(model_bottom)
    ),
    class = c("inventory_calculator", "pb210_fit")
  )
}

#' @rdname pb210_inventory
#' @export
predict.inventory_calculator <- function(object, cumulative_dry_mass, excess, ...) {
  check_mass_and_activity(cumulative_dry_mass, without_errors(excess))

  use_errors <- any(is.finite(extract_errors(excess)))
  if (use_errors) {
    excess <- with_errors(excess)
    na <- with_errors(NA_real_)
  } else {
    excess <- without_errors(excess)
    na <- NA_real_
  }

  # model bottom/top functions expect errors as input
  # always need zero in cumulative dry mass so that the surface is approximated
  if (0 %in% cumulative_dry_mass) {
    data <- tibble::tibble(cumulative_dry_mass, excess = excess)
  } else {
    data <- tibble::tibble(
      cumulative_dry_mass = c(0, cumulative_dry_mass),
      excess = c(na, excess)
    )
  }

  model_bottom <- pb210_as_fit(object$model_bottom, data = data)
  model_top <- pb210_as_fit(object$model_top, data = data)

  # the rest of this function needs errors separate from values
  data$excess_sd <- extract_errors(data$excess)
  data$excess <-  without_errors(data$excess)

  finite_pb210_indices <- which(
    is.finite(data$excess) & (data$excess > 0)
  )

  first_finite_mass <- data$cumulative_dry_mass[min(finite_pb210_indices)] # kg
  last_finite_mass <- data$cumulative_dry_mass[max(finite_pb210_indices)] # kg

  # approximate the surface activities first
  data$excess[data$cumulative_dry_mass < first_finite_mass] <- stats::predict(
    model_top,
    tibble::tibble(x = data$cumulative_dry_mass[data$cumulative_dry_mass < first_finite_mass])
  )

  # recalculate finite indices, so the surface can be used in the trapezoidal approximation
  finite_pb210_indices <- which(
    is.finite(data$excess) & (data$excess > 0)
  )

  # the model exp(m*x + b),
  # integrated, is exp(m*x + b) / m, and because at x = infinity the integral is 0
  # the definite integral from [background] to infinity is exp(m * [background] + b) / -m
  coeffs <- stats::coefficients(model_bottom)
  deep_pb210 <- function(mass) unname(exp(coeffs["m"] * mass  + coeffs["b"]) / -coeffs["m"])

  # in the middle, approximate the cumulative sum as trapezoids using
  # finite values
  finite_pb210 <- data$excess[finite_pb210_indices]
  finite_mass <- data$cumulative_dry_mass[finite_pb210_indices]
  inventory_middle_interp <- integrate_trapezoid_no_error(
    finite_mass, finite_pb210,
    xout = data$cumulative_dry_mass,
    from = "right"
  ) + without_errors(deep_pb210(last_finite_mass))

  # combine the two methods to calculate inventory
  inventory <- ifelse(
    data$cumulative_dry_mass <= last_finite_mass,
    inventory_middle_interp,
    deep_pb210(data$cumulative_dry_mass)
  ) # Bq

  # Zero inventories aren't allowed because they get logged
  inventory[inventory <= 0] <- NA_real_

  # this error propogation is from R. Jack Cornett's spreadsheet, and should
  # be examined with more scrutiny
  if (use_errors) {
    pb210_relative_sd <- data$excess_sd / data$excess
    inventory_sd <- inventory * sqrt(2 * pb210_relative_sd^2 + 0.02^2)
    inventory <- with_errors(inventory, inventory_sd)
  }

  if (0 %in% cumulative_dry_mass) {
    inventory
  } else {
    inventory[-1]
  }
}