ananke coverage - 63.81%

Files
Source

# DESCRIBE
#' @include AllGenerics.R
NULL

#' @export
#' @rdname describe
#' @aliases describe,CalibratedAges-method
setMethod(
  f = "describe",
  signature = signature(x = "CalibratedAges"),
  definition = function(x, calendar = get_calendar(),
                        level = 0.954, ...) {
    ## Get data
    lab <- labels(x)
    val <- x@values
    err <- x@errors
    crv <- x@curves
    reservoir_off <- x@reservoir_offsets
    reservoir_err <- x@reservoir_errors
    F14C <- x@F14C

    ## Laboratory code
    if (F14C) {
      txt_uncal <- "Sample %s contains %.0f +/- %.0f F14C,"
    } else {
      txt_uncal <- "Sample %s is dated to %.0f +/- %.0f BP,"
    }
    msg_uncal <- sprintf(txt_uncal, lab, val, err)

    ## Calibration results
    hdr <- as.list(interval_hdr(x, level = level), calendar = calendar)
    msg_cal <- lapply(
      X = hdr,
      FUN = function(x, calendar, level) {
        if (is.null(x)) return("but is out of the calibration range of")
        p <- if (NROW(x) > 1) sprintf(" (%.1f%%)", x$p * 100) else ""
        msg_hdr <- sprintf("[%.0f,%.0f]%s", x$start, x$end, p)
        txt_cal <- "calibrated to %s %s (%.1f%% HPD interval) with"
        sprintf(txt_cal, paste0(msg_hdr, collapse = " or "), calendar, level)
      },
      calendar = calendar@label,
      level = level * 100
    )

    ## Calibration curve
    txt_marine <- "marine reservoir offset: %.0f +/- %.0f; "
    msg_marine <- sprintf(txt_marine, reservoir_off, reservoir_err)
    msg_marine[reservoir_off == 0] <- ""
    txt_curve <- "%s (%s%s)."
    msg_curve <- sprintf(txt_curve, crv, msg_marine, cite_curve(crv))

    ## Text
    msg <- paste(msg_uncal, msg_cal, msg_curve, sep = " ")

    ## Software
    txt_soft <- "Calibration was computed with R %s.%s (R Core Team %s) and package ananke %s (Frerebeau %s)."
    date_soft <- utils::packageDate("ananke")
    date_soft <- if (is.na(date_soft)) Sys.Date() else date_soft
    msg_soft <- sprintf(txt_soft, R.version$major, R.version$minor,
                        R.version$year, utils::packageVersion("ananke"),
                        format(date_soft, format = "%Y"))

    ## Split string
    fill <- list(...)$fill
    if (!is.null(fill) && isTRUE(fill)) {
      msg <- lapply(
        X = msg,
        FUN = function(x, width) {
          split <- lapply(
            X = seq(from = 1, to = nchar(x), by = width),
            FUN = function(i) substr(x, i, i + width - 1)
          )
          paste0(split, collapse = "\n")
        },
        width = getOption("width")
      )
    }

    cat(unlist(msg), msg_soft, sep = "\n\n", ...)

    invisible(x)
  }
)

# PLOT
#' @include AllGenerics.R
NULL

## 14C =========================================================================
#' @export
#' @method plot CalibratedAges
plot.CalibratedAges <- function(x, calendar = get_calendar(),
                                interval = c("hdr", "credible"),
                                level = 0.954, sort = TRUE, decreasing = FALSE,
                                col.density = "grey", col.interval = "#77AADD",
                                main = NULL, sub = NULL,
                                ann = graphics::par("ann"),
                                axes = TRUE, frame.plot = TRUE,
                                panel.first = NULL, panel.last = NULL, ...) {
  ## Check
  c14_validate(x)
  interval <- match.arg(interval, several.ok = FALSE)

  ## Sort
  if (isTRUE(sort)) {
    mid <- order(mean(x, calendar = NULL), decreasing = !decreasing)
    x <- x[, mid, , drop = FALSE]
  }

  ## Graphical parameters
  fill.density <- grDevices::adjustcolor(col.density, alpha.f = 0.5)
  fill.interval <- grDevices::adjustcolor(col.interval, alpha.f = 0.5)

  ## Plot
  panel_density <- function(x, y, ...) {
    force(interval)
    force(level)
    force(fill.interval)

    tick_bottom <- min(y, na.rm = TRUE)
    tick_height <- tick_bottom + graphics::par("tcl") * graphics::strheight("M") * -1

    ## Keep only density > 0
    d0 <- which(y > tick_bottom)
    x <- x[d0]
    y <- y[d0]

    ## Draw density
    graphics::polygon(
      x = c(x, rev(x)),
      y = c(y, rep(tick_bottom, length(y))),
      border = NA,
      ...
    )

    ## Add interval
    if (level > 0) {
      y0 <- arkhe::scale_range(y)
      int <- switch(
        interval,
        hdr = arkhe::interval_hdr(x, y0, level = level),
        credible = {
          spl <- sample(x, size = length(x), replace = TRUE, prob = y0)
          arkhe::interval_credible(spl, level = level)
        }
      )

      for (i in seq_len(nrow(int))) {
        debut <- int[i, "start"]
        fin <- int[i, "end"]
        if (debut < fin) is_in_int <- which(x >= debut & x <= fin)
        else is_in_int <- which(x <= debut & x >= fin)
        xi <- x[is_in_int]
        yi <- y[is_in_int]
        graphics::polygon(
          x = c(xi[1], xi, xi[length(xi)]),
          y = c(tick_bottom, yi, tick_bottom),
          border = NA, col = fill.interval
        )
        graphics::segments(
          x0 = c(debut, debut, fin),
          x1 = c(debut, fin, fin),
          y0 = c(tick_bottom, tick_bottom, tick_bottom),
          y1 = c(tick_height, tick_bottom, tick_height),
          lend = 1
        )
      }
    }

    graphics::lines(x, y, col = "black")
  }

  ## Method for TimeSeries
  methods::callNextMethod(
    x, facet = "multiple",
    calendar = calendar,
    panel = panel_density,
    main = main, sub = sub, ann = ann, axes = axes,
    frame.plot = frame.plot,
    panel.first = panel.first,
    panel.last = panel.last,
    col = fill.density,
    ...
  )

  invisible(x)
}

#' @export
#' @rdname c14_plot
#' @aliases plot,CalibratedAges,missing-method
setMethod("plot", c(x = "CalibratedAges", y = "missing"), plot.CalibratedAges)

#' @export
#' @rdname ridgelines
#' @aliases ridgelines,CalibratedAges-method
setMethod(
  f = "ridgelines",
  signature = c(x = "CalibratedAges"),
  definition = function(x, calendar = get_calendar(),
                        interval = c("hdr", "credible"),
                        level = 0.954, fixed = TRUE, decreasing = FALSE,
                        col.density = "grey", col.interval = "#77AADD", ...) {
    ## Get data
    lab <- labels(x)

    ## Y scale
    if (isTRUE(fixed)) {
      mid <- mean(x, calendar = NULL)
      dy <- order(order(mid, decreasing = decreasing))
    } else {
      dy <- order(order(x@positions, decreasing = !decreasing))
    }
    for (j in seq_len(ncol(x))) x[, j, ] <- arkhe::scale_range(x[, j, ]) + dy[j]

    ## Permute
    x <- aion::flip(x)

    ## Plot
    plot(x, calendar = calendar, interval = interval, level = level,
         sort = FALSE, col.density = col.density,
         col.interval = col.interval,
         axes = FALSE, frame.plot = FALSE, yaxt = "n", ...)

    ## Construct Axis
    aion::year_axis(side = 1, format = TRUE, calendar = calendar,
                    xpd = NA)
    graphics::axis(side = 2, at = dy, labels = lab, las = 2,
                   lty = 0)
  }
)

## SPD =========================================================================
#' @export
#' @method plot CalibratedSPD
plot.CalibratedSPD <- function(x, calendar = get_calendar(),
                               main = NULL, sub = NULL,
                               ann = graphics::par("ann"),
                               axes = TRUE, frame.plot = TRUE,
                               panel.first = NULL, panel.last = NULL, ...) {
  ## Graphical parameters
  n <- NCOL(x)
  col <- list(...)$col %||% c("grey")
  if (length(col) != n) col <- rep(col, length.out = n)
  col <- grDevices::adjustcolor(col, alpha.f = 0.5)

  ## Plot
  panel_density <- function(x, y, ...) {
    graphics::polygon(
      x = c(x, rev(x)),
      y = c(y, rep(0, length(y))),
      border = NA,
      ...
    )
    graphics::lines(x, y, col = "black")
  }

  ## Method for TimeSeries
  methods::callNextMethod(
    x, facet = "multiple",
    calendar = calendar,
    panel = panel_density,
    main = main, sub = sub, ann = ann, axes = axes,
    frame.plot = frame.plot,
    panel.first = panel.first,
    panel.last = panel.last,
    col = col,
    ...
  )

  invisible(x)
}

#' @export
#' @rdname c14_plot
#' @aliases plot,CalibratedSPD,missing-method
setMethod("plot", c(x = "CalibratedSPD", y = "missing"), plot.CalibratedSPD)

## RECE ========================================================================
#' @export
#' @method plot RECE
plot.RECE <- function(x, calendar = get_calendar(), ...) {
  ## Binary array
  bin <- array(FALSE, dim = c(nrow(x), max(x), ncol(x)))
  for (j in seq_len(ncol(x))) {
    z <- x[, j, , drop = TRUE]
    for (i in seq_along(z)) {
      bin[i, z[i], j] <- z[i] > 0
    }
  }
  bin <- apply(X = bin, MARGIN = c(1, 2), FUN = sum)
  bin[bin == 0] <- NA

  ## Add annotation
  years <- aion::time(x, calendar = NULL)

  ## Plot
  graphics::image(
    x = years,
    y = seq_len(max(x)),
    z = log(bin),
    xlab = format(calendar),
    ylab = "Count",
    xaxt = "n",
    yaxt = "n",
    ...
  )

  ## Construct axes
  aion::year_axis(side = 1, format = TRUE, calendar = calendar,
                  current_calendar = NULL)
  graphics::axis(side = 2, at = seq_len(max(x)), las = 1)

  invisible(x)
}

#' @export
#' @rdname rec_plot
#' @aliases plot,RECE,missing-method
setMethod("plot", c(x = "RECE", y = "missing"), plot.RECE)


## Proxy =======================================================================
#' @export
#' @method plot ProxyRecord
plot.ProxyRecord <- function(x, calendar = get_calendar(),
                             iqr = TRUE,
                             xlab = NULL, ylab = NULL,
                             col = grDevices::hcl.colors(12, "YlOrRd", rev = TRUE),
                             col.mean = "black", col.iqr = col.mean,
                             lty.mean = 1, lty.iqr = 3,
                             lwd.mean = 2, lwd.iqr = lwd.mean, ...) {
  ## Get data
  years <- aion::time(x, calendar = NULL)
  z <- apply(
    X = x@density,
    MARGIN = 1,
    FUN = function(d) (d - min(d)) / max(d - min(d))
  )
  z[z == 0] <- NA

  ## Plot
  graphics::image(
    x = years,
    y = x@proxy,
    z = t(z),
    col = col,
    xaxt = "n",
    yaxt = "n",
    xlab = xlab %||% format(calendar),
    ylab = ylab %||% "Proxy",
    ...
  )

  ## Construct axes
  aion::year_axis(side = 1, format = TRUE, calendar = calendar,
                  current_calendar = NULL)
  graphics::axis(side = 2, las = 1)

  if (isTRUE(iqr)) {
    m <- mean(x)
    graphics::lines(x = years, y = m, col = col.mean,
                    lty = lty.mean, lwd = lwd.mean)

    q <- quantile(x, probs = c(0.25, 0.75))
    graphics::lines(x = years, y = q[, 1], col = col.iqr,
                    lty = lty.iqr, lwd = lwd.iqr)
    graphics::lines(x = years, y = q[, 2], col = col.iqr,
                    lty = lty.iqr, lwd = lwd.iqr)
  }

  invisible(x)
}

#' @export
#' @rdname proxy_plot
#' @aliases plot,ProxyRecord,missing-method
setMethod("plot", c(x = "ProxyRecord", y = "missing"), plot.ProxyRecord)

# 14C SPD
#' @include AllGenerics.R
NULL

#' @export
#' @rdname c14_spd
#' @aliases c14_spd,CalibratedAges-method
setMethod(
  f = "c14_spd",
  signature = "CalibratedAges",
  definition = function(object, normalize_date = FALSE, normalize_spd = FALSE) {
    ## Check
    c14_validate(object)

    dens <- t(object[, , 1, drop = TRUE])
    if (normalize_date) dens <- dens / rowSums(dens, na.rm = TRUE)
    spd <- colSums(dens, na.rm = TRUE)
    if (normalize_spd) spd <- spd / sum(spd, na.rm = TRUE)

    time_series <- aion::series(
      object = spd,
      time = aion::time(object, calendar = NULL)
    )
    .CalibratedSPD(time_series)
  }
)

# COERCION
#' @include AllGenerics.R
NULL

# To list ======================================================================
#' @export
#' @method as.list CalibratedIntervals
as.list.CalibratedIntervals <- function(x, ..., calendar = get_calendar()) {
  z <- as.data.frame(x, calendar = calendar)
  f <- factor(z$label, levels = unique(z$label))
  z$label <- NULL
  split(x = z, f = f)
}

#' @export
#' @rdname as.list
#' @aliases as.list,CalibratedIntervals-method
setMethod("as.list", "CalibratedIntervals", as.list.CalibratedIntervals)

# To data.frame ================================================================
#' @export
#' @method as.data.frame CalibratedAges
as.data.frame.CalibratedAges <- function(x, ...,
                                         calendar = get_calendar()) {
  dens <- x[, , 1, drop = TRUE]
  z <- data.frame(aion::time(x, calendar = calendar), dens)
  colnames(z) <- c("time", colnames(x) %||% paste0("X", seq_len(NCOL(x))))
  z
}

#' @export
#' @rdname as.data.frame
#' @aliases as.data.frame,CalibratedAges-method
setMethod("as.data.frame", "CalibratedAges", as.data.frame.CalibratedAges)

#' @export
#' @method as.data.frame CalibratedIntervals
as.data.frame.CalibratedIntervals <- function(x, ...,
                                              calendar = get_calendar()) {
  ## Build a data frame
  data.frame(
    label = labels(x),
    start = start(x, calendar = calendar),
    end = end(x, calendar = calendar),
    p = x@p
  )
}

#' @export
#' @rdname as.data.frame
#' @aliases as.data.frame,CalibratedIntervals-method
setMethod("as.data.frame", "CalibratedIntervals", as.data.frame.CalibratedIntervals)

#' @export
#' @method as.data.frame RECE
as.data.frame.RECE <- function(x, ..., calendar = get_calendar()) {
  dens <- x[, , 1, drop = TRUE]
  z <- data.frame(aion::time(x, calendar = calendar), dens)
  colnames(z) <- c("time", colnames(x) %||% paste0("X", seq_len(NCOL(x))))
  z
}

#' @export
#' @rdname as.data.frame
#' @aliases as.data.frame,RECE-method
setMethod("as.data.frame", "RECE", as.data.frame.RECE)

#' @export
#' @method as.data.frame ProxyRecord
as.data.frame.ProxyRecord <- function(x, ...,
                                      calendar = get_calendar()) {
  dens <- x[, , 1, drop = TRUE]
  z <- data.frame(aion::time(x, calendar = calendar), dens)
  colnames(z) <- c("time", colnames(x) %||% paste0("X", seq_len(NCOL(x))))
  z
}

#' @export
#' @rdname as.data.frame
#' @aliases as.data.frame,ProxyRecord-method
setMethod("as.data.frame", "ProxyRecord", as.data.frame.ProxyRecord)

# 14C UNCALIBRATION
#' @include AllGenerics.R
NULL

# Uncalibrate ==================================================================
#' @export
#' @rdname c14_uncalibrate
#' @aliases c14_uncalibrate,numeric-method
setMethod(
  f = "c14_uncalibrate",
  signature = c(object = "numeric"),
  definition = function(object, curves = "intcal20") {
    ## Validation
    n <- length(object)
    if (length(curves) != n) curves <- rep(curves, n)
    uncal <- rep(NA_real_, n)

    ## Calibration curve
    curve_data <- c14_curve(tolower(curves))

    for (i in seq_len(n)) {
      z <- stats::approx(
        x = curve_data[[i]][, 1],
        y = curve_data[[i]][, 2],
        xout = object[i],
        rule = 1
      )
      uncal[i] <- z$y
    }
    uncal
  }
)

#' @export
#' @rdname c14_uncalibrate
#' @aliases c14_uncalibrate,CalibratedAges-method
setMethod(
  f = "c14_uncalibrate",
  signature = c(object = "CalibratedAges"),
  definition = function(object, n = 10000, rounding = getOption("ananke.round"),
                        ...) {

    method <- list(...)$method %||% c("L-BFGS-B")

    ## Function to be optimized
    opt_fun <- function(param, curve, samples, ...) {
      cal <- c14_calibrate(
        values = round(param[1]),
        errors = round(param[2]),
        curves = curve
      )
      s1 <- samples
      s2 <- c14_sample(cal, n = n, calendar = BP())

      dens1 <- stats::density(s1, from = min(c(s1, s2)), to = max(c(s1, s2)))$y
      dens2 <- stats::density(s2, from = min(c(s1, s2)), to = max(c(s1, s2)))$y
      dens1 <- dens1 / sum(dens1)
      dens2 <- dens2 / sum(dens2)
      int <- dens1 * log(dens1 / dens2)
      int <- int[dens1 > 0]
      int <- int[!is.infinite(int)]
      sum(int)
    }

    n <- NCOL(object)
    curves <- object@curves

    ## Initial values for the parameters to be optimized over
    spl <- c14_sample(object, n = n, calendar = BP())
    med <- apply(X = spl, MARGIN = 2, FUN = stats::median)
    init_mean <- c14_uncalibrate(med, curves = curves)
    init_sd <- apply(X = spl, MARGIN = 2, FUN = stats::sd)

    opt_mean <- opt_sd <- rep(NA_real_, n)
    for (i in seq_len(n)) {
      cal_spl <- spl[, i]

      opt_run <- stats::optim(
        par = c(init_mean[[i]], init_sd[[i]]),
        fn = opt_fun,
        method = method,
        curve = curves[[i]],
        samples = cal_spl,
        lower = c(-Inf, 1),
        ...
      )
      opt_mean[[i]] <- opt_run$par[1]
      opt_sd[[i]] <- opt_run$par[2]
    }

    ## Rounding
    if (identical(rounding, "stuiver")) {
      opt_mean <- round_values_stuiver(opt_mean)
      opt_sd <- round_errors_stuiver(opt_sd)
    }

    data.frame(mean = opt_mean, sd = opt_sd)
  }
)

# REC
#' @include AllGenerics.R
NULL

#' @export
#' @rdname c14_ensemble
#' @aliases c14_ensemble,CalibratedAges-method
setMethod(
  f = "c14_ensemble",
  signature = "CalibratedAges",
  definition = function(object, from = NULL, to = NULL, by = 10, n = 100,
                        calendar = BP(), progress = getOption("ananke.progress")) {
    ## Check
    c14_validate(object)

    ## Get data
    rd <- aion::time(object, calendar = NULL)
    if (is.null(from)) from <- aion::start(object, calendar = calendar)
    if (is.null(to)) to <- aion::end(object, calendar = calendar)
    grid_years <- seq(from = from, to = to, by = by * sign(to - from))
    grid_rd <- aion::fixed(grid_years, calendar = calendar)

    ## Align 14C date densities onto the grid
    c14_dens <- object[, , 1, drop = TRUE]
    c14_aligned <- apply(
      X = c14_dens,
      MARGIN = 2,
      FUN = function(y, x, grid) {
        fun <- stats::approxfun(x = x, y = y)
        fun(grid)
      },
      x = rd,
      grid = grid_rd
    )
    c14_aligned[is.na(c14_aligned)] <- 0

    ## Build matrix to store the RECE
    count <- matrix(data = 0, nrow = n, ncol = length(grid_rd))
    colnames(count) <- grid_rd

    progress_bar <- progress
    if (progress_bar) pbar <- utils::txtProgressBar(max = n, style = 3)

    n_seq <- seq_len(n)
    for (i in n_seq) {
      ## Sample
      spl <- apply(
        X = c14_aligned,
        MARGIN = 2,
        FUN = function(x, grid) {
          if (sum(x) == 0) return(NA)
          sample(grid, size = 1, prob = x)
        },
        grid = grid_rd
      )
      ## Count
      tbl <- unclass(table(spl)) # Named integer vector
      count[i, names(tbl)] <- tbl
      if (progress_bar) utils::setTxtProgressBar(pbar, i)
    }
    count[is.na(count)] <- 0
    if (progress_bar) close(pbar)

    ## Return an RECE object
    ts <- aion::series(
      object = t(count),
      time = grid_rd
      # names = colnames(object)
    )
    .RECE(ts)
  }
)

# AGE-DEPTH MODELING
#' @include AllGenerics.R
NULL

#' @export
#' @rdname proxy_ensemble
#' @aliases proxy_ensemble,numeric-method
setMethod(
  f = "proxy_ensemble",
  signature = c("numeric"),
  definition = function(positions, proxy_values, proxy_errors, proxy_step,
                        time_values, time_errors, calendar,
                        from = NULL, to = NULL, by = NULL, n = 30,
                        progress = getOption("ananke.progress"),
                        verbose = getOption("ananke.verbose")) {
    ## Validation
    k <- length(positions)
    if (length(proxy_errors) == 1) proxy_errors <- rep(proxy_errors, k)
    arkhe::assert_positive(positions)
    arkhe::assert_decreasing(positions)
    arkhe::assert_length(proxy_values, k)
    arkhe::assert_length(proxy_errors, k)
    arkhe::assert_length(proxy_step, 1)
    arkhe::assert_length(time_values, k)
    arkhe::assert_length(time_errors, k)

    ## Missing values
    if (is.null(from)) from <- time_values[[1L]]
    if (is.null(to))   to <- time_values[[k]]
    if (is.null(by)) {
      grid <- getOption("ananke.grid")
      by <- ((to - from) / (grid - 1))
    }
    if (from > to && by > 0) by <- by * -1

    ## Build a matrix to contain the p(t|zi) densities
    ## Rows will refer to depth
    ## Columns will refer to the time density
    if (verbose) cat(tr_("Computing p(t|zi) densities..."), sep = "\n")
    t_grid <- seq(from = from, to = to, by = by)

    t_z <- .mapply(
      FUN = function(min, max, x) {
        stats::dunif(x = x, min = min, max = max)
      },
      dots = list(
        min = time_values - 1 * time_errors,
        max = time_values + 1 * time_errors
      ),
      MoreArgs = list(x = t_grid)
    )
    t_z <- do.call(rbind, t_z)
    t_z <- t(t_z)

    ## Build a matrix to contain the p(x|zi) densities
    ## Rows will refer to depth
    ## Columns will refer to the proxy density
    if (verbose) cat(tr_("Computing p(x|zi) densities..."), sep = "\n")
    d <- 2 * max(proxy_errors)
    x_range <- range(c(range(proxy_values) - d, range(proxy_values) + d))
    x_grid <- seq(from = x_range[[1L]], to = x_range[[2L]], by = proxy_step)

    x_z <- .mapply(
      FUN = function(mean, sd, x) {
        stats::dnorm(x = x, mean = mean, sd = sd)
      },
      dots = list(
        mean = proxy_values,
        sd = proxy_errors
      ),
      MoreArgs = list(x = x_grid)
    )
    x_z <- do.call(rbind, x_z)
    x_z <- t(x_z)

    ## Estimate the weighted average density function
    ## (for a given proxy measurement at a given time)
    ## Eq. 4 of Boers et al. 2017
    if (verbose) cat(tr_("Computing p(x|t) densities..."), sep = "\n")

    z <- length(positions)
    ri <- vapply(
      X = 2:(z - 1),
      FUN = function(x, y) (y[x + 1] - y[x - 1]) / 2,
      FUN.VALUE = numeric(1),
      y = positions
    )
    r <- abs(c(positions[2] - positions[1], ri, positions[z] - positions[z - 1]))
    r_t <- matrix(data = r, nrow = nrow(t_z), ncol = ncol(t_z), byrow = TRUE)
    r_x <- matrix(data = r, nrow = nrow(x_z), ncol = ncol(x_z), byrow = TRUE)

    x_t <- t(tcrossprod(r_x * x_z, t_z)) / rowSums(r_t * t_z)

    ## Create an ensemble of potential proxy records
    if (verbose) cat(tr_("Sampling proxy records..."), sep = "\n")
    Y <- matrix(data = 0, nrow = length(t_grid), ncol = n)

    if (progress) pb <- utils::txtProgressBar(min = 0, max = n, style = 3)

    n_seq <- seq_len(n)
    for (i in n_seq) {
      ## Sample
      Y[, i] <- apply(
        X = x_t,
        MARGIN = 1,
        FUN = function(x, g) sample(g, size = 1, prob = x),
        g = x_grid
      )
      if (progress) utils::setTxtProgressBar(pb, i)
    }

    if (progress) close(pb)

    ts <- aion::series(object = Y, time = t_grid, calendar = calendar)
    .ProxyRecord(ts, density = x_t, proxy = x_grid)
  }
)

# STATISTICS
#' @include AllGenerics.R
NULL

# Quantile =====================================================================
#' Quantiles of a Density Estimate
#'
#' @param x A [`numeric`] vector giving the n coordinates of the points where
#'  the density is estimated.
#' @param y A `numeric` [`matrix`] of density estimates (each column is a
#'  density estimate).
#' @param probs A [`numeric`] vector of probabilities with values in
#'  \eqn{[0,1]}.
#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
#'  computation proceeds?
#' @param ... Currently not used.
#' @return
#'  A `numeric` [`matrix`] containing the quantiles.
#' @keywords internal
#' @noRd
quantile_density <- function(x, y, probs = seq(0, 1, 0.25), na.rm = FALSE, ...) {
  eps <- 100 * .Machine$double.eps
  if (anyNA(probs) || any(probs < -eps | probs > 1 + eps))
    stop(sprintf("%s outside [0,1]", sQuote("probs")))

  q <- apply(
    X = y,
    MARGIN = 2,
    FUN = function(y, x, probs, na.rm) {
      np <- length(probs)
      qs <- rep(NA_real_, np)
      if (na.rm) {
        i <- !is.na(x) & !is.na(y)
        x <- x[i]
        y <- y[i]
      }
      if (np > 0) {
        nn <- length(x)
        Fx <- cumsum(y * c(0, diff(x)))
        Fx <- Fx / Fx[nn]
        for (j in seq_len(np)) {
          ii <- min(which(Fx >= probs[j]))
          if (!is.na(ii) && ii >= 1 && ii <= nn) qs[j] <- x[ii]
        }
        qs
      }
    },
    x = x,
    probs = probs,
    na.rm = na.rm
  )

  if (!is.null(dim(q))) {
    q <- t(q)
    colnames(q) <- paste0(round(probs * 100, digits = 0), "%")
  }
  q
}

#' @export
#' @method quantile CalibratedAges
quantile.CalibratedAges <- function(x, probs = seq(0, 1, 0.25), na.rm = FALSE,
                                    ..., calendar = get_calendar()) {
  quantile_density(x = aion::time(x, calendar = calendar), y = x,
                   probs = probs, na.rm = na.rm, ...)
}

#' @export
#' @rdname quantile
#' @aliases quantile,CalibratedAges,missing-method
setMethod("quantile", c(x = "CalibratedAges"), quantile.CalibratedAges)

#' @export
#' @method quantile ProxyRecord
quantile.ProxyRecord <- function(x, probs = seq(0, 1, 0.25), na.rm = FALSE, ...) {
  quantile_density(x = x@proxy, y = t(x@density),
                   probs = probs, na.rm = na.rm, ...)
}

#' @export
#' @rdname quantile
#' @aliases quantile,ProxyRecord,missing-method
setMethod("quantile", c(x = "ProxyRecord"), quantile.ProxyRecord)

# Median =======================================================================
#' @export
#' @method median CalibratedAges
median.CalibratedAges <- function(x, na.rm = FALSE, ...,
                                  calendar = get_calendar()) {
  quantile_density(x = aion::time(x, calendar = calendar), y = x,
                   probs = 0.5, na.rm = na.rm, ...)
}

#' @export
#' @rdname median
#' @aliases median,CalibratedAges,missing-method
setMethod("median", c(x = "CalibratedAges"), median.CalibratedAges)

# Mean =========================================================================
#' Mean of a Density Estimate
#'
#' @param x A [`numeric`] vector giving the n coordinates of the points where
#'  the density is estimated.
#' @param y A `numeric` [`matrix`] of density estimates (each column is a
#'  density estimate).
#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
#'  computation proceeds?
#' @param ... Currently not used.
#' @return A [`numeric`] vector.
#' @keywords internal
#' @noRd
mean_density <- function(x, y, na.rm = FALSE, ...) {
  apply(
    X = y,
    MARGIN = 2,
    FUN = function(w, x, na.rm) {
      if (na.rm) {
        i <- !is.na(w) & !is.na(x)
        x <- x[i]
        w <- w[i]
      }
      stats::weighted.mean(x = x, w = w)
    },
    x = x,
    na.rm = na.rm
  )
}

#' @export
#' @method mean CalibratedAges
mean.CalibratedAges <- function(x, na.rm = FALSE, ...,
                                calendar = get_calendar()) {
  mean_density(x = aion::time(x, calendar = calendar), y = x, na.rm = na.rm, ...)
}

#' @export
#' @rdname mean
#' @aliases mean,CalibratedAges,missing-method
setMethod("mean", c(x = "CalibratedAges"), mean.CalibratedAges)

#' @export
#' @method mean ProxyRecord
mean.ProxyRecord <- function(x, na.rm = FALSE, ...) {
  mean_density(x = x@proxy, y = t(x@density), na.rm = na.rm, ...)
}

#' @export
#' @rdname mean
#' @aliases mean,ProxyRecord,missing-method
setMethod("mean", c(x = "ProxyRecord"), mean.ProxyRecord)

# Summary ======================================================================
#' @export
#' @method summary CalibratedAges
summary.CalibratedAges <- function(object, ..., digits = NULL,
                                   calendar = get_calendar()) {
  mid <- matrix(mean(object, calendar = get_calendar()), ncol = 1)
  qq <- quantile(object, probs = c(0.25, 0.50, 0.75), calendar = get_calendar())

  sms <- t(cbind(qq[, 1L:2L, drop = FALSE], mid, qq[, 3L, drop = FALSE]))
  sms <- format(sms, digits = digits)

  lbs <- format(c(tr_("1st Qu."), tr_("Median"), tr_("Mean"), tr_("3rd Qu.")))
  lbs <- matrix(lbs, nrow = nrow(sms), ncol = ncol(sms))
  z <- paste0(lbs, ":", sms, "  ")
  nm <- format(colnames(sms), justify = "centre", width = max(nchar(z)))

  dim(z) <- dim(sms)
  dimnames(z) <- list(rep("", nrow(sms)), nm)
  attr(z, "class") <- c("table")

  z
}

#' @export
#' @rdname summary
#' @aliases summary,CalibratedAges-method
setMethod("summary", c(object = "CalibratedAges"), summary.CalibratedAges)

# COMBINE 14C
#' @include AllGenerics.R
NULL

#' @export
#' @rdname c14_combine
#' @aliases c14_combine,numeric,numeric-method
setMethod(
  f = "c14_combine",
  signature = c(values = "numeric", errors = "numeric"),
  definition = function(values, errors, groups = NULL) {
    ## Validation
    n <- length(values)
    if (is.null(groups)) groups <- "X"
    if (length(groups) == 1) groups <- rep(groups, n)

    ## Empty groups must be treated as NA
    groups[which(!nzchar(groups))] <- NA
    groups <- factor(x = groups, levels = unique(groups))

    arkhe::assert_missing(values)
    arkhe::assert_missing(errors)
    arkhe::assert_length(errors, n)
    arkhe::assert_length(groups, n)

    ## Groups with only one date must be treated as NA
    counts <- table(groups)
    one <- groups %in% names(counts)[counts == 1]

    # NA group will be removed
    # We need to keep isolated dates
    k <- one | is.na(groups)
    solo <- NULL
    if (any(k)) {
      solo <- data.frame(
        groups = as.character(groups[k]),
        values = values[k],
        errors = errors[k],
        chi2 = NA_real_,
        p = NA_real_
      )
    }

    combined <- NULL
    if (!all(k)) {
      groups[k] <- NA
      groups <- droplevels(groups)

      ## split() removes NA group
      values <- split(values, f = groups)
      errors <- split(errors, f = groups)
      cmbn <- mapply(
        FUN = combine,
        values = values,
        errors = errors,
        SIMPLIFY = FALSE
      )
      combined <- data.frame(names(cmbn), do.call(rbind, cmbn))
      colnames(combined) <- c("groups", "values", "errors", "chi2", "p")
    }

    final <- rbind(solo, combined, make.row.names = FALSE)
    final
  }
)

combine <- function(values, errors) {
  ## On calcule la moyenne pondérée
  w <- 1 / errors^2 # Facteur de pondération
  moy <- stats::weighted.mean(x = values, w = w)

  ## On calcule l'incertitude associée à la moyenne pondérée
  err <- sum(1 / errors^2)^(-1 / 2)

  ## On calcule la statistique du test
  chi2 <- sum(((values - moy) / errors)^2)

  ## On calcule la valeur-p
  p <- 1 - stats::pchisq(chi2, df = length(values))

  ## On stocke les résultats
  c(moy, err, chi2, p)
}

# F14C <> BP14C
#' @include AllGenerics.R
NULL

#' @export
#' @rdname F14C
#' @aliases c14_f14c,numeric,numeric-method
setMethod(
  f = "c14_f14c",
  signature = c(values = "numeric", errors = "numeric"),
  definition = function(values, errors, tau = 8033) {
    values <- exp(values / -tau)
    sigma <- values * errors / tau
    data.frame(value = values, error = sigma)
  }
)

#' @export
#' @rdname F14C
#' @aliases f14c_c14,numeric,numeric-method
setMethod(
  f = "f14c_c14",
  signature = c(values = "numeric", errors = "numeric"),
  definition = function(values, errors, tau = 8033, asymmetric = FALSE,
                        rounding = getOption("ananke.round")) {
    ## Validation
    choices <- c("none", "stuiver")
    rounding <- match.arg(rounding, choices = choices, several.ok = FALSE)

    z <- values

    ## van der Plicht and Hogg 2006, p. 239
    below_2sigma <- z < 2 * errors
    values[below_2sigma] <- values[below_2sigma] + 2 * errors[below_2sigma]

    below_zero <- z < 0
    values[below_zero] <- 2 * errors[below_zero]

    ## van der Plicht and Hogg 2006, eq. 6
    ## Bronk Ramsey 2008, p. 260
    ages <- -tau * log(values)
    sigma <- tau * errors / values

    sigma_plus <- sigma_minus <- sigma
    if (asymmetric) {
      sigma_plus <- - tau * log(values - errors) - ages
      sigma_minus <- ages + tau * log(values + errors)
    }

    sigma_plus[below_zero | below_2sigma] <- NA_real_
    sigma_minus[below_zero | below_2sigma] <- NA_real_

    ## The reported age can be no older than the radiocarbon age of the
    ## fraction modern of the sample plus it's error
    # limit <- -tau * log(values + errors)
    # ages[ages > limit] <- limit[ages > limit]

    ## Rounding
    if (identical(rounding, "stuiver")) {
      ages <- round_values_stuiver(ages)
      sigma_plus <- round_errors_stuiver(sigma_plus)
      sigma_minus <- round_errors_stuiver(sigma_minus)
    }

    data.frame(age = ages, plus = sigma_plus, minus = sigma_minus)
  }
)

## Stuiver & Polach (1977), p. 362
round_any <- function(x, accuracy, f = round) {
  f(x / accuracy) * accuracy
}
round_values_stuiver <- function(x) {
  y <- x[!is.na(x)]
  y[y < 1000] <- round_any(y[y < 1000], 5)
  y[1000 <= y & y <= 9999] <- round_any(y[1000 <= y & y < 9999], 10)
  y[10000 <= y & y <= 20000] <- round_any(y[10000 <= y & y <= 20000], 50)
  y[y > 20000] <- round_any(y[y > 20000], 100)
  x[!is.na(x)] <- y
  x
}
round_errors_stuiver <- function(x) {
  y <- x[!is.na(x)]
  y[y < 50] <- round_any(y[y < 50], 5)
  y[50 <= y & y <= 1000] <- round_any(y[50 <= y & y <= 1000], 10)
  y[y > 1000] <- round_any(y[y > 1000], 100)
  x[!is.na(x)] <- y
  x
}

# 14C CALIBRATION
#' @include AllGenerics.R
NULL

#' @export
#' @rdname c14_calibrate
#' @aliases c14_calibrate,numeric,numeric-method
setMethod(
  f = "c14_calibrate",
  signature = c(values = "numeric", errors = "numeric"),
  definition = function(values, errors, curves = "intcal20",
                        names = NULL, positions = NULL,
                        reservoir_offsets = 0, reservoir_errors = 0,
                        from = 55000, to = 0, resolution = 1,
                        normalize = TRUE, F14C = FALSE,
                        method = c("student", "normal"),
                        dfs = 100, drop = TRUE, eps = 1e-06,
                        verbose = getOption("ananke.verbose")) {
    ## Validation
    method <- match.arg(method)
    n <- length(values)
    if (is.null(names)) names <- paste0("X", seq_len(n))
    if (is.null(positions)) positions <- values
    if (length(curves) == 1)
      curves <- rep(curves, n)
    if (length(reservoir_offsets) == 1)
      reservoir_offsets <- rep(reservoir_offsets, n)
    if (length(reservoir_errors) == 1)
      reservoir_errors <- rep(reservoir_errors, n)
    if (length(dfs) == 1)
      dfs <- rep(dfs, n)

    arkhe::assert_missing(values)
    arkhe::assert_missing(errors)
    arkhe::assert_unique(names)
    arkhe::assert_length(errors, n)
    arkhe::assert_length(names, n)
    arkhe::assert_length(curves, n)
    arkhe::assert_length(reservoir_offsets, n)
    arkhe::assert_length(reservoir_errors, n)
    arkhe::assert_length(dfs, n)

    ## Calibration time range
    cal_range <- seq(from = from, to = to, by = -resolution)

    ## Calibration curve
    curves <- tolower(curves)
    curve_range <- approx_curve(unique(curves), out = cal_range, F14C = F14C)

    ## Marine reservoir offset
    values <- values - reservoir_offsets
    errors <- sqrt(errors^2 + reservoir_errors^2)

    ## Calibrate
    calibrate_fun <- if (F14C) calibrate_F14C else calibrate_BP14C
    dens <- vector(mode = "list", length = n)
    status <- integer(n)
    for (i in seq_len(n)) {
      d <- calibrate_fun(
        x = values[i],
        error = errors[i],
        mu = curve_range[[curves[i]]]$mu,
        tau = curve_range[[curves[i]]]$tau,
        df = dfs[i],
        method = method
      )

      ## Check
      if (anyNA(d)) {
        d[is.na(d)] <- 0
        msg <- tr_("Consider changing the calibration time range to a narrower interval.")
        if (verbose) message(msg)
      }

      max_cal <- curve_range[[curves[i]]]$max
      min_cal <- curve_range[[curves[i]]]$min

      if (values[[i]] >= max_cal || values[[i]] <= min_cal) {
        ## L'age à calibrer est hors de l'étendue de la courbe de calibration
        if (verbose) warning(print_out(names[[i]], maybe = FALSE), call. = FALSE)
        status[i] <- 1L
      } else if (d[1] > eps || d[length(d)] > eps) {
        ## L'age calibré est en partie hors de l'étendue de la courbe
        if (verbose) warning(print_out(names[[i]], maybe = TRUE), call. = FALSE)
        status[i] <- 2L
      }

      d[d < eps] <- 0
      dens[[i]] <- d
    }

    ## Build matrix
    dens <- do.call(rbind, dens)

    ## Normalize
    if (F14C && !normalize) normalize <- TRUE
    if (normalize) {
      dens <- dens / rowSums(dens, na.rm = TRUE)
      dens[dens < eps] <- 0
      dens <- dens / rowSums(dens, na.rm = TRUE)
    }

    ## Drop
    if (drop) {
      keep_zero <- colSums(dens, na.rm = TRUE) > 0
      keep_from <- which.max(keep_zero) # First TRUE
      keep_to <- length(keep_zero) - which.max(rev(keep_zero)) + 1 # Last TRUE
      keep <- seq(from = keep_from, to = keep_to, by = 1)
      dens <- dens[, keep, drop = FALSE]
      cal_range <- cal_range[keep]
    }

    time_series <- aion::series(
      object = t(dens),
      time = cal_range,
      calendar = BP(),
      names = names
    )
    .CalibratedAges(
      time_series,
      values = values,
      errors = errors,
      curves = curves,
      reservoir_offsets = reservoir_offsets,
      reservoir_errors = reservoir_errors,
      F14C = F14C,
      positions = positions,
      status = status
    )
  }
)

calibrate_BP14C <- function(x, error, mu, tau, df = 100, method = "student") {
  tau <- error^2 + tau^2
  dens <- switch(
    method,
    student = stats::dt(x = (x - mu) / sqrt(tau), df = df),
    normal = stats::dnorm(x, mean = mu, sd = sqrt(tau)),
    stop(sprintf("Unknown method %s.", dQuote(method)))
  )
  dens
}
calibrate_F14C <- function(x, error, mu, tau, ...) {
  p1 <- (x - mu)^2
  p2 <- 2 * (error^2 + tau^2)
  p3 <- sqrt(error^2 + tau^2)
  dens <- exp(-p1 / p2) / p3
  dens
}

# 14C CALIBRATION CHECK
#' @include AllGenerics.R
NULL

#' Check Calibrated Radiocarbon Dates
#'
#' @param object A [`CalibratedAges-class`] object.
#' @param verbose A [`logical`] scalar: should extra information be reported?
#' @return
#'  `c14_validate()` is called it for its side-effects: it prints
#'  [warning messages][warnings()] if calibrated ages are (partially) out of
#'  calibration range. Invisibly returns `x`.
#' @example inst/examples/ex-14c-calibrate.R
#' @author N. Frerebeau
#' @keywords internal
c14_validate <- function(object, verbose = getOption("ananke.verbose")) {
  status <- object@status
  lab <- labels(object)

  if (isTRUE(verbose)) {
    if (any(status == 1L)) {
      is_out <- which(status == 1L)
      warn <- print_out(lab[is_out], maybe = FALSE)
      for (w in warn) warning(w, call. = FALSE)
    }
    if (any(status == 2L)) {
      is_out <- which(status == 2L)
      warn <- print_out(lab[is_out], maybe = TRUE)
      for (w in warn) warning(w, call. = FALSE)
    }
  }

  invisible(object)
}

print_out <- function(label, maybe = FALSE) {
  if (maybe) {
    status <- tr_("Date %s may extent out of calibration range.")
  } else {
    status <- tr_("Date %s is out of calibration range.")
  }
  sprintf(status, dQuote(label))
}

# GEOLOGICAL MODEL AGE
#' @include AllGenerics.R
NULL

#' @export
#' @rdname pb_age
#' @aliases pb_age,numeric,numeric,numeric-method
setMethod(
  f = "pb_age",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, t0 = 3.8,
                        x_star = 18.75, y_star = 15.63, z_star = 38.86,
                        mu = 9.66, kappa = 3.90, th232 = 0.049475,
                        u238 = 0.155125, u235 = 0.98485, u238_235 = 137.79,
                        tolerance = sqrt(.Machine$double.eps),
                        stop = 100) {
    ## Validation
    n <- length(x)
    arkhe::assert_length(y, n)
    arkhe::assert_length(z, n)

    x0 <- x_star - mu * (exp(u238 * t0) - 1) # 206Pb/204Pb at t0
    y0 <- y_star - mu / u238_235 * (exp(u235 * t0) - 1) # 207Pb/204Pb at t0
    z0 <- z_star - mu * kappa * (exp(th232 * t0) - 1) # 208Pb/204Pb at t0

    t_model <- vector(mode = "numeric", length = n)
    f_x <- vector(mode = "numeric", length = n)
    for (i in seq_along(t_model)) {
      tt <- 0
      f <- 1

      iter <- 0
      while (abs(f) > tolerance) { # Newton loop
        slope <- (y[i] - y0) / (x[i] - x0)
        u <- exp(u235 * t0) - exp(u235 * tt)
        v <- exp(u238 * t0) - exp(u238 * tt)
        uprime <- -u235 * exp(u235 * tt)
        vprime <- -u238 * exp(u238 * tt)

        f <- u / v - slope * u238_235
        df <- (1 / v) * (uprime - u / v * vprime) # df/dtt
        dtt <- f / df # Newton-Raphson increment
        tt <- tt - dtt

        # Loop counter
        iter <- iter + 1
        if (iter >= stop) {
          msg <- tr_("Convergence not reached (possible infinite loop).")
          warning(msg, call. = FALSE)
          break
        }
      }

      t_model[i] <- tt
      f_x[i] <- f # Residual (should be zero)
    }

    mu_i <- (x - x0) / (exp(u238 * t0) - exp(u238 * t_model))
    kappa_i <- (z - z0) / (exp(th232 * t0) - exp(th232 * t_model)) / mu_i

    data.frame(
      age = t_model * 1000, # Use Ma instead of Ga
      mu = mu_i,
      kappa = kappa_i,
      residual = f_x
    )
  }
)

#' @export
#' @rdname pb_age
#' @aliases pb_age,list,missing,missing-method
setMethod(
  f = "pb_age",
  signature = c(x = "list", y = "missing", z = "missing"),
  definition = function(x, t0 = 3.8,
                        x_star = 18.75, y_star = 15.63, z_star = 38.86,
                        mu = 9.66, kappa = 3.90, th232 = 0.049475,
                        u238 = 0.155125, u235 = 0.98485, u238_235 = 137.79,
                        tolerance = sqrt(.Machine$double.eps), stop = 100) {
    ## Validation
    x <- grDevices::xyz.coords(x)
    methods::callGeneric(x = x$x, y = x$y, z = x$z, t0 = t0,
                         x_star = x_star, y_star = y_star, z_star = z_star,
                         mu = mu, kappa = kappa, th232 = th232,
                         u238 = u238, u235 = u235, u238_235 = u238_235,
                         tolerance = tolerance, stop = stop)
  }
)

# INTERVAL ESTIMATION
#' @include AllGenerics.R
NULL

# HDR ==========================================================================
#' @export
#' @rdname interval_hdr
#' @aliases interval_hdr,CalibratedAges,missing-method
setMethod(
  f = "interval_hdr",
  signature = c(x = "CalibratedAges", y = "missing"),
  definition = function(x, level = 0.954, ...) {
    ## Check
    c14_validate(x)

    hdr <- apply(
      X = x,
      MARGIN = 2,
      FUN = function(y, x, level, ...) {
        if (all(is.na(y))) return(NULL)
        arkhe::interval_hdr(x, as.numeric(y), level, ...)
      },
      x = aion::time(x, calendar = NULL),
      level = level,
      simplify = FALSE
    )
    n <- vapply(X = hdr, FUN = nrow, FUN.VALUE = integer(1))
    lab <- rep(labels(x), n)

    hdr <- do.call(rbind, hdr)
    .CalibratedIntervals(
      .Id = lab,
      .Start = aion::as_fixed(hdr[, 1]),
      .End = aion::as_fixed(hdr[, 2]),
      p = unname(hdr[, 3])
    )
  }
)

# Credible =====================================================================
#' @export
#' @rdname interval_credible
#' @aliases interval_credible,CalibratedAges,missing-method
setMethod(
  f = "interval_credible",
  signature = c(x = "CalibratedAges"),
  definition = function(x, level = 0.954, n = 100, ...) {
    ## Check
    c14_validate(x)

    spl <- c14_sample(x, n = n, calendar = NULL)
    crd <- apply(
      X = spl,
      MARGIN = 2,
      FUN = arkhe::interval_credible,
      level = level,
      simplify = FALSE
    )

    crd <- do.call(rbind, crd)
    .CalibratedIntervals(
      .Id = labels(x),
      .Start = aion::as_fixed(crd[, 1]),
      .End = aion::as_fixed(crd[, 2]),
      p = unname(crd[, 3])
    )
  }
)

# 14C CALIBRATION CURVE
#' @include AllGenerics.R
NULL

# Get curve ====================================================================
#' @export
#' @rdname c14_curve
#' @aliases c14_curve,character-method
setMethod(
  f = "c14_curve",
  signature = c(name = "character"),
  definition = function(name) {
    curve_ok <- c("intcal20", "intcal13", "intcal09", "intcal04",
                  "marine20", "marine13", "marine09", "marine04",
                  "shcal20",  "shcal13",  "shcal04")
    name <- match.arg(name, choices = curve_ok, several.ok = TRUE)

    curves <- lapply(X = name, FUN = read_curve)
    names(curves) <- name
    curves
  }
)

#' @export
#' @rdname c14_curve
#' @aliases c14_curve,CalibratedAges-method
setMethod(
  f = "c14_curve",
  signature = c(name = "CalibratedAges"),
  definition = function(name) {
    methods::callGeneric(name@curves)
  }
)

read_curve <- function(x) {
  ## Read data
  curve_dir <- system.file("extdata", package = "ananke")
  curve_path <- file.path(curve_dir, paste0(x, ".14c"))
  curve <- utils::read.table(curve_path, header = FALSE, sep = ",",
                             dec = ".", strip.white = TRUE,
                             comment.char = "#")

  curve <- curve[, c(1, 2, 3)]
  colnames(curve) <- c("CALBP", "AGE", "ERROR")
  curve
}

cite_curve <- function(x) {
  curve <- c(
    bomb04nh1 = "Hua and Barbetti 2004",
    bomb04nh2 = "Hua and Barbetti 2004",
    bomb04nh3 = "Hua and Barbetti 2004",
    bomb04sh = "Hua and Barbetti 2004",
    bomb13nh1 = "Hua, Berbetti and Rakowski 2013",
    bomb13nh2 = "Hua, Berbetti and Rakowski 2013",
    bomb13nh3 = "Hua, Berbetti and Rakowski 2013",
    bomb13sh12 = "Hua, Berbetti and Rakowski 2013",
    bomb13sh3 = "Hua, Berbetti and Rakowski 2013",
    bomb21nh1 = "Hua et al. 2022",
    bomb21nh2 = "Hua et al. 2022",
    bomb21nh3 = "Hua et al. 2022",
    bomb21sh12 = "Hua et al. 2022",
    bomb21sh3 = "Hua et al. 2022",
    cariaco04 = "Hughen et al. 2004",
    intcal98 = "Stuiver et al. 1998",
    intcal04 = "Reimer et al. 2004",
    intcal09 = "Reimer et al. 2009",
    intcal13 = "Reimer et al. 2013",
    intcal20 = "Reimer et al. 2020",
    kueppers04 = "Kueppers et al. 2004",
    marine98 = "Stuiver, Reimer and Braziunas 1998",
    marine04 = "Hughen et al. 2004",
    marine09 = "Reimer et al. 2009",
    marine13 = "Reimer et al. 2013",
    marine20 = "Heaton et al. 2020",
    shcal04 = "McCormac et al. 2004",
    shcal13 = "Hogg et al. 2013",
    shcal20 = "Hogg et al. 2020"
  )

  curve[x]
}

# Approximate curve ============================================================
#' Interpolate 14C Calibration Curve
#'
#' @param name A [`character`] vector naming calibration curves.
#' @param out A [`numeric`] vector specifying where interpolation is to take
#'  place.
#' @param F14C A [`logical`] scalar: should estimated F14C values be used
#'  instead of radiocarbon ages?
#' @return
#'  A `list` of `list` with the following elements:
#'  \tabular{ll}{
#'  `mu`  \tab Interpolated values                    \cr
#'  `tau` \tab Interpolated errors                    \cr
#'  `max` \tab Maximum value of the calibration curve \cr
#'  `min` \tab Minimum value of the calibration curve \cr
#'  }
#' @keywords internal
#' @noRd
approx_curve <- function(name, out, F14C = FALSE) {
  ## Get data
  curve_data <- c14_curve(name)

  ## Interpolate
  lapply(
    X = curve_data,
    FUN = function(x, xout, F14C) {
      if (F14C) {
        x_f14c <- c14_f14c(x[, 2], x[, 3])
        x[, 2] <- x_f14c$value
        x[, 3] <- x_f14c$error
      }

      list(
        mu = stats::approx(x[, 1], x[, 2], xout = xout)$y,
        tau = stats::approx(x[, 1], x[, 3], xout = xout)$y,
        max = max(x[, 2]),
        min = min(x[, 2])
      )
    },
    xout = out,
    F14C = F14C
  )
}

# SUBSET
#' @include AllGenerics.R
NULL

# Extract ======================================================================
## [ ---------------------------------------------------------------------------
#' @export
#' @rdname subset
#' @aliases [,CalibratedAges-method
setMethod(
  f = "[",
  signature = c(x = "CalibratedAges"),
  function(x, i, j, k, drop = FALSE) {
    z <- x@.Data
    time <- x@.Time
    values <- x@values
    errors <- x@errors
    curves <- x@curves
    reservoir_offsets <- x@reservoir_offsets
    reservoir_errors <- x@reservoir_errors
    positions <- x@positions
    status <- x@status

    z <- z[i, j, k, drop = drop]
    if (!missing(i)) {
      if (is.character(i)) i <- match(i, dimnames(x)[1L])
      time <- time[i]
    }
    if (!missing(j)) {
      if (is.character(j)) j <- match(j, dimnames(x)[2L])
      values <- values[j]
      errors <- errors[j]
      curves <- curves[j]
      reservoir_offsets <- reservoir_offsets[j]
      reservoir_errors <- reservoir_errors[j]
      positions <- positions[j]
      status <- status[j]
    }

    if (isTRUE(drop)) return(z)
    methods::initialize(x, z, .Time = time, values = values, errors = errors,
                        curves = curves, reservoir_offsets = reservoir_offsets,
                        reservoir_errors = reservoir_errors,
                        positions = positions, status = status)
  }
)

# SAMPLE 14C
#' @include AllGenerics.R
NULL

#' @export
#' @rdname c14_sample
#' @aliases c14_sample,CalibratedAges-method
setMethod(
  f = "c14_sample",
  signature = "CalibratedAges",
  definition = function(object, n = 100, calendar = get_calendar()) {
    apply(
      X = object,
      MARGIN = 2,
      FUN = function(prob, size, x) {
        sample(x, size = size, replace = TRUE, prob = prob)
      },
      x = aion::time(object, calendar = calendar),
      size = n
    )
  }
)

# MUTATORS
#' @include AllGenerics.R
NULL

# Get ==========================================================================
#' @export
#' @method labels CalibratedAges
labels.CalibratedAges <- function(object, ...) {
  colnames(object) %||% paste0("X", NCOL(object))
}

#' @export
#' @rdname labels
#' @aliases labels,CalibratedAges-method
setMethod("labels", "CalibratedAges", labels.CalibratedAges)

# GENERIC METHODS
#' @include AllClasses.R
NULL

# Import S4 generics ===========================================================
#' @importMethodsFrom arkhe describe
#' @importMethodsFrom arkhe interval_credible
#' @importMethodsFrom arkhe interval_hdr
NULL

# Tools ========================================================================
## Mutators --------------------------------------------------------------------
#' Get or Set Parts of an Object
#'
#' Getters and setters to extract or replace parts of an object.
#' @param x An object from which to get or set element(s).
#' @param value A possible value for the element(s) of `x`.
#' @return
#'  An object of the same sort as `x` with the new values assigned.
# @example inst/examples/ex-mutator.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name mutators
#' @rdname mutators
#' @aliases get set
NULL

#' Find Labels from Object
#'
#' Find a suitable set of labels from an object for use in printing or plotting,
#' for example.
#' @param object An object from which to find labels.
#' @param ... Currently not used.
#' @return
#'  A [`character`] vector.
# @example inst/examples/ex-mutator.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name labels
#' @rdname labels
NULL

## Subset ----------------------------------------------------------------------
#' Extract or Replace Parts of an Object
#'
#' Operators acting on objects to extract or replace parts.
#' @param x An object from which to extract element(s) or in which to replace
#'  element(s).
#' @param i,j,k Indices specifying elements to extract or replace.
#' @param drop A [`logical`] scalar: should the result be coerced to
#'  the lowest possible dimension? This only works for extracting elements,
#'  not for the replacement.
# @param value A possible value for the element(s) of `x`.
# @param ... Currently not used.
#' @return
#'  A subsetted object.
# @example inst/examples/ex-subset.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name subset
#' @rdname subset
NULL

## Coerce ----------------------------------------------------------------------
#' Coerce to a Data Frame
#'
#' @inheritParams as.list
#' @return
#'  A [`data.frame`] with an extra `time` column.
#' @example inst/examples/ex-coerce.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name as.data.frame
#' @rdname as.data.frame
NULL

#' Coerce to a list
#'
#' @param x An object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param ... Currently not used.
#' @return
#'  A [`list`].
#' @example inst/examples/ex-coerce.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name as.list
#' @rdname as.list
NULL

# Radiocarbon ==================================================================
## Calibration curve -----------------------------------------------------------
#' 14C Calibration Curve
#'
#' @param name A [`character`] vector naming calibration curves (see details).
#' @param ... Currently not used.
#' @details
#'  The following calibration curves are available:
#'
#' \tabular{ll}{
#'  **Curve**    \tab **Reference** \cr
#  `bomb04nh1`  \tab Hua and Barbetti 2004 \cr
#  `bomb04nh2`  \tab Hua and Barbetti 2004 \cr
#  `bomb04nh3`  \tab Hua and Barbetti 2004 \cr
#  `bomb04sh`   \tab Hua and Barbetti 2004 \cr
#  `bomb13nh1`  \tab Hua, Berbetti and Rakowski 2013 \cr
#  `bomb13nh2`  \tab Hua, Berbetti and Rakowski 2013 \cr
#  `bomb13nh3`  \tab Hua, Berbetti and Rakowski 2013 \cr
#  `bomb13sh12` \tab Hua, Berbetti and Rakowski 2013 \cr
#  `bomb13sh3`  \tab Hua, Berbetti and Rakowski 2013 \cr
#  `bomb21nh1`  \tab Hua et al. 2022 \cr
#  `bomb21nh2`  \tab Hua et al. 2022 \cr
#  `bomb21nh3`  \tab Hua et al. 2022 \cr
#  `bomb21sh12` \tab Hua et al. 2022 \cr
#  `bomb21sh3`  \tab Hua et al. 2022 \cr
#  `cariaco04`  \tab Hughen et al. 2004 \cr
#  `intcal98`   \tab Stuiver et al. 1998 \cr
#'  `intcal04`   \tab Reimer et al. 2004 \cr
#'  `intcal09`   \tab Reimer et al. 2009 \cr
#'  `intcal13`   \tab Reimer et al. 2013 \cr
#'  `intcal20`   \tab Reimer et al. 2020 \cr
#  `kueppers04` \tab Kueppers et al. 2004 \cr
#d  `marine98`   \tab Stuiver, Reimer and Braziunas 1998 \cr
#'  `marine04`   \tab Hughen et al. 2004 \cr
#'  `marine09`   \tab Reimer et al. 2009 \cr
#'  `marine13`   \tab Reimer et al. 2013 \cr
#'  `marine20`   \tab Heaton et al. 2020 \cr
#'  `shcal04`    \tab McCormac et al. 2004 \cr
#'  `shcal13`    \tab Hogg et al. 2013 \cr
#'  `shcal20`    \tab Hogg et al. 2020 \cr
#' }
#'
#' @return
#'  A `list` of three-column [`data.frame`]:
#'  \tabular{ll}{
#'  `CALBP` \tab Calibrated age BP            \cr
#'  `AGE`   \tab Uncalibrated radiocarbon age \cr
#'  `ERROR` \tab Standard deviation           \cr
#'  }
#' @references
#'  Heaton, Timothy J, Peter Köhler, Martin Butzin, Edouard Bard, Ron W Reimer,
#'  William E N Austin, Christopher Bronk Ramsey, et al. (2020). Marine20 The
#'  Marine Radiocarbon Age Calibration Curve (0-55,000 Cal BP).
#'  *Radiocarbon*, 62(4): 779-820. \doi{10.1017/RDC.2020.68}.
#'
#'  Hogg, Alan G, Timothy J Heaton, Quan Hua, Jonathan G Palmer, Chris SM
#'  Turney, John Southon, Alex Bayliss, et al. (2020). SHCal20 Southern
#'  Hemisphere Calibration, 0-55,000 Years Cal BP. *Radiocarbon*, 62(4): 759-78.
#'  \doi{10.1017/RDC.2020.59}.
#'
#'  Hogg, Alan G, Quan Hua, Paul G Blackwell, Mu Niu, Caitlin E Buck, Thomas P
#'  Guilderson, Timothy J Heaton, et al. (2013). SHCal13 Southern Hemisphere
#'  Calibration, 0-50,000 Years Cal BP. *Radiocarbon*, 55(4): 1889-1903.
#'  \doi{10.2458/azu_js_rc.55.16783}.
#'
#'  Hua, Quan, and Mike Barbetti (2004). Review of Tropospheric Bomb 14C Data
#'  for Carbon Cycle Modeling and Age Calibration Purposes. *Radiocarbon*,
#'  46(3): 1273-1298. \doi{10.1017/S0033822200033142}.
#'
#'  Hua, Quan, Mike Barbetti, and Andrzej Z Rakowski (2013). Atmospheric
#'  Radiocarbon for the Period 1950-2010. *Radiocarbon*, 55(4): 2059‑2072.
#'  \doi{10.2458/azu_js_rc.v55i2.16177}.
#'
#'  Hua, Quan, Jocelyn C Turnbull, Guaciara M Santos, Andrzej Z Rakowski,
#'  Santiago Ancapichún, Ricardo De Pol-Holz, Samuel Hammer, et al. (2022).
#'  Atmospheric Radiocarbon for the Period 1950-2019. *Radiocarbon*,
#'  64(4): 723‑745. \doi{10.1017/RDC.2021.95}.
#'
#'  Hughen, K., S. Lehman, J. Southon, J. Overpeck, O. Marchal, C. Herring,
#'  and J. Turnbull (2004). 14C Activity and Global Carbon Cycle Changes over
#'  the Past 50,000 Years. *Science*, 303(5655): 202‑207.
#'  \doi{10.1126/science.1090300}.
#'
#'  Hughen, Konrad A, Mike G L Baillie, Edouard Bard, J Warren Beck, Chanda J H
#'  Bertrand, Paul G Blackwell, Caitlin E Buck, et al. (2004). Marine04 Marine
#'  Radiocarbon Age Calibration, 0-26 cal kyr BP. *Radiocarbon*,
#'  46(3): 1059‑1086. \doi{10.1017/S0033822200033002}.
#'
#'  Kueppers, Lara M., John Southon, Paul Baer, and John Harte (2004). Dead Wood
#'  Biomass and Turnover Time, Measured by Radiocarbon, along a Subalpine
#'  Elevation Gradient. *Oecologia*, 141(4): 641‑651.
#'  \doi{10.1007/s00442-004-1689-x}.
#'
#'  McCormac, F G, A G Hogg, P G Blackwell, C E Buck, T F G Higham, and P J
#'  Reimer (2004). Shcal04 Southern Hemisphere Calibration, 0-11.0 cal kyr BP.
#'  *Radiocarbon*, 46(3): 1087‑1092. \doi{10.1017/S0033822200033014}.
#'
#'  Reimer, P J, M G L Baillie, E Bard, A Bayliss, J W Beck, P G Blackwell,
#'  C Bronk Ramsey, et al. (2009). IntCal09 and Marine09 Radiocarbon Age
#'  Calibration Curves, 0-50,000 Years cal BP. *Radiocarbon*, 51(4): 1111‑1150.
#'  \doi{10.1017/S0033822200034202}.
#'
#'  Reimer, Paula J, William E N Austin, Edouard Bard, Alex Bayliss, Paul G
#'  Blackwell, Christopher Bronk Ramsey, Martin Butzin, et al. (2020).
#'  The IntCal20 Northern Hemisphere Radiocarbon Age Calibration Curve
#'  (0-55 cal kBP). *Radiocarbon*, 62(4): 725‑757. \doi{10.1017/RDC.2020.41}.
#'
#'  Reimer, Paula J, Mike G L Baillie, Edouard Bard, Alex Bayliss,
#'  J Warren Beck, Chanda J H Bertrand, Paul G Blackwell, et al. (2004).
#'  Intcal04 Terrestrial Radiocarbon Age Calibration, 0-26 cal kyr BP.
#'  *Radiocarbon*, 46(3): 1029‑1058. \doi{10.1017/S0033822200032999}.
#'
#'  Reimer, Paula J, Edouard Bard, Alex Bayliss, J Warren Beck, Paul G
#'  Blackwell, Christopher Bronk Ramsey, Caitlin E Buck, et al. (2013).
#'  IntCal13 and Marine13 Radiocarbon Age Calibration Curves 0-50,000
#'  Years cal BP. *Radiocarbon*, 55(4): 1869‑1887.
#'  \doi{10.2458/azu_js_rc.55.16947}.
#'
#'  Stuiver, Minze, Paula J. Reimer, Edouard Bard, J. Warren Beck, G. S. Burr,
#'  Konrad A. Hughen, Bernd Kromer, Gerry McCormac, Johannes van der Plicht, and
#'  Marco Spurk (1998). INTCAL98 Radiocarbon Age Calibration, 24,000-0 cal BP.
#'  *Radiocarbon*, 40(3): 1041‑1083. \doi{10.1017/S0033822200019123}.
#'
#'  Stuiver, Minze, Paula J. Reimer, and Thomas F. Braziunas. (1998).
#'  High-Precision Radiocarbon Age Calibration for Terrestrial and Marine
#'  Samples. *Radiocarbon*, 40(3): 1127‑1151. \doi{10.1017/S0033822200019172}.
#' @example inst/examples/ex-14c-curve.R
#' @author N. Frerebeau
#' @docType methods
#' @family radiocarbon tools
#' @aliases c14_curve-method
setGeneric(
  name = "c14_curve",
  def = function(name, ...) standardGeneric("c14_curve")
)

## Calibration -----------------------------------------------------------------
#' 14C Calibration
#'
#' Calibrates radiocarbon ages.
#' @param values A [`numeric`] vector giving the BP ages or F14C values to be
#'  calibrated (conventional ages).
#' @param errors A [`numeric`] vector giving the errors associated to the
#'  values to be calibrated.
#' @param curves A [`character`] vector specifying the calibration curve to be
#'  used. Different curves can be specified per sample.
#' @param names A [`character`] vector specifying the names of the samples (e.g.
#'  laboratory codes).
#' @param positions A [`numeric`] vector giving the position values (e.g.
#'  depths) for each age.
#' @param reservoir_offsets A [`numeric`] vector giving the offset values for
#'  any marine reservoir effect (defaults to 0; i.e. no offset).
#' @param reservoir_errors A [`numeric`] vector giving the offset value errors
#'  for any marine reservoir effect (defaults to 0; i.e. no offset).
#' @param from length-one [`numeric`] vector specifying the earliest data to
#'  calibrate for, in cal. BP years.
#' @param to A length-one [`numeric`] vector specifying the latest data to
#'  calibrate for, in cal. BP years.
#' @param resolution A length-one [`numeric`] vector specifying the temporal
#'  resolution (in years) of the calibration.
#' @param normalize A [`logical`] scalar: should the calibration be normalized?
#' @param F14C A [`logical`] scalar: should the calibration be carried out in
#'  F14C space? If `TRUE`, `values` must be expressed as F14C.
#' @param method A [`character`] string specifying the distribution assumed for
#'  the 14C ages. It must be one of "`student`" (the default) or "`normal`.
#'  Only used if `F14C` is `FALSE`.
#' @param dfs A [`character`] vector giving the degrees-of-freedom values for
#'  the student t-distribution associated with the calibration calculation.
#'  Only used if `method` is "`student`".
#' @param drop A [`logical`] scalar: should years with zero probability be
#'  discarded? If `TRUE` (the default), results in a narrower time range.
#' @param eps A length-one [`numeric`] value giving the cutoff below which
#'  calibration values will be removed.
#' @param verbose A [`logical`] scalar: should extra information be reported
#'  (e.g. warning message for dates out of calibration range)?
#' @param ... Currently not used.
#' @return
#'  A [`CalibratedAges-class`] object.
#' @references
#'  Bronk Ramsey, C. (2008). Radiocarbon Dating: Revolutions in Understanding.
#'  *Archaeometry*, 50:249-275. \doi{10.1111/j.1475-4754.2008.00394.x}.
#' @note
#'  Adapted from \pkg{Bchron} `BchronCalibrate()` by Andrew Parnell and
#'  \pkg{rcarbon} `calibrate()` by Andrew Bevan and Enrico Crema.
#' @example inst/examples/ex-14c-calibrate.R
#' @author N. Frerebeau
#' @docType methods
#' @family radiocarbon tools
#' @aliases c14_calibrate-method
setGeneric(
  name = "c14_calibrate",
  def = function(values, errors, ...) standardGeneric("c14_calibrate"),
  valueClass = "CalibratedAges"
)

#' Uncalibrate a Radiocarbon Date
#'
#' @param object A [`CalibratedAges-class`] object or a [`numeric`] vector of
#'  calibrated ages (in years BP).
#' @param n An [`integer`] specifying the number of random samples.
#' @param curves A [`character`] vector specifying the calibration curve to be
#'  used. Different curves can be specified.
#' @param rounding A [`character`] string specifying the rounding convention.
#'  It can be one of "`none`" (the default, no rounding) or "`stuiver`".
#'  Any unambiguous substring can be given.
#' @param ... Currently not used.
#' @example inst/examples/ex-14c-uncalibrate.R
#' @author N. Frerebeau
#' @docType methods
#' @family radiocarbon tools
#' @aliases c14_uncalibrate-method
setGeneric(
  name = "c14_uncalibrate",
  def = function(object, ...) standardGeneric("c14_uncalibrate")
)

## F14C ------------------------------------------------------------------------
#' F14C
#'
#' Converts F14C values to 14C ages.
#' @param values A [`numeric`] vector giving the radiocarbon ages or the F14C
#'  values.
#' @param errors A [`numeric`] vector giving the standard deviations.
#' @param tau A length-one [`numeric`] vector specifying the mean-life of
#'  radiocarbon (defaults to 14C half-life value as introduced by Libby 1952).
#' @param asymmetric A [`logical`] scalar: should asymmetric 14C errors be
#'  returned (van der Plicht & Hogg, 2006)?
#' @param rounding A [`character`] string specifying the rounding convention.
#'  It can be one of "`none`" (the default, no rounding) or "`stuiver`".
#'  Any unambiguous substring can be given.
#' @param ... Currently not used.
#' @details
#'  `f14c_c14()` calculates the conventional radiocarbon age:
#'
#'  \eqn{t = -\frac{1}{\lambda}\ln{F^{14}C} = -\tau\ln{F^{14}C}}
#'
#'  `c14_f14c()` calculates the inverse:
#'
#'  \eqn{F^{14}C = e^{-\lambda t} = e^{-\frac{t}{\tau}}}
#' @return
#'  A [`data.frame`].
#' @references
#'  Bronk Ramsey, C. (2008). Radiocarbon Dating: Revolutions in Understanding.
#'  *Archaeometry*, 50:249-275. \doi{10.1111/j.1475-4754.2008.00394.x}.
#'
#'  Stuiver, M., Polach, H. A. (1977). Discussion Reporting of 14C Data.
#'  *Radiocarbon*, 19(3): 355-363. \doi{10.1017/S0033822200003672}.
#'
#'  van der Plicht, J., Hogg, A. (2006). A Note on Reporting Radiocarbon.
#'  *Quaternary Geochronology*, 1(4): 237-240.
#'  \doi{10.1016/j.quageo.2006.07.001}.
#' @example inst/examples/ex-f14c.R
#' @author N. Frerebeau
#' @docType methods
#' @family radiocarbon tools
#' @name F14C
#' @rdname F14C
NULL

#' @rdname F14C
#' @aliases c14_f14c-method
setGeneric(
  name = "c14_f14c",
  def = function(values, errors, ...) standardGeneric("c14_f14c"),
  valueClass = "data.frame"
)

#' @rdname F14C
#' @aliases f14c_c14-method
setGeneric(
  name = "f14c_c14",
  def = function(values, errors, ...) standardGeneric("f14c_c14"),
  valueClass = "data.frame"
)

## Combine ---------------------------------------------------------------------
#' Combine 14C
#'
#' Combines radiocarbon dates.
#' @param values A [`numeric`] vector giving the BP ages to be calibrated.
#' @param errors A [`numeric`] vector giving the standard deviation of the ages
#'  to be calibrated.
#' @param groups A [`factor`] in the sense that `as.factor(groups)` defines the
#'  the groups to combine with. If `NULL` (the default), all dates are combined.
#'  `NA`s will be treated as isolated dates.
#' @param ... Currently not used.
#' @return
#'  A [`data.frame`] with the following columns:
#'  \tabular{ll}{
#'  `groups` \tab Group names                      \cr
#'  `ages`   \tab Combined 14C ages                \cr
#'  `errors` \tab Combined 14C standard deviations \cr
#'  `chi2`   \tab Chi-squared test statistic       \cr
#'  `p`      \tab Chi-squared test p-value         \cr
#'  }
#' @references
#'  Ward, G. K. and Wilson, S. R. (1978). Procedures for Comparing and Combining
#'  Radiocarbon Age Determinations: A Critique. *Archaeometry* 20(1): 19‑31.
#'  \doi{10.1111/j.1475-4754.1978.tb00208.x}.
#' @example inst/examples/ex-14c-combine.R
#' @author N. Frerebeau
#' @docType methods
#' @family radiocarbon tools
#' @aliases c14_combine-method
setGeneric(
  name = "c14_combine",
  def = function(values, errors, ...) standardGeneric("c14_combine"),
  valueClass = "data.frame"
)

## Plot ------------------------------------------------------------------------
#' Plot Calibrated Radiocarbon Ages
#'
#' @param x A [`CalibratedAges-class`] or [`CalibratedSPD-class`] object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param density A [`logical`] scalar: should density be drawn?
#' @param interval A [`character`] string specifying the intervals to be drawn.
#'  It must be one of "`hrd`" (the default), "`credible`" or "`none`".
#'  Any unambiguous substring can be given.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#'  Only used if `interval` is `TRUE`.
#' @param sort A [`logical`] scalar: should ages be arranged in chronological
#'  order?
#' @param decreasing A [`logical`] scalar: should the sort order be decreasing?
#' @param main A [`character`] string giving a main title for the plot.
#' @param sub A [`character`] string giving a subtitle for the plot.
#' @param ann A [`logical`] scalar: should the default annotation (title and x
#'  and y labels) appear on the plot?
#' @param axes A [`logical`] scalar: should axes be drawn on the plot?
#' @param frame.plot A [`logical`] scalar: should a box be drawn around the
#'  plot?
#' @param panel.first An an `expression` to be evaluated after the plot axes are
#'  set up but before any plotting takes place. This can be useful for drawing
#'  background grids.
#' @param panel.last An `expression` to be evaluated after plotting has taken
#'  place but before the axes, title and box are added.
#' @param col.density,col.interval A specification for the plotting colors.
#' @param ... Further parameters to be passed to [aion::plot()].
#' @return
#'  `plot()` is called it for its side-effects: it results in a graphic
#'  being displayed. Invisibly returns `x`.
#' @example inst/examples/ex-14c-plot.R
#' @author N. Frerebeau
#' @family radiocarbon tools
#' @docType methods
#' @name c14_plot
#' @rdname c14_plot
NULL

#' Ridgeline Plot of Calibrated Radiocarbon Ages
#'
#' @inheritParams c14_plot
#' @param fixed A [`logical`] scalar: should a fixed y scale be used?
#'  If `TRUE` (the default), ages are equally spaced along the y axis and
#'  arranged in chronological order. If `FALSE`, age `positions` are used
#'  (see [c14_calibrate()]).
#' @return
#'  `ridgelines()` is called it for its side-effects: it results in a graphic
#'  being displayed (invisibly returns `x`).
#' @example inst/examples/ex-14c-ridgelines.R
#' @author N. Frerebeau
#' @family radiocarbon tools
#' @docType methods
#' @aliases ridgelines-method
setGeneric(
  name = "ridgelines",
  def = function(x, ...) standardGeneric("ridgelines")
)

#' Plot a Radiocarbon Event Count Ensemble
#'
#' @param x An [`RECE-class`] object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param ... Further parameters to be passed to [graphics::image()].
#' @return
#'  `plot()` is called it for its side-effects: it results in a graphic being
#'  displayed (invisibly returns `x`).
#' @references
#'  Carleton, W. C. (2021). Evaluating Bayesian Radiocarbon‐dated Event Count
#'  (REC) Models for the Study of Long‐term Human and Environmental Processes.
#'  *Journal of Quaternary Science*, 36(1): 110‑23. \doi{10.1002/jqs.3256}.
#' @author N. Frerebeau
#' @family radiocarbon tools
#' @docType methods
#' @name rec_plot
#' @rdname rec_plot
NULL

## Sample ----------------------------------------------------------------------
#' Sample Calibrated Ages
#'
#' @param object A [`CalibratedAges-class`] object.
#' @param n An [`integer`] specifying the number of random samples.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). Defaults to [aion::CE()]. If `NULL`,
#'  *rata die* are returned.
#' @param ... Currently not used.
#' @example inst/examples/ex-14c-sample.R
#' @return An `numeric` matrix.
#' @author N. Frerebeau
#' @family radiocarbon tools
#' @docType methods
#' @aliases c14_sample-method
setGeneric(
  name = "c14_sample",
  def = function(object, ...) standardGeneric("c14_sample")
)

## RECE ------------------------------------------------------------------------
#' Radiocarbon Event Count
#'
#' @param object A [`CalibratedAges-class`] object.
#' @param from length-one [`numeric`] vector specifying the earliest data to
#'  calibrate for (in cal BP years).
#' @param to A length-one [`numeric`] vector specifying the latest data to
#'  calibrate for (in cal BP years).
#' @param by A length-one [`numeric`] vector specifying the temporal
#'  resolution (in years) of the calibration.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). Defaults to [aion::CE()]. If `NULL`,
#'  *rata die* are returned.
#' @param n An [`integer`] specifying the number of item to choose randomly.
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @param ... Currently not used.
#' @return
#'  An [`RECE-class`] object.
#' @note
#'  This function is currently *experimental*.
#' @references
#'  Carleton, W. C. (2021). Evaluating Bayesian Radiocarbon‐dated Event Count
#'  (REC) Models for the Study of Long‐term Human and Environmental Processes.
#'  *Journal of Quaternary Science*, 36(1): 110‑23. \doi{10.1002/jqs.3256}.
#' @author N. Frerebeau
#' @family radiocarbon tools
#' @docType methods
#' @aliases c14_ensemble-method
setGeneric(
  name = "c14_ensemble",
  def = function(object, ...) standardGeneric("c14_ensemble"),
  valueClass = "RECE"
)

## SPD -------------------------------------------------------------------------
#' Summed Probability Distributions
#'
#' Computes summed probability distributions (SPD) of radiocarbon dates.
#' @param object A [`CalibratedAges-class`] object.
#' @param normalize_date A [`logical`] scalar: should the total probability mass
#'  of the calibrated dates be normalised (to sum to unity within the time-span
#'  of analysis)?
#' @param normalize_spd A [`logical`] scalar: should the total probability mass
#'  of the SPD be normalised (to sum to unity)?
#' @param ... Currently not used.
#' @details
#'  Summed probability distributions (SPD) are not statistically valid
#'  estimators of the calendar age of a potential future sample. They should not
#'  be used in any dates-as-data approach to provide a population proxy.
#' @return
#'  A [`CalibratedSPD-class`] object.
#' @example inst/examples/ex-14c-spd.R
#' @author N. Frerebeau
#' @docType methods
#' @family radiocarbon tools
#' @aliases c14_spd-method
setGeneric(
  name = "c14_spd",
  def = function(object, ...) standardGeneric("c14_spd"),
  valueClass = "CalibratedSPD"
)

# Interval =====================================================================
## HDR -------------------------------------------------------------------------
#' Highest Density Regions
#'
#' @param x A [`CalibratedAges-class`] object.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @param ... Currently not used.
#' @return
#'  A [`CalibratedIntervals-class`] object.
#' @references
#'  Hyndman, R. J. (1996). Computing and graphing highest density regions.
#'  *American Statistician*, 50: 120-126. \doi{10.2307/2684423}.
#' @example inst/examples/ex-14c-hdr.R
#' @seealso [stats::density()], [arkhe::interval_hdr()]
#' @author N. Frerebeau
#' @family statistics
#' @docType methods
#' @rdname interval_hdr
#' @name interval_hdr
NULL

## Credible --------------------------------------------------------------------
#' Bayesian Credible Interval
#'
#' @param x A [`CalibratedAges-class`] object.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @param n An [`integer`] specifying the number of random samples.
#' @param ... Currently not used.
#' @return
#'  A [`CalibratedIntervals-class`] object.
#' @example inst/examples/ex-14c-credible.R
#' @seealso [arkhe::interval_credible()]
#' @author N. Frerebeau
#' @family statistics
#' @docType methods
#' @rdname interval_credible
#' @name interval_credible
NULL

# Proxy Records ================================================================
#' Layer-Counted Proxy Records Uncertainties
#'
#' Represents layer-counted proxy records as sequences of probability
#' distributions on absolute, error-free time axes.
#' @param positions A positive [`numeric`] vector giving the positions (e.g.
#'  depths) at which proxy values and calendar ages were measured. In the case
#'  of layers of non-zero thickness, this should be the middle value of the
#'  slice. It must be in decreasing order (i.e. in chronological order).
#' @param proxy_values A [`numeric`] vector giving the proxy values.
#' @param proxy_errors A [`numeric`] vector giving the proxy uncertainties.
#' @param proxy_step A length-one [`numeric`] vector specifying the step size
#'  (in units of `proxy_values`) at which proxy records densities are to be
#'  estimated.
#' @param time_values A [`numeric`] vector giving the calendar ages (in years).
#' @param time_errors A [`numeric`] vector giving the calendar age uncertainties
#'  (in years).
#' @param calendar An [`aion::TimeScale-class`] object specifying the calendar
#'  of `time` (see [aion::calendar()]).
#' @param from A length-one [`numeric`] vector specifying the starting value of
#'  the temporal sequence at which densities are to be estimated (in years).
#' @param to A length-one [`numeric`] vector specifying the end value of the
#'  temporal sequence at which densities are to be estimated (in cal BP years).
#' @param by A length-one [`numeric`] vector specifying the increment of
#'  the temporal sequence at which densities are to be estimated (in years).
#' @param n An [`integer`] specifying the number of item to choose randomly.
#' @param verbose A [`logical`] scalar: should extra information be reported?
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @param ... Currently not used.
#' @return
#'  A [`ProxyRecord-class`] object.
#' @note
#'  This function is currently *experimental*.
#' @references
#'  Boers, N., Goswami, B. & Ghil, M. (2017). A Complete Representation of
#'  Uncertainties in Layer-Counted Paleoclimatic Archives. *Climate of the
#'  Past*, 13(9): 1169-1180. \doi{10.5194/cp-13-1169-2017}.
#' @example inst/examples/ex-proxy.R
#' @author N. Frerebeau
#' @family proxy tools
#' @docType methods
#' @aliases proxy_ensemble-method
setGeneric(
  name = "proxy_ensemble",
  def = function(positions, ...) standardGeneric("proxy_ensemble"),
  valueClass = "ProxyRecord"
)

#' Plot Layer-Counted Proxy Records Uncertainties
#'
#' @param x A [`ProxyRecord-class`] object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param iqr A [`logical`] scalar: should the mean and IQR be displayed?
#' @param xlab,ylab A [`character`] string giving a label for the x and y axis.
#' @param col A list of colors such as that generated by [grDevices::hcl.colors()].
#' @param col.mean,col.iqr A specification for the line colors. Only used if
#'  `iqr` is `TRUE`.
#' @param lty.mean,lty.iqr A specification for the line types. Only used if
#'  `iqr` is `TRUE`.
#' @param lwd.mean,lwd.iqr A specification for the line widths. Only used if
#'  `iqr` is `TRUE`.
#' @param ... Further parameters to be passed to [graphics::image()].
#' @return
#'  `plot()` is called it for its side-effects: it results in a graphic
#'  being displayed. Invisibly returns `x`.
#' @example inst/examples/ex-proxy.R
#' @author N. Frerebeau
#' @docType methods
#' @family proxy tools
#' @name proxy_plot
#' @rdname proxy_plot
NULL

# Isotopes =====================================================================
#' Geological Model Age from Lead Isotope Analysis
#'
#' Compute geological model age (T) and U/Pb (mu) and Th/U (kappa) ratios from
#' lead isotopic measurements.
#' @param x A [`numeric`] vector of 206Pb/204Pb ratios. If `y` and `z` are
#'  missing, must be a [`list`] (or a [`data.frame`]) with `numeric` components
#'  (columns) `x`, `y` and `z`.
#' @param y A [`numeric`] vector of 207Pb/204Pb ratios.  If missing, an attempt
#'  is made to interpret `x` in a suitable way.
#' @param z A [`numeric`] vector of 208Pb/204Pb ratios.  If missing, an attempt
#'  is made to interpret `x` in a suitable way.
#' @param t0 A [`numeric`] value giving the time of the second stage of the
#'  reference model.
#' @param x_star A [`numeric`] value giving the 206Pb/204Pb ratio at
#'  \eqn{t = 0}.
#' @param y_star A [`numeric`] value giving the 207Pb/204Pb ratio at
#'  \eqn{t = 0}.
#' @param z_star A [`numeric`] value giving the 208Pb/204Pb ratio at
#'  \eqn{t = 0}.
#' @param mu A [`numeric`] value giving the 238U/204Pb ratio of the
#'  reference model.
#' @param kappa A [`numeric`] value giving the 232Th/238U ratio of the
#'  reference model.
#' @param th232 A [`numeric`] value giving the decay constants of 232Th.
#' @param u238 A [`numeric`] value giving the decay constants of 238U.
#' @param u235 A [`numeric`] value giving the decay constants of 235U.
#' @param u238_235 A [`numeric`] value giving the actual 238U/235U ratio.
#' @param tolerance A [`numeric`] value specifying the tolerance (stopping
#'  criteria for the Newton–Raphson method).
#' @param stop An [`integer`] giving the stopping rule (i.e. maximum number of
#'  iterations) to avoid infinite loop.
#' @param ... Currently not used.
#' @note
#'  Reference values from Albarede & Juteau (1984).
#' @return
#'  A four columns [`data.frame`]:
#'  \describe{
#'   \item{`age`}{Geological model age (in Ma).}
#'   \item{`mu`}{238U/204Pb ratio.}
#'   \item{`kappa`}{232Th/238U ratio.}
#'   \item{`residual`}{Newton loop residual.}
#'  }
#' @references
#'  Albarède, F., Desaulty, A.-M. & Blichert-Toft, J. (2012). A Geological
#'  Perspective on the Use of Pb Isotopes in Archaeometry. *Archaeometry*, 54:
#'  853-867. \doi{10.1111/j.1475-4754.2011.00653.x}.
#'
#'  Albarède, F. & Juteau, M. (1984). Unscrambling the Lead Model Ages.
#'  *Geochimica et Cosmochimica Acta*, 48(1): 207-12.
#'  \doi{10.1016/0016-7037(84)90364-8}.
#'
#'  Allègre, C. (2005). *Géologie isotopique*. Belin sup. Paris: Belin.
#' @example inst/examples/ex-lia.R
#' @author N. Frerebeau, F. Albarede (original Matlab code)
#' @docType methods
#' @family isotope analysis
#' @export
setGeneric(
  name = "pb_age",
  def = function(x, y, z, ...) standardGeneric("pb_age"),
  valueClass = "data.frame"
)

# Statistics ===================================================================
#' Quantiles of a Density Estimate
#'
#' @param x A [`CalibratedAges-class`] object.
#' @param probs A [`numeric`] vector of probabilities with values in
#'  \eqn{[0,1]}.
#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
#'  computation proceeds?
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param ... Currently not used.
#' @return
#'  A `numeric` [`matrix`] containing the quantiles.
#' @example inst/examples/ex-14c-statistics.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name quantile
#' @rdname quantile
NULL

#' Mean
#'
#' @param x A [`CalibratedAges-class`] object.
#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
#'  computation proceeds?
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param ... Currently not used.
#' @return A [`numeric`] vector.
#' @example inst/examples/ex-14c-statistics.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name mean
#' @rdname mean
NULL

#' Median
#'
#' @param x A [`CalibratedAges-class`] object.
#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
#'  computation proceeds?
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param ... Currently not used.
#' @return A [`numeric`] vector.
#' @example inst/examples/ex-14c-statistics.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name median
#' @rdname median
NULL

# Summary ======================================================================
#' Object Summaries
#'
#' @param object A [`CalibratedAges-class`] object.
#' @param digits An [`integer`] (or `NULL`) used for number formatting (see
#'  [base::format()]).
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param ... Currently not used.
#' @return
#'  An object of class [`table`][base::table].
#' @example inst/examples/ex-summary.R
#' @author N. Frerebeau
#' @docType methods
#' @family summary
#' @name summary
#' @rdname summary
NULL

#' Data Description
#'
#' @param x A [`CalibratedAges-class`] object.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]).
#' @param ... Further parameters to be passed to [cat()].
#' @return
#'  `describe()` is called for its side-effects. Invisibly returns `x`.
#' @references
#'  Millard, A. R. (2014). Conventions for Reporting Radiocarbon Determinations.
#'  *Radiocarbon*, 56(2): 555-559. \doi{10.2458/56.17455}.
#' @example inst/examples/ex-describe.R
#' @author N. Frerebeau
#' @family summary
#' @docType methods
#' @rdname describe
#' @name describe
NULL

1		# HELPERS
2
3		## https://michaelchirico.github.io/potools/articles/developers.html
4		tr_ <- function(...) {
5	8x	enc2utf8(gettext(paste0(...), domain = "R-ananke"))
6		}

1		# DESCRIBE
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname describe
7		#' @aliases describe,CalibratedAges-method
8		setMethod(
9		f = "describe",
10		signature = signature(x = "CalibratedAges"),
11		definition = function(x, calendar = get_calendar(),
12		level = 0.954, ...) {
13		## Get data
14	1x	lab <- labels(x)
15	1x	val <- x@values
16	1x	err <- x@errors
17	1x	crv <- x@curves
18	1x	reservoir_off <- x@reservoir_offsets
19	1x	reservoir_err <- x@reservoir_errors
20	1x	F14C <- x@F14C
21
22		## Laboratory code
23	1x	if (F14C) {
24	!	txt_uncal <- "Sample %s contains %.0f +/- %.0f F14C,"
25		} else {
26	1x	txt_uncal <- "Sample %s is dated to %.0f +/- %.0f BP,"
27		}
28	1x	msg_uncal <- sprintf(txt_uncal, lab, val, err)
29
30		## Calibration results
31	1x	hdr <- as.list(interval_hdr(x, level = level), calendar = calendar)
32	1x	msg_cal <- lapply(
33	1x	X = hdr,
34	1x	FUN = function(x, calendar, level) {
35	!	if (is.null(x)) return("but is out of the calibration range of")
36	2x	p <- if (NROW(x) > 1) sprintf(" (%.1f%%)", x$p * 100) else ""
37	2x	msg_hdr <- sprintf("[%.0f,%.0f]%s", x$start, x$end, p)
38	2x	txt_cal <- "calibrated to %s %s (%.1f%% HPD interval) with"
39	2x	sprintf(txt_cal, paste0(msg_hdr, collapse = " or "), calendar, level)
40		},
41	1x	calendar = calendar@label,
42	1x	level = level * 100
43		)
44
45		## Calibration curve
46	1x	txt_marine <- "marine reservoir offset: %.0f +/- %.0f; "
47	1x	msg_marine <- sprintf(txt_marine, reservoir_off, reservoir_err)
48	1x	msg_marine[reservoir_off == 0] <- ""
49	1x	txt_curve <- "%s (%s%s)."
50	1x	msg_curve <- sprintf(txt_curve, crv, msg_marine, cite_curve(crv))
51
52		## Text
53	1x	msg <- paste(msg_uncal, msg_cal, msg_curve, sep = " ")
54
55		## Software
56	1x	txt_soft <- "Calibration was computed with R %s.%s (R Core Team %s) and package ananke %s (Frerebeau %s)."
57	1x	date_soft <- utils::packageDate("ananke")
58	1x	date_soft <- if (is.na(date_soft)) Sys.Date() else date_soft
59	1x	msg_soft <- sprintf(txt_soft, R.version$major, R.version$minor,
60	1x	R.version$year, utils::packageVersion("ananke"),
61	1x	format(date_soft, format = "%Y"))
62
63		## Split string
64	1x	fill <- list(...)$fill
65	1x	if (!is.null(fill) && isTRUE(fill)) {
66	!	msg <- lapply(
67	!	X = msg,
68	!	FUN = function(x, width) {
69	!	split <- lapply(
70	!	X = seq(from = 1, to = nchar(x), by = width),
71	!	FUN = function(i) substr(x, i, i + width - 1)
72		)
73	!	paste0(split, collapse = "\n")
74		},
75	!	width = getOption("width")
76		)
77		}
78
79	1x	cat(unlist(msg), msg_soft, sep = "\n\n", ...)
80
81	1x	invisible(x)
82		}
83		)

1		# PLOT
2		#' @include AllGenerics.R
3		NULL
4
5		## 14C =========================================================================
6		#' @export
7		#' @method plot CalibratedAges
8		plot.CalibratedAges <- function(x, calendar = get_calendar(),
9		interval = c("hdr", "credible"),
10		level = 0.954, sort = TRUE, decreasing = FALSE,
11		col.density = "grey", col.interval = "#77AADD",
12		main = NULL, sub = NULL,
13		ann = graphics::par("ann"),
14		axes = TRUE, frame.plot = TRUE,
15		panel.first = NULL, panel.last = NULL, ...) {
16		## Check
17	6x	c14_validate(x)
18	6x	interval <- match.arg(interval, several.ok = FALSE)
19
20		## Sort
21	6x	if (isTRUE(sort)) {
22	4x	mid <- order(mean(x, calendar = NULL), decreasing = !decreasing)
23	4x	x <- x[, mid, , drop = FALSE]
24		}
25
26		## Graphical parameters
27	6x	fill.density <- grDevices::adjustcolor(col.density, alpha.f = 0.5)
28	6x	fill.interval <- grDevices::adjustcolor(col.interval, alpha.f = 0.5)
29
30		## Plot
31	6x	panel_density <- function(x, y, ...) {
32	18x	force(interval)
33	18x	force(level)
34	18x	force(fill.interval)
35
36	18x	tick_bottom <- min(y, na.rm = TRUE)
37	18x	tick_height <- tick_bottom + graphics::par("tcl") * graphics::strheight("M") * -1
38
39		## Keep only density > 0
40	18x	d0 <- which(y > tick_bottom)
41	18x	x <- x[d0]
42	18x	y <- y[d0]
43
44		## Draw density
45	18x	graphics::polygon(
46	18x	x = c(x, rev(x)),
47	18x	y = c(y, rep(tick_bottom, length(y))),
48	18x	border = NA,
49		...
50		)
51
52		## Add interval
53	18x	if (level > 0) {
54	18x	y0 <- arkhe::scale_range(y)
55	18x	int <- switch(
56	18x	interval,
57	18x	hdr = arkhe::interval_hdr(x, y0, level = level),
58	18x	credible = {
59	6x	spl <- sample(x, size = length(x), replace = TRUE, prob = y0)
60	6x	arkhe::interval_credible(spl, level = level)
61		}
62		)
63
64	18x	for (i in seq_len(nrow(int))) {
65	58x	debut <- int[i, "start"]
66	58x	fin <- int[i, "end"]
67	54x	if (debut < fin) is_in_int <- which(x >= debut & x <= fin)
68	4x	else is_in_int <- which(x <= debut & x >= fin)
69	58x	xi <- x[is_in_int]
70	58x	yi <- y[is_in_int]
71	58x	graphics::polygon(
72	58x	x = c(xi[1], xi, xi[length(xi)]),
73	58x	y = c(tick_bottom, yi, tick_bottom),
74	58x	border = NA, col = fill.interval
75		)
76	58x	graphics::segments(
77	58x	x0 = c(debut, debut, fin),
78	58x	x1 = c(debut, fin, fin),
79	58x	y0 = c(tick_bottom, tick_bottom, tick_bottom),
80	58x	y1 = c(tick_height, tick_bottom, tick_height),
81	58x	lend = 1
82		)
83		}
84		}
85
86	18x	graphics::lines(x, y, col = "black")
87		}
88
89		## Method for TimeSeries
90	6x	methods::callNextMethod(
91	6x	x, facet = "multiple",
92	6x	calendar = calendar,
93	6x	panel = panel_density,
94	6x	main = main, sub = sub, ann = ann, axes = axes,
95	6x	frame.plot = frame.plot,
96	6x	panel.first = panel.first,
97	6x	panel.last = panel.last,
98	6x	col = fill.density,
99		...
100		)
101
102	6x	invisible(x)
103		}
104
105		#' @export
106		#' @rdname c14_plot
107		#' @aliases plot,CalibratedAges,missing-method
108		setMethod("plot", c(x = "CalibratedAges", y = "missing"), plot.CalibratedAges)
109
110		#' @export
111		#' @rdname ridgelines
112		#' @aliases ridgelines,CalibratedAges-method
113		setMethod(
114		f = "ridgelines",
115		signature = c(x = "CalibratedAges"),
116		definition = function(x, calendar = get_calendar(),
117		interval = c("hdr", "credible"),
118		level = 0.954, fixed = TRUE, decreasing = FALSE,
119		col.density = "grey", col.interval = "#77AADD", ...) {
120		## Get data
121	2x	lab <- labels(x)
122
123		## Y scale
124	2x	if (isTRUE(fixed)) {
125	2x	mid <- mean(x, calendar = NULL)
126	2x	dy <- order(order(mid, decreasing = decreasing))
127		} else {
128	!	dy <- order(order(x@positions, decreasing = !decreasing))
129		}
130	2x	for (j in seq_len(ncol(x))) x[, j, ] <- arkhe::scale_range(x[, j, ]) + dy[j]
131
132		## Permute
133	2x	x <- aion::flip(x)
134
135		## Plot
136	2x	plot(x, calendar = calendar, interval = interval, level = level,
137	2x	sort = FALSE, col.density = col.density,
138	2x	col.interval = col.interval,
139	2x	axes = FALSE, frame.plot = FALSE, yaxt = "n", ...)
140
141		## Construct Axis
142	2x	aion::year_axis(side = 1, format = TRUE, calendar = calendar,
143	2x	xpd = NA)
144	2x	graphics::axis(side = 2, at = dy, labels = lab, las = 2,
145	2x	lty = 0)
146		}
147		)
148
149		## SPD =========================================================================
150		#' @export
151		#' @method plot CalibratedSPD
152		plot.CalibratedSPD <- function(x, calendar = get_calendar(),
153		main = NULL, sub = NULL,
154		ann = graphics::par("ann"),
155		axes = TRUE, frame.plot = TRUE,
156		panel.first = NULL, panel.last = NULL, ...) {
157		## Graphical parameters
158	!	n <- NCOL(x)
159	!	col <- list(...)$col %\|\|% c("grey")
160	!	if (length(col) != n) col <- rep(col, length.out = n)
161	!	col <- grDevices::adjustcolor(col, alpha.f = 0.5)
162
163		## Plot
164	!	panel_density <- function(x, y, ...) {
165	!	graphics::polygon(
166	!	x = c(x, rev(x)),
167	!	y = c(y, rep(0, length(y))),
168	!	border = NA,
169		...
170		)
171	!	graphics::lines(x, y, col = "black")
172		}
173
174		## Method for TimeSeries
175	!	methods::callNextMethod(
176	!	x, facet = "multiple",
177	!	calendar = calendar,
178	!	panel = panel_density,
179	!	main = main, sub = sub, ann = ann, axes = axes,
180	!	frame.plot = frame.plot,
181	!	panel.first = panel.first,
182	!	panel.last = panel.last,
183	!	col = col,
184		...
185		)
186
187	!	invisible(x)
188		}
189
190		#' @export
191		#' @rdname c14_plot
192		#' @aliases plot,CalibratedSPD,missing-method
193		setMethod("plot", c(x = "CalibratedSPD", y = "missing"), plot.CalibratedSPD)
194
195		## RECE ========================================================================
196		#' @export
197		#' @method plot RECE
198		plot.RECE <- function(x, calendar = get_calendar(), ...) {
199		## Binary array
200	!	bin <- array(FALSE, dim = c(nrow(x), max(x), ncol(x)))
201	!	for (j in seq_len(ncol(x))) {
202	!	z <- x[, j, , drop = TRUE]
203	!	for (i in seq_along(z)) {
204	!	bin[i, z[i], j] <- z[i] > 0
205		}
206		}
207	!	bin <- apply(X = bin, MARGIN = c(1, 2), FUN = sum)
208	!	bin[bin == 0] <- NA
209
210		## Add annotation
211	!	years <- aion::time(x, calendar = NULL)
212
213		## Plot
214	!	graphics::image(
215	!	x = years,
216	!	y = seq_len(max(x)),
217	!	z = log(bin),
218	!	xlab = format(calendar),
219	!	ylab = "Count",
220	!	xaxt = "n",
221	!	yaxt = "n",
222		...
223		)
224
225		## Construct axes
226	!	aion::year_axis(side = 1, format = TRUE, calendar = calendar,
227	!	current_calendar = NULL)
228	!	graphics::axis(side = 2, at = seq_len(max(x)), las = 1)
229
230	!	invisible(x)
231		}
232
233		#' @export
234		#' @rdname rec_plot
235		#' @aliases plot,RECE,missing-method
236		setMethod("plot", c(x = "RECE", y = "missing"), plot.RECE)
237
238
239		## Proxy =======================================================================
240		#' @export
241		#' @method plot ProxyRecord
242		plot.ProxyRecord <- function(x, calendar = get_calendar(),
243		iqr = TRUE,
244		xlab = NULL, ylab = NULL,
245		col = grDevices::hcl.colors(12, "YlOrRd", rev = TRUE),
246		col.mean = "black", col.iqr = col.mean,
247		lty.mean = 1, lty.iqr = 3,
248		lwd.mean = 2, lwd.iqr = lwd.mean, ...) {
249		## Get data
250	!	years <- aion::time(x, calendar = NULL)
251	!	z <- apply(
252	!	X = x@density,
253	!	MARGIN = 1,
254	!	FUN = function(d) (d - min(d)) / max(d - min(d))
255		)
256	!	z[z == 0] <- NA
257
258		## Plot
259	!	graphics::image(
260	!	x = years,
261	!	y = x@proxy,
262	!	z = t(z),
263	!	col = col,
264	!	xaxt = "n",
265	!	yaxt = "n",
266	!	xlab = xlab %\|\|% format(calendar),
267	!	ylab = ylab %\|\|% "Proxy",
268		...
269		)
270
271		## Construct axes
272	!	aion::year_axis(side = 1, format = TRUE, calendar = calendar,
273	!	current_calendar = NULL)
274	!	graphics::axis(side = 2, las = 1)
275
276	!	if (isTRUE(iqr)) {
277	!	m <- mean(x)
278	!	graphics::lines(x = years, y = m, col = col.mean,
279	!	lty = lty.mean, lwd = lwd.mean)
280
281	!	q <- quantile(x, probs = c(0.25, 0.75))
282	!	graphics::lines(x = years, y = q[, 1], col = col.iqr,
283	!	lty = lty.iqr, lwd = lwd.iqr)
284	!	graphics::lines(x = years, y = q[, 2], col = col.iqr,
285	!	lty = lty.iqr, lwd = lwd.iqr)
286		}
287
288	!	invisible(x)
289		}
290
291		#' @export
292		#' @rdname proxy_plot
293		#' @aliases plot,ProxyRecord,missing-method
294		setMethod("plot", c(x = "ProxyRecord", y = "missing"), plot.ProxyRecord)

1		# 14C SPD
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname c14_spd
7		#' @aliases c14_spd,CalibratedAges-method
8		setMethod(
9		f = "c14_spd",
10		signature = "CalibratedAges",
11		definition = function(object, normalize_date = FALSE, normalize_spd = FALSE) {
12		## Check
13	1x	c14_validate(object)
14
15	1x	dens <- t(object[, , 1, drop = TRUE])
16	!	if (normalize_date) dens <- dens / rowSums(dens, na.rm = TRUE)
17	1x	spd <- colSums(dens, na.rm = TRUE)
18	!	if (normalize_spd) spd <- spd / sum(spd, na.rm = TRUE)
19
20	1x	time_series <- aion::series(
21	1x	object = spd,
22	1x	time = aion::time(object, calendar = NULL)
23		)
24	1x	.CalibratedSPD(time_series)
25		}
26		)

1		# COERCION
2		#' @include AllGenerics.R
3		NULL
4
5		# To list ======================================================================
6		#' @export
7		#' @method as.list CalibratedIntervals
8		as.list.CalibratedIntervals <- function(x, ..., calendar = get_calendar()) {
9	4x	z <- as.data.frame(x, calendar = calendar)
10	4x	f <- factor(z$label, levels = unique(z$label))
11	4x	z$label <- NULL
12	4x	split(x = z, f = f)
13		}
14
15		#' @export
16		#' @rdname as.list
17		#' @aliases as.list,CalibratedIntervals-method
18		setMethod("as.list", "CalibratedIntervals", as.list.CalibratedIntervals)
19
20		# To data.frame ================================================================
21		#' @export
22		#' @method as.data.frame CalibratedAges
23		as.data.frame.CalibratedAges <- function(x, ...,
24		calendar = get_calendar()) {
25	!	dens <- x[, , 1, drop = TRUE]
26	!	z <- data.frame(aion::time(x, calendar = calendar), dens)
27	!	colnames(z) <- c("time", colnames(x) %\|\|% paste0("X", seq_len(NCOL(x))))
28	!	z
29		}
30
31		#' @export
32		#' @rdname as.data.frame
33		#' @aliases as.data.frame,CalibratedAges-method
34		setMethod("as.data.frame", "CalibratedAges", as.data.frame.CalibratedAges)
35
36		#' @export
37		#' @method as.data.frame CalibratedIntervals
38		as.data.frame.CalibratedIntervals <- function(x, ...,
39		calendar = get_calendar()) {
40		## Build a data frame
41	4x	data.frame(
42	4x	label = labels(x),
43	4x	start = start(x, calendar = calendar),
44	4x	end = end(x, calendar = calendar),
45	4x	p = x@p
46		)
47		}
48
49		#' @export
50		#' @rdname as.data.frame
51		#' @aliases as.data.frame,CalibratedIntervals-method
52		setMethod("as.data.frame", "CalibratedIntervals", as.data.frame.CalibratedIntervals)
53
54		#' @export
55		#' @method as.data.frame RECE
56		as.data.frame.RECE <- function(x, ..., calendar = get_calendar()) {
57	!	dens <- x[, , 1, drop = TRUE]
58	!	z <- data.frame(aion::time(x, calendar = calendar), dens)
59	!	colnames(z) <- c("time", colnames(x) %\|\|% paste0("X", seq_len(NCOL(x))))
60	!	z
61		}
62
63		#' @export
64		#' @rdname as.data.frame
65		#' @aliases as.data.frame,RECE-method
66		setMethod("as.data.frame", "RECE", as.data.frame.RECE)
67
68		#' @export
69		#' @method as.data.frame ProxyRecord
70		as.data.frame.ProxyRecord <- function(x, ...,
71		calendar = get_calendar()) {
72	!	dens <- x[, , 1, drop = TRUE]
73	!	z <- data.frame(aion::time(x, calendar = calendar), dens)
74	!	colnames(z) <- c("time", colnames(x) %\|\|% paste0("X", seq_len(NCOL(x))))
75	!	z
76		}
77
78		#' @export
79		#' @rdname as.data.frame
80		#' @aliases as.data.frame,ProxyRecord-method
81		setMethod("as.data.frame", "ProxyRecord", as.data.frame.ProxyRecord)

1		# 14C UNCALIBRATION
2		#' @include AllGenerics.R
3		NULL
4
5		# Uncalibrate ==================================================================
6		#' @export
7		#' @rdname c14_uncalibrate
8		#' @aliases c14_uncalibrate,numeric-method
9		setMethod(
10		f = "c14_uncalibrate",
11		signature = c(object = "numeric"),
12		definition = function(object, curves = "intcal20") {
13		## Validation
14	1x	n <- length(object)
15	!	if (length(curves) != n) curves <- rep(curves, n)
16	1x	uncal <- rep(NA_real_, n)
17
18		## Calibration curve
19	1x	curve_data <- c14_curve(tolower(curves))
20
21	1x	for (i in seq_len(n)) {
22	1x	z <- stats::approx(
23	1x	x = curve_data[[i]][, 1],
24	1x	y = curve_data[[i]][, 2],
25	1x	xout = object[i],
26	1x	rule = 1
27		)
28	1x	uncal[i] <- z$y
29		}
30	1x	uncal
31		}
32		)
33
34		#' @export
35		#' @rdname c14_uncalibrate
36		#' @aliases c14_uncalibrate,CalibratedAges-method
37		setMethod(
38		f = "c14_uncalibrate",
39		signature = c(object = "CalibratedAges"),
40		definition = function(object, n = 10000, rounding = getOption("ananke.round"),
41		...) {
42
43	!	method <- list(...)$method %\|\|% c("L-BFGS-B")
44
45		## Function to be optimized
46	!	opt_fun <- function(param, curve, samples, ...) {
47	!	cal <- c14_calibrate(
48	!	values = round(param[1]),
49	!	errors = round(param[2]),
50	!	curves = curve
51		)
52	!	s1 <- samples
53	!	s2 <- c14_sample(cal, n = n, calendar = BP())
54
55	!	dens1 <- stats::density(s1, from = min(c(s1, s2)), to = max(c(s1, s2)))$y
56	!	dens2 <- stats::density(s2, from = min(c(s1, s2)), to = max(c(s1, s2)))$y
57	!	dens1 <- dens1 / sum(dens1)
58	!	dens2 <- dens2 / sum(dens2)
59	!	int <- dens1 * log(dens1 / dens2)
60	!	int <- int[dens1 > 0]
61	!	int <- int[!is.infinite(int)]
62	!	sum(int)
63		}
64
65	!	n <- NCOL(object)
66	!	curves <- object@curves
67
68		## Initial values for the parameters to be optimized over
69	!	spl <- c14_sample(object, n = n, calendar = BP())
70	!	med <- apply(X = spl, MARGIN = 2, FUN = stats::median)
71	!	init_mean <- c14_uncalibrate(med, curves = curves)
72	!	init_sd <- apply(X = spl, MARGIN = 2, FUN = stats::sd)
73
74	!	opt_mean <- opt_sd <- rep(NA_real_, n)
75	!	for (i in seq_len(n)) {
76	!	cal_spl <- spl[, i]
77
78	!	opt_run <- stats::optim(
79	!	par = c(init_mean[[i]], init_sd[[i]]),
80	!	fn = opt_fun,
81	!	method = method,
82	!	curve = curves[[i]],
83	!	samples = cal_spl,
84	!	lower = c(-Inf, 1),
85		...
86		)
87	!	opt_mean[[i]] <- opt_run$par[1]
88	!	opt_sd[[i]] <- opt_run$par[2]
89		}
90
91		## Rounding
92	!	if (identical(rounding, "stuiver")) {
93	!	opt_mean <- round_values_stuiver(opt_mean)
94	!	opt_sd <- round_errors_stuiver(opt_sd)
95		}
96
97	!	data.frame(mean = opt_mean, sd = opt_sd)
98		}
99		)

1		# REC
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname c14_ensemble
7		#' @aliases c14_ensemble,CalibratedAges-method
8		setMethod(
9		f = "c14_ensemble",
10		signature = "CalibratedAges",
11		definition = function(object, from = NULL, to = NULL, by = 10, n = 100,
12		calendar = BP(), progress = getOption("ananke.progress")) {
13		## Check
14	!	c14_validate(object)
15
16		## Get data
17	!	rd <- aion::time(object, calendar = NULL)
18	!	if (is.null(from)) from <- aion::start(object, calendar = calendar)
19	!	if (is.null(to)) to <- aion::end(object, calendar = calendar)
20	!	grid_years <- seq(from = from, to = to, by = by * sign(to - from))
21	!	grid_rd <- aion::fixed(grid_years, calendar = calendar)
22
23		## Align 14C date densities onto the grid
24	!	c14_dens <- object[, , 1, drop = TRUE]
25	!	c14_aligned <- apply(
26	!	X = c14_dens,
27	!	MARGIN = 2,
28	!	FUN = function(y, x, grid) {
29	!	fun <- stats::approxfun(x = x, y = y)
30	!	fun(grid)
31		},
32	!	x = rd,
33	!	grid = grid_rd
34		)
35	!	c14_aligned[is.na(c14_aligned)] <- 0
36
37		## Build matrix to store the RECE
38	!	count <- matrix(data = 0, nrow = n, ncol = length(grid_rd))
39	!	colnames(count) <- grid_rd
40
41	!	progress_bar <- progress
42	!	if (progress_bar) pbar <- utils::txtProgressBar(max = n, style = 3)
43
44	!	n_seq <- seq_len(n)
45	!	for (i in n_seq) {
46		## Sample
47	!	spl <- apply(
48	!	X = c14_aligned,
49	!	MARGIN = 2,
50	!	FUN = function(x, grid) {
51	!	if (sum(x) == 0) return(NA)
52	!	sample(grid, size = 1, prob = x)
53		},
54	!	grid = grid_rd
55		)
56		## Count
57	!	tbl <- unclass(table(spl)) # Named integer vector
58	!	count[i, names(tbl)] <- tbl
59	!	if (progress_bar) utils::setTxtProgressBar(pbar, i)
60		}
61	!	count[is.na(count)] <- 0
62	!	if (progress_bar) close(pbar)
63
64		## Return an RECE object
65	!	ts <- aion::series(
66	!	object = t(count),
67	!	time = grid_rd
68		# names = colnames(object)
69		)
70	!	.RECE(ts)
71		}
72		)

1		# AGE-DEPTH MODELING
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname proxy_ensemble
7		#' @aliases proxy_ensemble,numeric-method
8		setMethod(
9		f = "proxy_ensemble",
10		signature = c("numeric"),
11		definition = function(positions, proxy_values, proxy_errors, proxy_step,
12		time_values, time_errors, calendar,
13		from = NULL, to = NULL, by = NULL, n = 30,
14		progress = getOption("ananke.progress"),
15		verbose = getOption("ananke.verbose")) {
16		## Validation
17	!	k <- length(positions)
18	!	if (length(proxy_errors) == 1) proxy_errors <- rep(proxy_errors, k)
19	!	arkhe::assert_positive(positions)
20	!	arkhe::assert_decreasing(positions)
21	!	arkhe::assert_length(proxy_values, k)
22	!	arkhe::assert_length(proxy_errors, k)
23	!	arkhe::assert_length(proxy_step, 1)
24	!	arkhe::assert_length(time_values, k)
25	!	arkhe::assert_length(time_errors, k)
26
27		## Missing values
28	!	if (is.null(from)) from <- time_values[[1L]]
29	!	if (is.null(to)) to <- time_values[[k]]
30	!	if (is.null(by)) {
31	!	grid <- getOption("ananke.grid")
32	!	by <- ((to - from) / (grid - 1))
33		}
34	!	if (from > to && by > 0) by <- by * -1
35
36		## Build a matrix to contain the p(t\|zi) densities
37		## Rows will refer to depth
38		## Columns will refer to the time density
39	!	if (verbose) cat(tr_("Computing p(t\|zi) densities..."), sep = "\n")
40	!	t_grid <- seq(from = from, to = to, by = by)
41
42	!	t_z <- .mapply(
43	!	FUN = function(min, max, x) {
44	!	stats::dunif(x = x, min = min, max = max)
45		},
46	!	dots = list(
47	!	min = time_values - 1 * time_errors,
48	!	max = time_values + 1 * time_errors
49		),
50	!	MoreArgs = list(x = t_grid)
51		)
52	!	t_z <- do.call(rbind, t_z)
53	!	t_z <- t(t_z)
54
55		## Build a matrix to contain the p(x\|zi) densities
56		## Rows will refer to depth
57		## Columns will refer to the proxy density
58	!	if (verbose) cat(tr_("Computing p(x\|zi) densities..."), sep = "\n")
59	!	d <- 2 * max(proxy_errors)
60	!	x_range <- range(c(range(proxy_values) - d, range(proxy_values) + d))
61	!	x_grid <- seq(from = x_range[[1L]], to = x_range[[2L]], by = proxy_step)
62
63	!	x_z <- .mapply(
64	!	FUN = function(mean, sd, x) {
65	!	stats::dnorm(x = x, mean = mean, sd = sd)
66		},
67	!	dots = list(
68	!	mean = proxy_values,
69	!	sd = proxy_errors
70		),
71	!	MoreArgs = list(x = x_grid)
72		)
73	!	x_z <- do.call(rbind, x_z)
74	!	x_z <- t(x_z)
75
76		## Estimate the weighted average density function
77		## (for a given proxy measurement at a given time)
78		## Eq. 4 of Boers et al. 2017
79	!	if (verbose) cat(tr_("Computing p(x\|t) densities..."), sep = "\n")
80
81	!	z <- length(positions)
82	!	ri <- vapply(
83	!	X = 2:(z - 1),
84	!	FUN = function(x, y) (y[x + 1] - y[x - 1]) / 2,
85	!	FUN.VALUE = numeric(1),
86	!	y = positions
87		)
88	!	r <- abs(c(positions[2] - positions[1], ri, positions[z] - positions[z - 1]))
89	!	r_t <- matrix(data = r, nrow = nrow(t_z), ncol = ncol(t_z), byrow = TRUE)
90	!	r_x <- matrix(data = r, nrow = nrow(x_z), ncol = ncol(x_z), byrow = TRUE)
91
92	!	x_t <- t(tcrossprod(r_x * x_z, t_z)) / rowSums(r_t * t_z)
93
94		## Create an ensemble of potential proxy records
95	!	if (verbose) cat(tr_("Sampling proxy records..."), sep = "\n")
96	!	Y <- matrix(data = 0, nrow = length(t_grid), ncol = n)
97
98	!	if (progress) pb <- utils::txtProgressBar(min = 0, max = n, style = 3)
99
100	!	n_seq <- seq_len(n)
101	!	for (i in n_seq) {
102		## Sample
103	!	Y[, i] <- apply(
104	!	X = x_t,
105	!	MARGIN = 1,
106	!	FUN = function(x, g) sample(g, size = 1, prob = x),
107	!	g = x_grid
108		)
109	!	if (progress) utils::setTxtProgressBar(pb, i)
110		}
111
112	!	if (progress) close(pb)
113
114	!	ts <- aion::series(object = Y, time = t_grid, calendar = calendar)
115	!	.ProxyRecord(ts, density = x_t, proxy = x_grid)
116		}
117		)

1		# STATISTICS
2		#' @include AllGenerics.R
3		NULL
4
5		# Quantile =====================================================================
6		#' Quantiles of a Density Estimate
7		#'
8		#' @param x A [`numeric`] vector giving the n coordinates of the points where
9		#' the density is estimated.
10		#' @param y A `numeric` [`matrix`] of density estimates (each column is a
11		#' density estimate).
12		#' @param probs A [`numeric`] vector of probabilities with values in
13		#' \eqn{[0,1]}.
14		#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
15		#' computation proceeds?
16		#' @param ... Currently not used.
17		#' @return
18		#' A `numeric` [`matrix`] containing the quantiles.
19		#' @keywords internal
20		#' @noRd
21		quantile_density <- function(x, y, probs = seq(0, 1, 0.25), na.rm = FALSE, ...) {
22	3x	eps <- 100 * .Machine$double.eps
23	3x	if (anyNA(probs) \|\| any(probs < -eps \| probs > 1 + eps))
24	!	stop(sprintf("%s outside [0,1]", sQuote("probs")))
25
26	3x	q <- apply(
27	3x	X = y,
28	3x	MARGIN = 2,
29	3x	FUN = function(y, x, probs, na.rm) {
30	6x	np <- length(probs)
31	6x	qs <- rep(NA_real_, np)
32	6x	if (na.rm) {
33	!	i <- !is.na(x) & !is.na(y)
34	!	x <- x[i]
35	!	y <- y[i]
36		}
37	6x	if (np > 0) {
38	6x	nn <- length(x)
39	6x	Fx <- cumsum(y * c(0, diff(x)))
40	6x	Fx <- Fx / Fx[nn]
41	6x	for (j in seq_len(np)) {
42	18x	ii <- min(which(Fx >= probs[j]))
43	18x	if (!is.na(ii) && ii >= 1 && ii <= nn) qs[j] <- x[ii]
44		}
45	6x	qs
46		}
47		},
48	3x	x = x,
49	3x	probs = probs,
50	3x	na.rm = na.rm
51		)
52
53	3x	if (!is.null(dim(q))) {
54	2x	q <- t(q)
55	2x	colnames(q) <- paste0(round(probs * 100, digits = 0), "%")
56		}
57	3x	q
58		}
59
60		#' @export
61		#' @method quantile CalibratedAges
62		quantile.CalibratedAges <- function(x, probs = seq(0, 1, 0.25), na.rm = FALSE,
63		..., calendar = get_calendar()) {
64	2x	quantile_density(x = aion::time(x, calendar = calendar), y = x,
65	2x	probs = probs, na.rm = na.rm, ...)
66		}
67
68		#' @export
69		#' @rdname quantile
70		#' @aliases quantile,CalibratedAges,missing-method
71		setMethod("quantile", c(x = "CalibratedAges"), quantile.CalibratedAges)
72
73		#' @export
74		#' @method quantile ProxyRecord
75		quantile.ProxyRecord <- function(x, probs = seq(0, 1, 0.25), na.rm = FALSE, ...) {
76	!	quantile_density(x = x@proxy, y = t(x@density),
77	!	probs = probs, na.rm = na.rm, ...)
78		}
79
80		#' @export
81		#' @rdname quantile
82		#' @aliases quantile,ProxyRecord,missing-method
83		setMethod("quantile", c(x = "ProxyRecord"), quantile.ProxyRecord)
84
85		# Median =======================================================================
86		#' @export
87		#' @method median CalibratedAges
88		median.CalibratedAges <- function(x, na.rm = FALSE, ...,
89		calendar = get_calendar()) {
90	1x	quantile_density(x = aion::time(x, calendar = calendar), y = x,
91	1x	probs = 0.5, na.rm = na.rm, ...)
92		}
93
94		#' @export
95		#' @rdname median
96		#' @aliases median,CalibratedAges,missing-method
97		setMethod("median", c(x = "CalibratedAges"), median.CalibratedAges)
98
99		# Mean =========================================================================
100		#' Mean of a Density Estimate
101		#'
102		#' @param x A [`numeric`] vector giving the n coordinates of the points where
103		#' the density is estimated.
104		#' @param y A `numeric` [`matrix`] of density estimates (each column is a
105		#' density estimate).
106		#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
107		#' computation proceeds?
108		#' @param ... Currently not used.
109		#' @return A [`numeric`] vector.
110		#' @keywords internal
111		#' @noRd
112		mean_density <- function(x, y, na.rm = FALSE, ...) {
113	8x	apply(
114	8x	X = y,
115	8x	MARGIN = 2,
116	8x	FUN = function(w, x, na.rm) {
117	22x	if (na.rm) {
118	!	i <- !is.na(w) & !is.na(x)
119	!	x <- x[i]
120	!	w <- w[i]
121		}
122	22x	stats::weighted.mean(x = x, w = w)
123		},
124	8x	x = x,
125	8x	na.rm = na.rm
126		)
127		}
128
129		#' @export
130		#' @method mean CalibratedAges
131		mean.CalibratedAges <- function(x, na.rm = FALSE, ...,
132		calendar = get_calendar()) {
133	8x	mean_density(x = aion::time(x, calendar = calendar), y = x, na.rm = na.rm, ...)
134		}
135
136		#' @export
137		#' @rdname mean
138		#' @aliases mean,CalibratedAges,missing-method
139		setMethod("mean", c(x = "CalibratedAges"), mean.CalibratedAges)
140
141		#' @export
142		#' @method mean ProxyRecord
143		mean.ProxyRecord <- function(x, na.rm = FALSE, ...) {
144	!	mean_density(x = x@proxy, y = t(x@density), na.rm = na.rm, ...)
145		}
146
147		#' @export
148		#' @rdname mean
149		#' @aliases mean,ProxyRecord,missing-method
150		setMethod("mean", c(x = "ProxyRecord"), mean.ProxyRecord)
151
152		# Summary ======================================================================
153		#' @export
154		#' @method summary CalibratedAges
155		summary.CalibratedAges <- function(object, ..., digits = NULL,
156		calendar = get_calendar()) {
157	1x	mid <- matrix(mean(object, calendar = get_calendar()), ncol = 1)
158	1x	qq <- quantile(object, probs = c(0.25, 0.50, 0.75), calendar = get_calendar())
159
160	1x	sms <- t(cbind(qq[, 1L:2L, drop = FALSE], mid, qq[, 3L, drop = FALSE]))
161	1x	sms <- format(sms, digits = digits)
162
163	1x	lbs <- format(c(tr_("1st Qu."), tr_("Median"), tr_("Mean"), tr_("3rd Qu.")))
164	1x	lbs <- matrix(lbs, nrow = nrow(sms), ncol = ncol(sms))
165	1x	z <- paste0(lbs, ":", sms, " ")
166	1x	nm <- format(colnames(sms), justify = "centre", width = max(nchar(z)))
167
168	1x	dim(z) <- dim(sms)
169	1x	dimnames(z) <- list(rep("", nrow(sms)), nm)
170	1x	attr(z, "class") <- c("table")
171
172	1x	z
173		}
174
175		#' @export
176		#' @rdname summary
177		#' @aliases summary,CalibratedAges-method
178		setMethod("summary", c(object = "CalibratedAges"), summary.CalibratedAges)

1		# COMBINE 14C
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname c14_combine
7		#' @aliases c14_combine,numeric,numeric-method
8		setMethod(
9		f = "c14_combine",
10		signature = c(values = "numeric", errors = "numeric"),
11		definition = function(values, errors, groups = NULL) {
12		## Validation
13	4x	n <- length(values)
14	1x	if (is.null(groups)) groups <- "X"
15	3x	if (length(groups) == 1) groups <- rep(groups, n)
16
17		## Empty groups must be treated as NA
18	4x	groups[which(!nzchar(groups))] <- NA
19	4x	groups <- factor(x = groups, levels = unique(groups))
20
21	4x	arkhe::assert_missing(values)
22	4x	arkhe::assert_missing(errors)
23	4x	arkhe::assert_length(errors, n)
24	4x	arkhe::assert_length(groups, n)
25
26		## Groups with only one date must be treated as NA
27	4x	counts <- table(groups)
28	4x	one <- groups %in% names(counts)[counts == 1]
29
30		# NA group will be removed
31		# We need to keep isolated dates
32	4x	k <- one \| is.na(groups)
33	4x	solo <- NULL
34	4x	if (any(k)) {
35	3x	solo <- data.frame(
36	3x	groups = as.character(groups[k]),
37	3x	values = values[k],
38	3x	errors = errors[k],
39	3x	chi2 = NA_real_,
40	3x	p = NA_real_
41		)
42		}
43
44	4x	combined <- NULL
45	4x	if (!all(k)) {
46	2x	groups[k] <- NA
47	2x	groups <- droplevels(groups)
48
49		## split() removes NA group
50	2x	values <- split(values, f = groups)
51	2x	errors <- split(errors, f = groups)
52	2x	cmbn <- mapply(
53	2x	FUN = combine,
54	2x	values = values,
55	2x	errors = errors,
56	2x	SIMPLIFY = FALSE
57		)
58	2x	combined <- data.frame(names(cmbn), do.call(rbind, cmbn))
59	2x	colnames(combined) <- c("groups", "values", "errors", "chi2", "p")
60		}
61
62	4x	final <- rbind(solo, combined, make.row.names = FALSE)
63	4x	final
64		}
65		)
66
67		combine <- function(values, errors) {
68		## On calcule la moyenne pondérée
69	2x	w <- 1 / errors^2 # Facteur de pondération
70	2x	moy <- stats::weighted.mean(x = values, w = w)
71
72		## On calcule l'incertitude associée à la moyenne pondérée
73	2x	err <- sum(1 / errors^2)^(-1 / 2)
74
75		## On calcule la statistique du test
76	2x	chi2 <- sum(((values - moy) / errors)^2)
77
78		## On calcule la valeur-p
79	2x	p <- 1 - stats::pchisq(chi2, df = length(values))
80
81		## On stocke les résultats
82	2x	c(moy, err, chi2, p)
83		}

1		# F14C <> BP14C
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname F14C
7		#' @aliases c14_f14c,numeric,numeric-method
8		setMethod(
9		f = "c14_f14c",
10		signature = c(values = "numeric", errors = "numeric"),
11		definition = function(values, errors, tau = 8033) {
12	1x	values <- exp(values / -tau)
13	1x	sigma <- values * errors / tau
14	1x	data.frame(value = values, error = sigma)
15		}
16		)
17
18		#' @export
19		#' @rdname F14C
20		#' @aliases f14c_c14,numeric,numeric-method
21		setMethod(
22		f = "f14c_c14",
23		signature = c(values = "numeric", errors = "numeric"),
24		definition = function(values, errors, tau = 8033, asymmetric = FALSE,
25		rounding = getOption("ananke.round")) {
26		## Validation
27	5x	choices <- c("none", "stuiver")
28	5x	rounding <- match.arg(rounding, choices = choices, several.ok = FALSE)
29
30	5x	z <- values
31
32		## van der Plicht and Hogg 2006, p. 239
33	5x	below_2sigma <- z < 2 * errors
34	5x	values[below_2sigma] <- values[below_2sigma] + 2 * errors[below_2sigma]
35
36	5x	below_zero <- z < 0
37	5x	values[below_zero] <- 2 * errors[below_zero]
38
39		## van der Plicht and Hogg 2006, eq. 6
40		## Bronk Ramsey 2008, p. 260
41	5x	ages <- -tau * log(values)
42	5x	sigma <- tau * errors / values
43
44	5x	sigma_plus <- sigma_minus <- sigma
45	5x	if (asymmetric) {
46	3x	sigma_plus <- - tau * log(values - errors) - ages
47	3x	sigma_minus <- ages + tau * log(values + errors)
48		}
49
50	5x	sigma_plus[below_zero \| below_2sigma] <- NA_real_
51	5x	sigma_minus[below_zero \| below_2sigma] <- NA_real_
52
53		## The reported age can be no older than the radiocarbon age of the
54		## fraction modern of the sample plus it's error
55		# limit <- -tau * log(values + errors)
56		# ages[ages > limit] <- limit[ages > limit]
57
58		## Rounding
59	5x	if (identical(rounding, "stuiver")) {
60	2x	ages <- round_values_stuiver(ages)
61	2x	sigma_plus <- round_errors_stuiver(sigma_plus)
62	2x	sigma_minus <- round_errors_stuiver(sigma_minus)
63		}
64
65	5x	data.frame(age = ages, plus = sigma_plus, minus = sigma_minus)
66		}
67		)
68
69		## Stuiver & Polach (1977), p. 362
70		round_any <- function(x, accuracy, f = round) {
71	31x	f(x / accuracy) * accuracy
72		}
73		round_values_stuiver <- function(x) {
74	4x	y <- x[!is.na(x)]
75	4x	y[y < 1000] <- round_any(y[y < 1000], 5)
76	4x	y[1000 <= y & y <= 9999] <- round_any(y[1000 <= y & y < 9999], 10)
77	4x	y[10000 <= y & y <= 20000] <- round_any(y[10000 <= y & y <= 20000], 50)
78	4x	y[y > 20000] <- round_any(y[y > 20000], 100)
79	4x	x[!is.na(x)] <- y
80	4x	x
81		}
82		round_errors_stuiver <- function(x) {
83	5x	y <- x[!is.na(x)]
84	5x	y[y < 50] <- round_any(y[y < 50], 5)
85	5x	y[50 <= y & y <= 1000] <- round_any(y[50 <= y & y <= 1000], 10)
86	5x	y[y > 1000] <- round_any(y[y > 1000], 100)
87	5x	x[!is.na(x)] <- y
88	5x	x
89		}

1		# 14C CALIBRATION
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname c14_calibrate
7		#' @aliases c14_calibrate,numeric,numeric-method
8		setMethod(
9		f = "c14_calibrate",
10		signature = c(values = "numeric", errors = "numeric"),
11		definition = function(values, errors, curves = "intcal20",
12		names = NULL, positions = NULL,
13		reservoir_offsets = 0, reservoir_errors = 0,
14		from = 55000, to = 0, resolution = 1,
15		normalize = TRUE, F14C = FALSE,
16		method = c("student", "normal"),
17		dfs = 100, drop = TRUE, eps = 1e-06,
18		verbose = getOption("ananke.verbose")) {
19		## Validation
20	12x	method <- match.arg(method)
21	12x	n <- length(values)
22	9x	if (is.null(names)) names <- paste0("X", seq_len(n))
23	12x	if (is.null(positions)) positions <- values
24	12x	if (length(curves) == 1)
25	12x	curves <- rep(curves, n)
26	12x	if (length(reservoir_offsets) == 1)
27	12x	reservoir_offsets <- rep(reservoir_offsets, n)
28	12x	if (length(reservoir_errors) == 1)
29	12x	reservoir_errors <- rep(reservoir_errors, n)
30	12x	if (length(dfs) == 1)
31	12x	dfs <- rep(dfs, n)
32
33	12x	arkhe::assert_missing(values)
34	12x	arkhe::assert_missing(errors)
35	12x	arkhe::assert_unique(names)
36	12x	arkhe::assert_length(errors, n)
37	12x	arkhe::assert_length(names, n)
38	12x	arkhe::assert_length(curves, n)
39	12x	arkhe::assert_length(reservoir_offsets, n)
40	12x	arkhe::assert_length(reservoir_errors, n)
41	12x	arkhe::assert_length(dfs, n)
42
43		## Calibration time range
44	12x	cal_range <- seq(from = from, to = to, by = -resolution)
45
46		## Calibration curve
47	12x	curves <- tolower(curves)
48	12x	curve_range <- approx_curve(unique(curves), out = cal_range, F14C = F14C)
49
50		## Marine reservoir offset
51	11x	values <- values - reservoir_offsets
52	11x	errors <- sqrt(errors^2 + reservoir_errors^2)
53
54		## Calibrate
55	11x	calibrate_fun <- if (F14C) calibrate_F14C else calibrate_BP14C
56	11x	dens <- vector(mode = "list", length = n)
57	11x	status <- integer(n)
58	11x	for (i in seq_len(n)) {
59	17x	d <- calibrate_fun(
60	17x	x = values[i],
61	17x	error = errors[i],
62	17x	mu = curve_range[[curves[i]]]$mu,
63	17x	tau = curve_range[[curves[i]]]$tau,
64	17x	df = dfs[i],
65	17x	method = method
66		)
67
68		## Check
69	17x	if (anyNA(d)) {
70	!	d[is.na(d)] <- 0
71	!	msg <- tr_("Consider changing the calibration time range to a narrower interval.")
72	!	if (verbose) message(msg)
73		}
74
75	17x	max_cal <- curve_range[[curves[i]]]$max
76	17x	min_cal <- curve_range[[curves[i]]]$min
77
78	17x	if (values[[i]] >= max_cal \|\| values[[i]] <= min_cal) {
79		## L'age à calibrer est hors de l'étendue de la courbe de calibration
80	!	if (verbose) warning(print_out(names[[i]], maybe = FALSE), call. = FALSE)
81	3x	status[i] <- 1L
82	14x	} else if (d[1] > eps \|\| d[length(d)] > eps) {
83		## L'age calibré est en partie hors de l'étendue de la courbe
84	!	if (verbose) warning(print_out(names[[i]], maybe = TRUE), call. = FALSE)
85	3x	status[i] <- 2L
86		}
87
88	17x	d[d < eps] <- 0
89	17x	dens[[i]] <- d
90		}
91
92		## Build matrix
93	11x	dens <- do.call(rbind, dens)
94
95		## Normalize
96	!	if (F14C && !normalize) normalize <- TRUE
97	11x	if (normalize) {
98	11x	dens <- dens / rowSums(dens, na.rm = TRUE)
99	11x	dens[dens < eps] <- 0
100	11x	dens <- dens / rowSums(dens, na.rm = TRUE)
101		}
102
103		## Drop
104	11x	if (drop) {
105	11x	keep_zero <- colSums(dens, na.rm = TRUE) > 0
106	11x	keep_from <- which.max(keep_zero) # First TRUE
107	11x	keep_to <- length(keep_zero) - which.max(rev(keep_zero)) + 1 # Last TRUE
108	11x	keep <- seq(from = keep_from, to = keep_to, by = 1)
109	11x	dens <- dens[, keep, drop = FALSE]
110	11x	cal_range <- cal_range[keep]
111		}
112
113	11x	time_series <- aion::series(
114	11x	object = t(dens),
115	11x	time = cal_range,
116	11x	calendar = BP(),
117	11x	names = names
118		)
119	11x	.CalibratedAges(
120	11x	time_series,
121	11x	values = values,
122	11x	errors = errors,
123	11x	curves = curves,
124	11x	reservoir_offsets = reservoir_offsets,
125	11x	reservoir_errors = reservoir_errors,
126	11x	F14C = F14C,
127	11x	positions = positions,
128	11x	status = status
129		)
130		}
131		)
132
133		calibrate_BP14C <- function(x, error, mu, tau, df = 100, method = "student") {
134	17x	tau <- error^2 + tau^2
135	17x	dens <- switch(
136	17x	method,
137	17x	student = stats::dt(x = (x - mu) / sqrt(tau), df = df),
138	17x	normal = stats::dnorm(x, mean = mu, sd = sqrt(tau)),
139	17x	stop(sprintf("Unknown method %s.", dQuote(method)))
140		)
141	17x	dens
142		}
143		calibrate_F14C <- function(x, error, mu, tau, ...) {
144	!	p1 <- (x - mu)^2
145	!	p2 <- 2 * (error^2 + tau^2)
146	!	p3 <- sqrt(error^2 + tau^2)
147	!	dens <- exp(-p1 / p2) / p3
148	!	dens
149		}

1		# 14C CALIBRATION CHECK
2		#' @include AllGenerics.R
3		NULL
4
5		#' Check Calibrated Radiocarbon Dates
6		#'
7		#' @param object A [`CalibratedAges-class`] object.
8		#' @param verbose A [`logical`] scalar: should extra information be reported?
9		#' @return
10		#' `c14_validate()` is called it for its side-effects: it prints
11		#' [warning messages][warnings()] if calibrated ages are (partially) out of
12		#' calibration range. Invisibly returns `x`.
13		#' @example inst/examples/ex-14c-calibrate.R
14		#' @author N. Frerebeau
15		#' @keywords internal
16		c14_validate <- function(object, verbose = getOption("ananke.verbose")) {
17	15x	status <- object@status
18	15x	lab <- labels(object)
19
20	15x	if (isTRUE(verbose)) {
21	4x	if (any(status == 1L)) {
22	2x	is_out <- which(status == 1L)
23	2x	warn <- print_out(lab[is_out], maybe = FALSE)
24	2x	for (w in warn) warning(w, call. = FALSE)
25		}
26	4x	if (any(status == 2L)) {
27	2x	is_out <- which(status == 2L)
28	2x	warn <- print_out(lab[is_out], maybe = TRUE)
29	2x	for (w in warn) warning(w, call. = FALSE)
30		}
31		}
32
33	15x	invisible(object)
34		}
35
36		print_out <- function(label, maybe = FALSE) {
37	4x	if (maybe) {
38	2x	status <- tr_("Date %s may extent out of calibration range.")
39		} else {
40	2x	status <- tr_("Date %s is out of calibration range.")
41		}
42	4x	sprintf(status, dQuote(label))
43		}

1		# GEOLOGICAL MODEL AGE
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname pb_age
7		#' @aliases pb_age,numeric,numeric,numeric-method
8		setMethod(
9		f = "pb_age",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, t0 = 3.8,
12		x_star = 18.75, y_star = 15.63, z_star = 38.86,
13		mu = 9.66, kappa = 3.90, th232 = 0.049475,
14		u238 = 0.155125, u235 = 0.98485, u238_235 = 137.79,
15		tolerance = sqrt(.Machine$double.eps),
16		stop = 100) {
17		## Validation
18	2x	n <- length(x)
19	2x	arkhe::assert_length(y, n)
20	2x	arkhe::assert_length(z, n)
21
22	2x	x0 <- x_star - mu * (exp(u238 * t0) - 1) # 206Pb/204Pb at t0
23	2x	y0 <- y_star - mu / u238_235 * (exp(u235 * t0) - 1) # 207Pb/204Pb at t0
24	2x	z0 <- z_star - mu * kappa * (exp(th232 * t0) - 1) # 208Pb/204Pb at t0
25
26	2x	t_model <- vector(mode = "numeric", length = n)
27	2x	f_x <- vector(mode = "numeric", length = n)
28	2x	for (i in seq_along(t_model)) {
29	6x	tt <- 0
30	6x	f <- 1
31
32	6x	iter <- 0
33	6x	while (abs(f) > tolerance) { # Newton loop
34	30x	slope <- (y[i] - y0) / (x[i] - x0)
35	30x	u <- exp(u235 * t0) - exp(u235 * tt)
36	30x	v <- exp(u238 * t0) - exp(u238 * tt)
37	30x	uprime <- -u235 * exp(u235 * tt)
38	30x	vprime <- -u238 * exp(u238 * tt)
39
40	30x	f <- u / v - slope * u238_235
41	30x	df <- (1 / v) * (uprime - u / v * vprime) # df/dtt
42	30x	dtt <- f / df # Newton-Raphson increment
43	30x	tt <- tt - dtt
44
45		# Loop counter
46	30x	iter <- iter + 1
47	30x	if (iter >= stop) {
48	!	msg <- tr_("Convergence not reached (possible infinite loop).")
49	!	warning(msg, call. = FALSE)
50	!	break
51		}
52		}
53
54	6x	t_model[i] <- tt
55	6x	f_x[i] <- f # Residual (should be zero)
56		}
57
58	2x	mu_i <- (x - x0) / (exp(u238 * t0) - exp(u238 * t_model))
59	2x	kappa_i <- (z - z0) / (exp(th232 * t0) - exp(th232 * t_model)) / mu_i
60
61	2x	data.frame(
62	2x	age = t_model * 1000, # Use Ma instead of Ga
63	2x	mu = mu_i,
64	2x	kappa = kappa_i,
65	2x	residual = f_x
66		)
67		}
68		)
69
70		#' @export
71		#' @rdname pb_age
72		#' @aliases pb_age,list,missing,missing-method
73		setMethod(
74		f = "pb_age",
75		signature = c(x = "list", y = "missing", z = "missing"),
76		definition = function(x, t0 = 3.8,
77		x_star = 18.75, y_star = 15.63, z_star = 38.86,
78		mu = 9.66, kappa = 3.90, th232 = 0.049475,
79		u238 = 0.155125, u235 = 0.98485, u238_235 = 137.79,
80		tolerance = sqrt(.Machine$double.eps), stop = 100) {
81		## Validation
82	3x	x <- grDevices::xyz.coords(x)
83	2x	methods::callGeneric(x = x$x, y = x$y, z = x$z, t0 = t0,
84	2x	x_star = x_star, y_star = y_star, z_star = z_star,
85	2x	mu = mu, kappa = kappa, th232 = th232,
86	2x	u238 = u238, u235 = u235, u238_235 = u238_235,
87	2x	tolerance = tolerance, stop = stop)
88		}
89		)

1		# INTERVAL ESTIMATION
2		#' @include AllGenerics.R
3		NULL
4
5		# HDR ==========================================================================
6		#' @export
7		#' @rdname interval_hdr
8		#' @aliases interval_hdr,CalibratedAges,missing-method
9		setMethod(
10		f = "interval_hdr",
11		signature = c(x = "CalibratedAges", y = "missing"),
12		definition = function(x, level = 0.954, ...) {
13		## Check
14	4x	c14_validate(x)
15
16	4x	hdr <- apply(
17	4x	X = x,
18	4x	MARGIN = 2,
19	4x	FUN = function(y, x, level, ...) {
20	!	if (all(is.na(y))) return(NULL)
21	5x	arkhe::interval_hdr(x, as.numeric(y), level, ...)
22		},
23	4x	x = aion::time(x, calendar = NULL),
24	4x	level = level,
25	4x	simplify = FALSE
26		)
27	4x	n <- vapply(X = hdr, FUN = nrow, FUN.VALUE = integer(1))
28	4x	lab <- rep(labels(x), n)
29
30	4x	hdr <- do.call(rbind, hdr)
31	4x	.CalibratedIntervals(
32	4x	.Id = lab,
33	4x	.Start = aion::as_fixed(hdr[, 1]),
34	4x	.End = aion::as_fixed(hdr[, 2]),
35	4x	p = unname(hdr[, 3])
36		)
37		}
38		)
39
40		# Credible =====================================================================
41		#' @export
42		#' @rdname interval_credible
43		#' @aliases interval_credible,CalibratedAges,missing-method
44		setMethod(
45		f = "interval_credible",
46		signature = c(x = "CalibratedAges"),
47		definition = function(x, level = 0.954, n = 100, ...) {
48		## Check
49	!	c14_validate(x)
50
51	!	spl <- c14_sample(x, n = n, calendar = NULL)
52	!	crd <- apply(
53	!	X = spl,
54	!	MARGIN = 2,
55	!	FUN = arkhe::interval_credible,
56	!	level = level,
57	!	simplify = FALSE
58		)
59
60	!	crd <- do.call(rbind, crd)
61	!	.CalibratedIntervals(
62	!	.Id = labels(x),
63	!	.Start = aion::as_fixed(crd[, 1]),
64	!	.End = aion::as_fixed(crd[, 2]),
65	!	p = unname(crd[, 3])
66		)
67		}
68		)

1		# 14C CALIBRATION CURVE
2		#' @include AllGenerics.R
3		NULL
4
5		# Get curve ====================================================================
6		#' @export
7		#' @rdname c14_curve
8		#' @aliases c14_curve,character-method
9		setMethod(
10		f = "c14_curve",
11		signature = c(name = "character"),
12		definition = function(name) {
13	13x	curve_ok <- c("intcal20", "intcal13", "intcal09", "intcal04",
14	13x	"marine20", "marine13", "marine09", "marine04",
15	13x	"shcal20", "shcal13", "shcal04")
16	13x	name <- match.arg(name, choices = curve_ok, several.ok = TRUE)
17
18	12x	curves <- lapply(X = name, FUN = read_curve)
19	12x	names(curves) <- name
20	12x	curves
21		}
22		)
23
24		#' @export
25		#' @rdname c14_curve
26		#' @aliases c14_curve,CalibratedAges-method
27		setMethod(
28		f = "c14_curve",
29		signature = c(name = "CalibratedAges"),
30		definition = function(name) {
31	!	methods::callGeneric(name@curves)
32		}
33		)
34
35		read_curve <- function(x) {
36		## Read data
37	12x	curve_dir <- system.file("extdata", package = "ananke")
38	12x	curve_path <- file.path(curve_dir, paste0(x, ".14c"))
39	12x	curve <- utils::read.table(curve_path, header = FALSE, sep = ",",
40	12x	dec = ".", strip.white = TRUE,
41	12x	comment.char = "#")
42
43	12x	curve <- curve[, c(1, 2, 3)]
44	12x	colnames(curve) <- c("CALBP", "AGE", "ERROR")
45	12x	curve
46		}
47
48		cite_curve <- function(x) {
49	1x	curve <- c(
50	1x	bomb04nh1 = "Hua and Barbetti 2004",
51	1x	bomb04nh2 = "Hua and Barbetti 2004",
52	1x	bomb04nh3 = "Hua and Barbetti 2004",
53	1x	bomb04sh = "Hua and Barbetti 2004",
54	1x	bomb13nh1 = "Hua, Berbetti and Rakowski 2013",
55	1x	bomb13nh2 = "Hua, Berbetti and Rakowski 2013",
56	1x	bomb13nh3 = "Hua, Berbetti and Rakowski 2013",
57	1x	bomb13sh12 = "Hua, Berbetti and Rakowski 2013",
58	1x	bomb13sh3 = "Hua, Berbetti and Rakowski 2013",
59	1x	bomb21nh1 = "Hua et al. 2022",
60	1x	bomb21nh2 = "Hua et al. 2022",
61	1x	bomb21nh3 = "Hua et al. 2022",
62	1x	bomb21sh12 = "Hua et al. 2022",
63	1x	bomb21sh3 = "Hua et al. 2022",
64	1x	cariaco04 = "Hughen et al. 2004",
65	1x	intcal98 = "Stuiver et al. 1998",
66	1x	intcal04 = "Reimer et al. 2004",
67	1x	intcal09 = "Reimer et al. 2009",
68	1x	intcal13 = "Reimer et al. 2013",
69	1x	intcal20 = "Reimer et al. 2020",
70	1x	kueppers04 = "Kueppers et al. 2004",
71	1x	marine98 = "Stuiver, Reimer and Braziunas 1998",
72	1x	marine04 = "Hughen et al. 2004",
73	1x	marine09 = "Reimer et al. 2009",
74	1x	marine13 = "Reimer et al. 2013",
75	1x	marine20 = "Heaton et al. 2020",
76	1x	shcal04 = "McCormac et al. 2004",
77	1x	shcal13 = "Hogg et al. 2013",
78	1x	shcal20 = "Hogg et al. 2020"
79		)
80
81	1x	curve[x]
82		}
83
84		# Approximate curve ============================================================
85		#' Interpolate 14C Calibration Curve
86		#'
87		#' @param name A [`character`] vector naming calibration curves.
88		#' @param out A [`numeric`] vector specifying where interpolation is to take
89		#' place.
90		#' @param F14C A [`logical`] scalar: should estimated F14C values be used
91		#' instead of radiocarbon ages?
92		#' @return
93		#' A `list` of `list` with the following elements:
94		#' \tabular{ll}{
95		#' `mu` \tab Interpolated values \cr
96		#' `tau` \tab Interpolated errors \cr
97		#' `max` \tab Maximum value of the calibration curve \cr
98		#' `min` \tab Minimum value of the calibration curve \cr
99		#' }
100		#' @keywords internal
101		#' @noRd
102		approx_curve <- function(name, out, F14C = FALSE) {
103		## Get data
104	12x	curve_data <- c14_curve(name)
105
106		## Interpolate
107	11x	lapply(
108	11x	X = curve_data,
109	11x	FUN = function(x, xout, F14C) {
110	11x	if (F14C) {
111	!	x_f14c <- c14_f14c(x[, 2], x[, 3])
112	!	x[, 2] <- x_f14c$value
113	!	x[, 3] <- x_f14c$error
114		}
115
116	11x	list(
117	11x	mu = stats::approx(x[, 1], x[, 2], xout = xout)$y,
118	11x	tau = stats::approx(x[, 1], x[, 3], xout = xout)$y,
119	11x	max = max(x[, 2]),
120	11x	min = min(x[, 2])
121		)
122		},
123	11x	xout = out,
124	11x	F14C = F14C
125		)
126		}

1		# SUBSET
2		#' @include AllGenerics.R
3		NULL
4
5		# Extract ======================================================================
6		## [ ---------------------------------------------------------------------------
7		#' @export
8		#' @rdname subset
9		#' @aliases [,CalibratedAges-method
10		setMethod(
11		f = "[",
12		signature = c(x = "CalibratedAges"),
13		function(x, i, j, k, drop = FALSE) {
14	41x	z <- x@.Data
15	41x	time <- x@.Time
16	41x	values <- x@values
17	41x	errors <- x@errors
18	41x	curves <- x@curves
19	41x	reservoir_offsets <- x@reservoir_offsets
20	41x	reservoir_errors <- x@reservoir_errors
21	41x	positions <- x@positions
22	41x	status <- x@status
23
24	41x	z <- z[i, j, k, drop = drop]
25	41x	if (!missing(i)) {
26	!	if (is.character(i)) i <- match(i, dimnames(x)[1L])
27	!	time <- time[i]
28		}
29	41x	if (!missing(j)) {
30	!	if (is.character(j)) j <- match(j, dimnames(x)[2L])
31	40x	values <- values[j]
32	40x	errors <- errors[j]
33	40x	curves <- curves[j]
34	40x	reservoir_offsets <- reservoir_offsets[j]
35	40x	reservoir_errors <- reservoir_errors[j]
36	40x	positions <- positions[j]
37	40x	status <- status[j]
38		}
39
40	19x	if (isTRUE(drop)) return(z)
41	22x	methods::initialize(x, z, .Time = time, values = values, errors = errors,
42	22x	curves = curves, reservoir_offsets = reservoir_offsets,
43	22x	reservoir_errors = reservoir_errors,
44	22x	positions = positions, status = status)
45		}
46		)

1		# SAMPLE 14C
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname c14_sample
7		#' @aliases c14_sample,CalibratedAges-method
8		setMethod(
9		f = "c14_sample",
10		signature = "CalibratedAges",
11		definition = function(object, n = 100, calendar = get_calendar()) {
12	1x	apply(
13	1x	X = object,
14	1x	MARGIN = 2,
15	1x	FUN = function(prob, size, x) {
16	2x	sample(x, size = size, replace = TRUE, prob = prob)
17		},
18	1x	x = aion::time(object, calendar = calendar),
19	1x	size = n
20		)
21		}
22		)

1		# MUTATORS
2		#' @include AllGenerics.R
3		NULL
4
5		# Get ==========================================================================
6		#' @export
7		#' @method labels CalibratedAges
8		labels.CalibratedAges <- function(object, ...) {
9	26x	colnames(object) %\|\|% paste0("X", NCOL(object))
10		}
11
12		#' @export
13		#' @rdname labels
14		#' @aliases labels,CalibratedAges-method
15		setMethod("labels", "CalibratedAges", labels.CalibratedAges)