kairos coverage - 73.16%

Files
Source

#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname seriate_rank
#' @aliases seriate_rank,data.frame-method
setMethod(
  f = "seriate_rank",
  signature = c(object = "data.frame"),
  definition = function(object, EPPM = FALSE, margin = c(1, 2), stop = 100) {
    object <- data.matrix(object)
    methods::callGeneric(object, EPPM = EPPM, margin = margin, stop = stop)
  }
)

#' @export
#' @rdname seriate_rank
#' @aliases seriate_rank,matrix-method
setMethod(
  f = "seriate_rank",
  signature = c(object = "matrix"),
  definition = function(object, EPPM = FALSE, margin = c(1, 2), stop = 100) {
    ## Validation
    margin <- as.integer(margin)
    stop <- as.integer(stop)

    data <- object
    if (EPPM) data <- tabula::eppm(object)

    ## Compute ranks
    ## margin = 1 : on rows
    ## margin = 2 : on columns
    reorder <- function(x, margin) {
      i <- seq_len(nrow(x))
      j <- seq_len(ncol(x))
      if (margin == 1) k <- colSums(t(x) * j) / rowSums(x)
      if (margin == 2) k <- colSums(x * i) / colSums(x)
      order(k)
    }

    start <- 0
    index <- list(rows = seq_len(nrow(data)), columns = seq_len(ncol(data)))
    convergence <- FALSE
    while (!convergence) {
      old_index <- index
      # Rearrange along margins
      for (k in margin) {
        index[[k]] <- index[[k]][reorder(data[index[[1]], index[[2]]], margin = k)]
      }
      # Loop counter
      convergence <- identical(index, old_index)
      start <- start + 1
      if (start >= stop) {
        msg <- tr_("Convergence not reached (possible infinite loop).")
        warning(msg, call. = FALSE)
        break
      }
    }

    # New PermutationOrder object
    .RankPermutationOrder(
      rows_order = as.integer(index[[1]]),
      columns_order = as.integer(index[[2]])
    )
  }
)

# SUBSET
#' @include AllGenerics.R
NULL

# Extract ======================================================================
## [ ---------------------------------------------------------------------------
#' @export
#' @rdname subset
#' @aliases [,MeanDate-method
setMethod(
  f = "[",
  signature = c(x = "MeanDate"),
  function(x, i, j, k, drop = FALSE) {
    z <- x@.Data
    time <- x@.Time
    dates <- x@dates

    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])
      dates <- dates[j]
    }

    if (isTRUE(drop)) return(z)
    methods::initialize(x, z, .Time = time, dates = dates)
  }
)

#' @export
#' @rdname subset
#' @aliases [,IncrementTest-method
setMethod(
  f = "[",
  signature = c(x = "IncrementTest"),
  function(x, i, j, k, drop = FALSE) {
    z <- x@.Data
    time <- x@.Time
    statistic <- x@statistic
    p_value <- x@p_value

    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])
      statistic <- statistic[j]
      p_value <- p_value[j]
    }

    if (isTRUE(drop)) return(z)
    methods::initialize(x, z, .Time = time, statistic = statistic, p_value = p_value)
  }
)

## [[ --------------------------------------------------------------------------
#' Extract Parts of an Object
#'
#' @inheritParams [[
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
extract_slot <- function(x, i) {
  class_name <- class(x)
  i <- match.arg(i, choices = methods::slotNames(class_name),
                 several.ok = FALSE)
  data <- methods::slot(x, i)
  data
}

#' @export
#' @rdname subset
#' @aliases [[,PermutationOrder,ANY,missing-method
setMethod(
  f = "[[",
  signature = c(x = "PermutationOrder", i = "ANY", j = "missing"),
  definition = extract_slot
)

# MEAN CERAMIC DATE
#' @include AllClasses.R AllGenerics.R
NULL

# MCD ==========================================================================
#' @export
#' @rdname mcd
#' @aliases mcd,numeric,numeric-method
setMethod(
  f = "mcd",
  signature = c(object = "numeric", dates = "numeric"),
  definition = function(object, dates, calendar = CE()) {
    object <- matrix(object, nrow = 1)
    methods::callGeneric(object, dates = dates, calendar = calendar)
  }
)

#' @export
#' @rdname mcd
#' @aliases mcd,data.frame,numeric-method
setMethod(
  f = "mcd",
  signature = c(object = "data.frame", dates = "numeric"),
  definition = function(object, dates, calendar = CE()) {
    object <- data.matrix(object)
    methods::callGeneric(object, dates = dates, calendar = calendar)
  }
)

#' @export
#' @rdname mcd
#' @aliases mcd,matrix,numeric-method
setMethod(
  f = "mcd",
  signature = c(object = "matrix", dates = "numeric"),
  definition = function(object, dates, calendar = CE()) {
    ## Validation
    arkhe::assert_length(dates, NCOL(object))

    if (is.null(calendar)) {
      dates <- aion::as_fixed(dates) # Assume rata die
    } else {
      dates <- aion::fixed(dates, calendar = calendar)
    }

    ## Calculate MCD
    mcd_dates <- apply(X = object, MARGIN = 1, FUN = .mcd, dates = dates)
    mcd_dates <- aion::as_fixed(mcd_dates)

    ts <- aion::series(object, time = mcd_dates)
    rownames(ts) <- rownames(object)
    .MeanDate(ts, dates = dates)
  }
)

.mcd <- function(x, dates) {
  stats::weighted.mean(x = dates, w = x)
}

# Resample =====================================================================
#' @export
#' @rdname resample_mcd
#' @aliases bootstrap,MeanDate-method
setMethod(
  f = "bootstrap",
  signature = c(object = "MeanDate"),
  definition = function(object, n = 1000, f = NULL,
                        calendar = get_calendar()) {

    w <- object
    m <- nrow(w)
    p <- ncol(w)
    seq_col <- seq_len(p)

    dates <- aion::as_year(object@dates, calendar = calendar)
    theta <- function(x, counts, dates) {
      .mcd(counts[x], dates[x])
    }

    results <- vector(mode = "list", length = m)
    for (i in seq_len(m)) {
      results[[i]] <- arkhe::bootstrap(
        object = seq_col,
        do = theta,
        n = n,
        counts = w[i, , 1, drop = TRUE],
        dates = dates,
        f = f
      )
    }

    results <- do.call(rbind, results)
    rownames(results) <- rownames(w)
    as.data.frame(results)
  }
)

#' @export
#' @rdname resample_mcd
#' @aliases jackknife,MeanDate-method
setMethod(
  f = "jackknife",
  signature = c(object = "MeanDate"),
  definition = function(object, f = NULL,
                        calendar = get_calendar()) {

    w <- object
    m <- nrow(w)
    p <- ncol(w)

    dates <- aion::as_year(object@dates, calendar = calendar)
    theta <- function(x, counts, dates) {
      .mcd(counts[x], dates[x])
    }

    results <- vector(mode = "list", length = m)
    for (i in seq_len(m)) {
      results[[i]] <- arkhe::jackknife(
        object = seq_len(p),
        do = theta,
        counts = w[i, , 1, drop = TRUE],
        dates = dates,
        f = f
      )
    }
    results <- do.call(rbind, results)
    rownames(results) <- rownames(w)
    results
  }
)

#' @export
#' @rdname resample_mcd
#' @aliases simulate,MeanDate-method
setMethod(
  f = "simulate",
  signature = c(object = "MeanDate"),
  definition = function(object, nsim = 1000) {

    results <- apply(
      X = object[, , 1, drop = TRUE],
      MARGIN = 1,
      FUN = resample,
      do = .mcd,
      n = nsim,
      dates = object@dates,
      f = function(x) x
    )

    .SimulationMeanDate(object, replications = t(results))
  }
)

# Plot =========================================================================
#' @export
#' @method plot MeanDate
plot.MeanDate <- function(x, calendar = get_calendar(),
                          decreasing = TRUE,
                          main = NULL, sub = NULL,
                          ann = graphics::par("ann"), axes = TRUE,
                          frame.plot = axes,
                          panel.first = NULL, panel.last = NULL, ...) {
  ## Get data
  n <- NROW(x)
  n_seq <- seq_len(n)
  sites <- rownames(x) %||% paste0("S", n_seq)

  ## Graphical parameters
  col <- list(...)$col %||% graphics::par("col")
  bg <- list(...)$bg %||% graphics::par("bg")
  pch <- list(...)$pch %||% c(16)
  cex <- list(...)$cex %||% graphics::par("cex")
  lwd <- list(...)$lwd %||% graphics::par("lwd")
  cex.axis <- list(...)$cex.axis %||% graphics::par("cex.axis")
  col.axis <- list(...)$col.axis %||% graphics::par("col.axis")
  font.axis <- list(...)$font.axis %||% graphics::par("font.axis")
  if (length(col) != n) col <- rep(col, length.out = n)
  if (length(bg) != n) bg <- rep(bg, length.out = n)
  if (length(pch) != n) pch <- rep(pch, length.out = n)
  if (length(cex) != n) cex <- rep(cex, length.out = n)
  if (length(lwd) != n) lwd <- rep(lwd, length.out = n)

  ## Save and restore
  mar <- graphics::par("mar")
  mar[2] <- inch2line(sites, cex = cex.axis) + 0.5
  old_par <- graphics::par(mar = mar)
  on.exit(graphics::par(old_par))

  ## Open new window
  grDevices::dev.hold()
  on.exit(grDevices::dev.flush(), add = TRUE)
  graphics::plot.new()

  ## Set plotting coordinates
  years <- time(x, calendar = calendar)
  xlim <- range(years)
  ylim <- c(1, n)
  graphics::plot.window(xlim = xlim, ylim = ylim)

  ## Evaluate pre-plot expressions
  panel.first

  ## Order data
  k <- order(years, decreasing = decreasing)

  ## Plot
  graphics::points(x = years[k], y = n_seq, col = col[k],
                   bg = bg[k], pch = pch[k], cex = cex[k], lwd = lwd[k])

  ## Evaluate post-plot and pre-axis expressions
  panel.last

  ## Construct Axis
  if (axes) {
    aion::year_axis(side = 1, format = TRUE, calendar = calendar,
                    current_calendar = calendar, cex.axis = cex.axis,
                    col.axis = col.axis, font.axis = font.axis)
    graphics::axis(side = 2, at = n_seq, labels = sites[k],
                   xpd = NA, cex.axis = cex.axis, las = 1,
                   col.axis = col.axis, font.axis = font.axis)
  }

  ## Plot frame
  if (frame.plot) {
    graphics::box()
  }

  ## Add annotation
  if (ann) {
    xlab <- format(calendar)
    ylab <- NULL
    graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
  }

  invisible(x)
}

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

#' @export
#' @method plot SimulationMeanDate
plot.SimulationMeanDate <- function(x, calendar = get_calendar(),
                                    interval = "student", level = 0.80,
                                    decreasing = TRUE,
                                    main = NULL, sub = NULL,
                                    ann = graphics::par("ann"), axes = TRUE,
                                    frame.plot = axes,
                                    panel.first = NULL, panel.last = NULL, ...) {
  ## Get data
  n <- NROW(x)
  n_seq <- seq_len(n)
  sites <- rownames(x) %||% paste0("S1", n_seq)

  ## Compute confidence interval
  fun <- function(x) conf(x, level = level, type = interval)
  inter <- apply(X = x@replications, MARGIN = 1, FUN = fun)
  inter <- apply(X = inter, MARGIN = 1, FUN = as_year, calendar = calendar)

  ## Graphical parameters
  col <- list(...)$col %||% graphics::par("col")
  bg <- list(...)$bg %||% graphics::par("bg")
  pch <- list(...)$pch %||% c(16)
  cex <- list(...)$cex %||% graphics::par("cex")
  lty <- list(...)$lty %||% graphics::par("lty")
  lwd <- list(...)$lwd %||% graphics::par("lwd")
  if (length(col) != n) col <- rep(col, length.out = n)
  if (length(bg) != n) bg <- rep(bg, length.out = n)
  if (length(pch) != n) pch <- rep(pch, length.out = n)
  if (length(cex) != n) cex <- rep(cex, length.out = n)
  if (length(lty) != n) lty <- rep(lty, length.out = n)
  if (length(lwd) != n) lwd <- rep(lwd, length.out = n)
  cex.axis <- list(...)$cex.axis %||% graphics::par("cex.axis")
  col.axis <- list(...)$col.axis %||% graphics::par("col.axis")
  font.axis <- list(...)$font.axis %||% graphics::par("font.axis")

  ## Save and restore
  mar <- graphics::par("mar")
  mar[2] <- inch2line(sites, cex = cex.axis) + 0.5
  old_par <- graphics::par(mar = mar)
  on.exit(graphics::par(old_par))

  ## Open new window
  grDevices::dev.hold()
  on.exit(grDevices::dev.flush(), add = TRUE)
  graphics::plot.new()

  ## Set plotting coordinates
  years <- time(x, calendar = calendar)
  xlim <- range(inter)
  ylim <- c(1, n)
  graphics::plot.window(xlim = xlim, ylim = ylim)

  ## Evaluate pre-plot expressions
  panel.first

  ## Order data
  k <- order(years, decreasing = decreasing)

  ## Plot
  graphics::segments(x0 = inter[k, 2], y0 = n_seq,
                     x1 = inter[k, 3], y1 = n_seq,
                     col = col[k], lty = lty[k], lwd = lwd[k])
  graphics::points(x = years[k], y = n_seq, col = col[k],
                   bg = bg[k], pch = pch[k], cex = cex[k], lwd = lwd[k])

  ## Evaluate post-plot and pre-axis expressions
  panel.last

  ## Construct Axis
  if (axes) {
    aion::year_axis(side = 1, format = TRUE, calendar = calendar,
                    current_calendar = calendar, cex.axis = cex.axis,
                    col.axis = col.axis, font.axis = font.axis)
    graphics::axis(side = 2, at = n_seq, labels = sites[k],
                   xpd = NA, cex.axis = cex.axis, las = 1,
                   col.axis = col.axis, font.axis = font.axis)
  }

  ## Plot frame
  if (frame.plot) {
    graphics::box()
  }

  ## Add annotation
  if (ann) {
    ## Caption
    cap <- sprintf(tr_("Simulated assemblages, %d replications."), NCOL(x@replications))
    xlab <- format(calendar)
    ylab <- NULL
    graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
  }

  invisible(x)
}

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

conf <- function(x, type = c("percentiles", "student", "normal", "range"),
                 level = 0.80) {
  type <- match.arg(type, several.ok = FALSE)
  if (type == "range") {
    conf <- range(x, na.rm = FALSE)
  } else if (type == "percentiles") {
    ## Confidence interval as described in Kintigh 1989
    k <- (1 - level) / 2
    conf <- stats::quantile(x, probs = c(k, 1 - k), names = FALSE)
  } else {
    ## Confidence interval
    conf <- arkhe::confidence_mean(x, level = level, type = type)
  }

  result <- c(mean(x), conf)
  names(result) <- c("mean", "lower", "upper")
  result
}

#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname assess
#' @aliases assess,AveragePermutationOrder-method
setMethod(
  f = "assess",
  signature = c(object = "AveragePermutationOrder"),
  definition = function(object, axes = 1, n = 1000,
                        progress = getOption("kairos.progress")) {
    ## Validation
    arkhe::assert_length(axes, 1)
    arkhe::assert_length(n, 1)

    ## Reorder along 'axes'
    data <- dimensio::get_data(object)
    perm <- permute(data, object)

    ## Number of local maxima
    freq <- data / rowSums(data)
    a <- sum(apply(X = freq, MARGIN = 2, FUN = local_maxima))

    b <- numeric(n)
    if (n > 0) {
      progress_bar <- interactive() && isTRUE(progress)
      if (progress_bar) pbar <- utils::txtProgressBar(max = n, style = 3)

      for (i in seq_len(n)) {
        ## Randomize original data
        i_obj <- apply(X = data, MARGIN = 2, FUN = sample, replace = FALSE)

        ## Number of local maxima
        i_ca <- seriate_average(i_obj, margin = c(1, 2), axes = axes)
        i_perm <- permute(i_obj, i_ca)
        i_freq <- i_perm / rowSums(i_perm)

        b[[i]] <- sum(apply(X = i_freq, MARGIN = 2, FUN = local_maxima))

        if (progress_bar) utils::setTxtProgressBar(pbar, i)
      }

      if (progress_bar) close(pbar)
    }

    ## Seriation coefficient
    E <- ncol(data)
    M <- nrow(data) * ifelse(E%%2 > 0, E + 1, E) / 2
    S <- (M - a) / (M - E)

    list(
      random = b,
      observed = a,
      expected = E,
      maximum = M,
      coef = S
    )
  }
)

local_maxima <- function(x) {
  n <- length(x)
  left <- c(0, x[-n])
  right <- c(x[-1L], 0)
  sum(x > left & x > right)
}

# HELPERS

## https://michaelchirico.github.io/potools/articles/developers.html
tr_ <- function(...) {
  enc2utf8(gettext(paste0(...), domain = "R-kairos"))
}

make_par <- function(params, x, n = 0) {
  p <- params[[x]] %||% graphics::par()[[x]]
  if (n > 0) p <- rep(p, length.out = n)
  p
}

#' Plotting Dimensions of Character Strings
#'
#' Convert string length in inch to number of (margin) lines.
#' @param x A [`character`] vector of string whose length is to be calculated.
#' @param ... Further parameter to be passed to [graphics::strwidth()]`, such as
#'  `cex`.
#' @return
#'  A [`numeric`] vector (maximum string width in units of margin lines).
#' @note For internal use only.
#' @family graphic tools
#' @keywords internal
#' @noRd
inch2line <- function(x, ...) {
  (max(graphics::strwidth(x, units = "inch", ...)) /
     graphics::par("cin")[2] + graphics::par("mgp")[2]) * graphics::par("cex")
}

#' Compute x Limits
#'
#' Computes x limits for a time series according to a given calendar.
#' This ensures that the x axis is always in chronological order.
#' @param x A [`TimeSeries-class`] object.
#' @param calendar A [`TimeScale-class`] object.
#' @return A length-two [`numeric`] vector.
#' @keywords internal
#' @noRd
xlim <- function(x, calendar) {
  x <- range(time(x, calendar = NULL))
  if (is.null(calendar)) return(x)
  as_year(x, calendar = calendar)
}

#' Indices of a rolling window
#'
#' @param x An object.
#' @param window An [`integer`] scalar giving the window size.
#' @return A [`list`] with the following components:
#'  \describe{
#'   \item{i}{An [`integer`] vector of indices.}
#'   \item{w}{An [`integer`] vector of indices giving the indice of the window
#'   mid-point.}
#'  }
#' @keywords internal
#' @noRd
roll <- function(x, window = 3) {
  ## Validation
  arkhe::assert_odd(window)

  n <- if (is.matrix(x) || is.data.frame(x)) nrow(x) else length(x)
  i <- seq_len(n) # Indices of the rows

  ## Matrix of rolling-window indices of length w
  w <- stats::embed(i, window)[, window:1]
  inds <- as.vector(t(w)) # Flatten indices

  ## Window mid-point
  m <- w[, ceiling(window / 2)]

  list(i = inds, w = rep(m, each = window))
}

#' Resample
#'
#' Simulates observations from a multinomial distribution.
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param do A [`function`] that takes `x` as an argument
#'  and returns a single numeric value.
#' @param n A non-negative [`integer`] specifying the number of bootstrap
#'  replications.
#' @param size A non-negative [`integer`] specifying the sample size.
#' @param f A [`function`] that takes a single numeric vector (the result of
#'  `do`) as argument.
#' @param ... Extra arguments passed to `do`.
#' @return
#'  If `f` is `NULL`, `resample()` returns the `n` values of `do`. Else,
#'  returns the result of `f` applied to the `n` values of `do`.
#' @seealso [stats::rmultinom()]
#' @keywords internal
#' @noRd
resample <- function(x, do, n, size = sum(x), ..., f = NULL) {
  ## Validation
  arkhe::assert_count(x)

  prob <- x / sum(x)
  replicates <- stats::rmultinom(n, size = size, prob = prob)
  values <- apply(X = replicates, MARGIN = 2, FUN = do, ...)
  if (is.function(f)) values <- f(values)
  values
}

#' 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
}

# CHRONOLOGICAL APPORTIONING
#' @include AllGenerics.R AllClasses.R
NULL

#' @export
#' @rdname apportion
#' @aliases apportion,data.frame-method
setMethod(
  f = "apportion",
  signature = c(object = "data.frame"),
  definition = function(object, s0, s1, t0, t1, from = min(s0), to = max(s1),
                        step = 25, method = c("uniform", "truncated"), z = 2,
                        progress = getOption("kairos.progress")) {
    object <- data.matrix(object)
    methods::callGeneric(object, s0, s1, t0, t1, from = from, to = to,
                         step = step, method = method, z = z,
                         progress = progress)
  }
)

#' @export
#' @rdname apportion
#' @aliases apportion,matrix-method
setMethod(
  f = "apportion",
  signature = c(object = "matrix"),
  definition = function(object, s0, s1, t0, t1, from = min(s0), to = max(s1),
                        step = 25, method = c("uniform", "truncated"), z = 2,
                        progress = getOption("kairos.progress")) {
    ## Validation
    method <- match.arg(method, several.ok = FALSE)
    arkhe::assert_lower(s0, s1, strict = FALSE)
    arkhe::assert_lower(t0, t1, strict = FALSE)
    arkhe::assert_lower(from, to, strict = TRUE)

    ## Get data
    n_site <- nrow(object)
    n_type <- ncol(object)
    span <- to - from

    ## Number of periods (rounded toward the smallest integer)
    n_periode <- ceiling(span / step)
    t <- seq_len(n_periode)
    periode <- paste(from + (t - 1) * step, from + t * step, sep = "_")

    ## Empty array to store apportioning probabilities
    p <- array(data = NA_real_, dim = c(nrow(object), ncol(object), n_periode))
    dimnames(p) <- c(dimnames(object), list(periode))
    a <- p

    ## Round the site dates to be multiples of the step size
    ## (rounded toward the smallest/largest multiple)
    s0[s0 %% step != 0] <- floor(s0[s0 %% step != 0] / step) * step
    s1[s1 %% step != 0] <- ceiling(s1[s1 %% step != 0] / step) * step

    ## Type midpoints and life spans
    m <- (t1 + t0) / 2
    g <- t1 - t0

    ## Distribution function
    fun <- switch (
      method,
      uniform = dist_uniform,
      truncated = dist_truncated
    )

    ## Apportion
    k_site <- seq_len(n_site)
    k_type <- seq_len(n_type)

    progress_bar <- interactive() && isTRUE(progress)
    if (progress_bar) pbar <- utils::txtProgressBar(max = n_site, style = 3)

    for (i in k_site) {
      for (j in k_type) {
        ## If the type lies outside the known site occupation or study interval
        if (t1[j] <= s0[i] | t0[j] >= s1[i]) next # Do not apportion

        ## Earliest and latest overlaps between the ware and site
        vj0 <- max(s0[i], t0[j])
        vj1 <- min(s1[i], t1[j])

        ## Get the apportioning probability for any period
        qjt0 <- vapply(X = s0[i] + (t - 1) * step, FUN = max,
                       FUN.VALUE = numeric(1), t0[j])
        qjt1 <- vapply(X = s0[i] + t * step, FUN = min,
                       FUN.VALUE = numeric(1), t1[j])

        p[i, j, ] <- fun(vj0, vj1, qjt0, qjt1, g[j], m[j], z)
      }

      if (progress_bar) utils::setTxtProgressBar(pbar, i)
    }

    if (progress_bar) close(pbar)

    ## Remove negative values
    p[p < 0] <- 0

    ## Apportion
    a[] <- apply(X = p, MARGIN = 3, FUN = function(x, counts) x * counts,
                 counts = object)

    .CountApportion(
      a,
      p = p,
      method = method,
      from = from,
      to = to,
      step = step
    )
  }
)

dist_uniform <- function(v_j0, v_j1, q_jt0, q_jt1, ...) {
  (q_jt1 - q_jt0) / (v_j1 - v_j0)
}
dist_truncated <- function(v_j0, v_j1, q_jt0, q_jt1, g, m, z) {
  z_j0 <- (v_j0 - m) / (g / (2 * z))
  z_j1 <- (v_j1 - m) / (g / (2 * z))

  z_jt0 <- (q_jt0 - m) / (g / (2 * z))
  z_jt1 <- (q_jt1 - m) / (g / (2 * z))

  (phi(z_jt1, z) - phi(z_jt0, z)) / (phi(z_j1, z) - phi(z_j0, z))
}
phi <- function(x, z) {
  extraDistr::ptnorm(x, mean = 0, sd = 1, a = -z, b = z)
}

# AORISTIC ANALYSIS
#' @include AllClasses.R AllGenerics.R
NULL

# Aoristic sum =================================================================
#' @export
#' @rdname aoristic
#' @aliases aoristic,numeric,numeric-method
setMethod(
  f = "aoristic",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, step = 1, start = min(x), end = max(y),
                        calendar = CE(), weight = TRUE, groups = NULL) {
    ## Validation
    n <- length(x)
    arkhe::assert_length(y, n)
    if (is.null(groups)) groups <- rep(NA, n)
    groups <- as.character(groups)
    arkhe::assert_length(groups, n)

    ## Grouping
    lvl <- unique(groups)
    grp <- factor(groups, levels = lvl, exclude = NULL)
    spl <- split(seq_along(x), f = grp)
    m <- length(spl)

    ## Convert to rata die
    ## NB: argument are not evaluated until they are used.
    ## ('start' and 'end' must go first)
    breaks <- seq(from = start, to = end, by = step)
    breaks <- aion::fixed(breaks, calendar = calendar)
    x <- aion::fixed(x, calendar = calendar)
    y <- aion::fixed(y, calendar = calendar)

    ## Set up breaks
    block_start <- utils::head(breaks, -1)
    block_end <- utils::tail(breaks, -1)
    block_mid <- (block_start + block_end) / 2

    ## Number of time blocks
    n_blocks <- length(breaks) - 1
    ao_probs <- array(data = 0, dim = c(n, n_blocks, m))

    ## Aoristic probability
    span <- abs(y - x)
    for (k in seq_len(m)) {
      s <- spl[[k]]
      for (j in seq_len(n_blocks)) {
        a <- block_start[[j]]
        b <- block_end[[j]]

        max_start <- pmax(x[s], a)
        min_end <- pmin(y[s], b)
        overlap <- (min_end - max_start) / span[s]
        overlap[overlap <= 0] <- 0
        if (!weight) overlap[overlap > 0] <- 1
        ao_probs[s, j, k] <- overlap
      }
    }

    ## Aoristic sum
    ao_sum <- apply(X = ao_probs, MARGIN = c(2, 3), FUN = sum)

    if (!anyNA(groups)) {
      dimnames(ao_probs)[[3]] <- colnames(ao_sum) <- lvl
    }

    ts <- aion::series(
      ao_sum,
      time = block_mid
    )
    .AoristicSum(
      ts,
      breaks = breaks,
      span = span,
      groups = groups,
      p = ao_probs
    )
  }
)

#' @export
#' @rdname aoristic
#' @aliases aoristic,ANY,missing-method
setMethod(
  f = "aoristic",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, step = 1, start = NULL, end = NULL,
                        calendar = CE(), weight = TRUE, groups = NULL) {
    xy <- grDevices::xy.coords(x)
    g <- groups
    if (is.list(x) && !is.null(g) && length(g) == 1) {
      g <- x[[g]]
      if (is.null(g)) {
        txt <- tr_("%s is a list, but does not have component %s.")
        msg <- sprintf(txt, sQuote("x"), sQuote(groups))
        stop(msg, call. = FALSE)
      }
    }

    ## Start/end
    if (is.null(start)) start <- min(xy$x, na.rm = TRUE)
    if (is.null(end)) end <- max(xy$y, na.rm = TRUE)

    methods::callGeneric(x = xy$x, y = xy$y, step = step, start = start,
                         end = end, calendar = calendar, weight = weight,
                         groups = g)
  }
)

# Rate of change ===============================================================
#' @export
#' @rdname roc
#' @aliases roc,AoristicSum-method
setMethod(
  f = "roc",
  signature = c(object = "AoristicSum"),
  definition = function(object, n = 100) {
    ## Get time span of the blocks
    mid <- aion::time(object, calendar = NULL)
    step <- mean(diff(mid))
    end <- utils::tail(mid, -1) + step / 2

    ## Get aoristic weights
    w <- object@p

    n <- as.integer(n)
    m <- dim(w)[[1]]
    p <- dim(w)[[2]]
    q <- dim(w)[[3]]

    ## Monte Carlo simulation
    roc <- array(data = 0, dim = c(p-1, n, q))
    for (j in seq_len(q)) {
      sim <- array(data = 0, dim = c(m, p, n))
      for (i in seq_len(m)) {
        probs <- w[i, , j]
        if (sum(probs) == 0) next
        counts <- replicate(n, {
          tmp <- numeric(p)
          k <- sample(seq_len(p), size = 1, prob = probs)
          tmp[k] <- 1
          tmp
        })
        sim[i, , ] <- counts
      }
      total <- apply(X = sim, MARGIN = c(2, 3), FUN = sum)
      roc[, , j] <- apply(X = total, MARGIN = 2, FUN = diff, lag = 1)
    }

    dimnames(roc)[[3]] <- dimnames(w)[[3]]

    ts <- aion::series(
      object = roc / step,
      time = end
    )
    .RateOfChange(
      ts,
      replicates = n,
      groups = unique(object@groups)
    )
  }
)

# Plot =========================================================================
## AoristicSum -----------------------------------------------------------------
#' @export
#' @method plot AoristicSum
plot.AoristicSum <- function(x, calendar = get_calendar(),
                             type = c("bar"),
                             flip = FALSE, ncol = NULL,
                             main = NULL, sub = NULL,
                             ann = graphics::par("ann"), axes = TRUE,
                             frame.plot = axes,
                             panel.first = NULL, panel.last = NULL, ...) {
  ## Validation
  type <- match.arg(type, several.ok = FALSE)

  ## Plot
  panel_bar <- function(x, y, ...) {
    step <- mean(diff(x))
    start <- x - step / 2
    end <- x + step / 2
    graphics::rect(xleft = start, xright = end, ybottom = 0, ytop = y, ...)
  }
  panel_density <- function(x, y, ...) {
    step <- mean(diff(x))
    s <- sample(x = x, size = 10^5, replace = TRUE, prob = y)
    d <- stats::density(x = s, bw = step)
    x <- d$x
    y <- d$y
    graphics::polygon(x = c(x, rev(x)), y = c(y, numeric(length(y))),
                      border = NA, ...)
    graphics::lines(x = x, y = y)
  }

  panel_fun <- switch (
    type,
    bar = panel_bar,
    density = panel_density
  )

  ## Method for TimeSeries
  methods::callNextMethod(x, facet = "multiple", calendar = calendar,
                          panel = panel_fun, flip = flip, ncol = ncol,
                          xlab = NULL, ylab = NULL,
                          main = main, sub = sub,
                          ann = ann, axes = axes,
                          frame.plot = frame.plot,
                          panel.first = panel.first,
                          panel.last = panel.last, ...)


  invisible(x)
}

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

#' @export
#' @method image AoristicSum
image.AoristicSum <- function(x, calendar = get_calendar(), ...) {
  methods::callNextMethod()
  invisible(x)
}

#' @export
#' @rdname plot_aoristic
#' @aliases image,AoristicSum-method
setMethod("image", c(x = "AoristicSum"), image.AoristicSum)

## RateOfChange ----------------------------------------------------------------
#' @export
#' @method plot RateOfChange
plot.RateOfChange <- function(x, calendar = get_calendar(),
                              level = 0.95, flip = FALSE, ncol = NULL,
                              main = NULL, sub = NULL,
                              ann = graphics::par("ann"), axes = TRUE,
                              frame.plot = axes,
                              panel.first = NULL, panel.last = NULL, ...) {
  ## Validation
  n <- dim(x)[[3L]]
  n_seq <- seq_len(n)
  if (is.null(ncol)) ncol <- if (n > 4) 2 else 1
  n_row <- ceiling(n / ncol)
  ylabs <- dimnames(x)[[3L]] %||% c(tr_("Rate of Change"))

  ## Compute mean and confidence interval
  ci <- apply(
    X = x,
    MARGIN = c(1, 3),
    FUN = function(x, level) {
      ci <- arkhe::confidence_mean(x, level = level)
      c(mean = mean(x), ci)
    },
    level = level
  )

  ## Graphical parameters
  ## Save and restore
  old_par <- graphics::par(
    mar = c(0, 5.1, 0, if (flip) 5.1 else 2.1),
    oma = c(6, 0, 5, 0),
    mfcol = c(n_row, ncol)
  )
  on.exit(graphics::par(old_par))

  cex.lab <- list(...)$cex.lab %||% graphics::par("cex.lab")
  col.lab <- list(...)$col.lab %||% graphics::par("col.lab")
  font.lab <- list(...)$font.lab %||% graphics::par("font.lab")
  cex.axis <- list(...)$cex.axis %||% graphics::par("cex.axis")
  col.axis <- list(...)$col.axis %||% graphics::par("col.axis")
  font.axis <- list(...)$font.axis %||% graphics::par("font.axis")
  cex.main <- list(...)$cex.main %||% graphics::par("cex.main")
  font.main <- list(...)$font.main %||% graphics::par("font.main")
  col.main <- list(...)$col.main %||% graphics::par("col.main")
  col <- list(...)$col %||% c("grey")
  if (length(col) != n) col <- rep(col, n)

  years <- aion::time(x, calendar = calendar)
  xlim <- range(years)
  for (i in n_seq) {
    xi <- ci[, , i, drop = TRUE]

    ## Open new window
    grDevices::dev.hold()
    on.exit(grDevices::dev.flush(), add = TRUE)
    graphics::plot.new()

    ## Set plotting coordinates
    ylim <- range(xi, na.rm = TRUE)
    graphics::plot.window(xlim = xlim, ylim = ylim)

    ## Evaluate pre-plot expressions
    panel.first

    ## Plot
    graphics::abline(h = 0, lty = "dotted", lwd = 1)
    graphics::polygon(x = c(years, rev(years)),
                      y = c(xi[2, ], rev(xi[3, ])),
                      col = col[i], border = NA)
    graphics::lines(x = years, y = xi[1, ])
    graphics::points(x = years, y = xi[1, ], pch = 16)

    ## Evaluate post-plot and pre-axis expressions
    panel.last

    ## Construct Axis
    do_x <- i %% n_row == 0 || i == n
    y_side <- if (i %% 2 || !flip) 2 else 4
    if (axes) {
      if (do_x) {
        year_axis(side = 1, format = TRUE, calendar = calendar,
                  current_calendar = calendar, cex.axis = cex.axis,
                  col.axis = col.axis, font.axis = font.axis)
      }
      graphics::axis(side = y_side, xpd = NA, cex.axis = cex.axis,
                     col.axis = col.axis, font.axis = font.axis, las = 1)
    }

    ## Plot frame
    if (frame.plot) {
      graphics::box()
    }

    ## Add annotation
    if (ann) {
      if (do_x) {
        xlab <- format(calendar)
        graphics::mtext(xlab, side = 1, line = 3, cex = cex.lab, col = col.lab,
                        font = font.lab)
      }
      graphics::mtext(ylabs[[i]], side = y_side, line = 4, cex = cex.lab,
                      col = col.lab, font = font.lab)
    }
  }

  ## Add annotation
  if (ann) {
    graphics::par(mfcol = c(1, 1))
    graphics::mtext(main, side = 3, line = 3, cex = cex.main, font = font.main,
                    col = col.main)
  }

  invisible(x)
}

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

# FREQUENCY INCREMENT TEST
#' @include AllGenerics.R AllClasses.R
NULL

# FIT ==========================================================================
#' @export
#' @rdname fit
#' @aliases fit,data.frame,numeric-method
setMethod(
  f = "fit",
  signature = c(object = "data.frame", dates = "numeric"),
  definition = function(object, dates, calendar = CE(),
                        level = 0.95, roll = FALSE, window = 3) {
    object <- data.matrix(object)
    methods::callGeneric(object, dates = dates, calendar = calendar,
                         level = level, roll = roll, window = window)
  }
)

#' @export
#' @rdname fit
#' @aliases fit,matrix,numeric-method
setMethod(
  f = "fit",
  signature = c(object = "matrix", dates = "numeric"),
  definition = function(object, dates, calendar = CE(),
                        level = 0.95, roll = FALSE, window = 3) {
    ## Validation
    alpha <- 1 - level
    half_window <- (window - 1) / 2
    arkhe::assert_length(dates, nrow(object))

    ## Convert to rata die
    if (is.null(calendar)) {
      dates <- aion::as_fixed(dates)
    } else {
      dates <- aion::fixed(dates, calendar = calendar)
    }

    ## Compute test
    results <- test_FIT(object, dates, roll = FALSE)
    results <- results[[1L]]

    ## Check results
    failed <- is.na(results$p.value)
    if (any(failed)) {
      txt <- ngettext(sum(failed), "%d element was skipped: %s.",
                      "%d elements were skipped: %s.")
      msg <- sprintf(txt, sum(failed),
                     paste0(rownames(results)[failed], collapse = ", "))
      warning(msg, call. = FALSE)
    }

    ## Rolling window
    fit_roll <- matrix(FALSE, nrow = nrow(object), ncol = ncol(object))
    if (roll) {
      results_roll <- test_FIT(object, dates, roll = roll, window = window)
      results_roll <- lapply(
        X = results_roll,
        FUN = function(x, alpha) { x$p.value <= alpha },
        alpha = alpha
      )
      fit_row <- seq(from = half_window + 1, to = nrow(fit_roll) - half_window)
      fit_roll[fit_row, ] <- do.call(rbind, results_roll)

      fit_roll <- apply(
        X = fit_roll,
        MARGIN = 2,
        FUN = function(x, half_window) {
          vapply(
            X = seq_along(x),
            FUN = function(x, y, k) {
              max <- length(y)
              lower <- x - k
              lower[lower < 1] <- 1
              upper <- x + k
              upper[upper > max] <- max
              any(y[lower:upper], na.rm = TRUE)
            },
            FUN.VALUE = logical(1),
            y = x,
            k = half_window
          )
        },
        half_window = half_window
      )
    }

    ## Build array
    select <- results$p.value <= alpha
    x <- y <- z <- object
    x[, !select] <- NA
    y[, select] <- NA
    z[!fit_roll] <- NA

    res <- array(
      data = c(z, x, y),
      dim = c(dim(object), 3),
      dimnames = c(dimnames(object), list(c("roll", "selection", "neutral")))
    )

    ts <- aion::series(object = res, time = dates)
    .IncrementTest(
      ts,
      statistic = results$t,
      parameter = 1L,
      p_value = results$p.value
    )
  }
)

#' @param x A \eqn{m x p} [`numeric`] matrix.
#' @param time A length-\eqn{m} [`numeric`] vector.
#' @param roll A [`logical`] scalar.
#' @param window An [`integer`].
#' @return A list of \eqn{p x 3} [`data.frame`].
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
test_FIT <- function(x, time, roll = FALSE, window = 3, ...) {
  ## Prepare data
  time <- as.numeric(time)

  ## Compute frequency as ratio of count of type of interest to all other types
  count_others <- lapply(
    X = seq_len(ncol(x)),
    FUN = function(i, data) { rowSums(data[, -i]) },
    data = x
  )
  freq <- x / do.call(cbind, count_others)

  ## Compute test
  if (!roll) {
    roll_list <- list(seq_len(nrow(freq)))
  } else {
    roll_index <- roll(freq, window = window)
    roll_list <- split(x = roll_index[["i"]], f = roll_index[["w"]])
  }

  k <- 1
  results <- vector(mode = "list", length = length(roll_list))
  for (i in roll_list) {
    roll_freq <- freq[i, ]
    roll_time <- time[i]
    roll_date <- roll_time[ceiling(window / 2)]

    fit <- apply(
      X = roll_freq,
      MARGIN = 2,
      FUN = function(v, t, ...) {
        tryCatch(
          error = function(cnd) c(t = NA, p.value = NA),
          { FIT(v, t, ...) }
        )
      },
      t = roll_time,
      ...
    )

    results[[k]] <- data.frame(
      column = colnames(fit),
      time = roll_date,
      t(fit),
      row.names = NULL
    )
    k <- k + 1
  }

  return(results)
}

#' Frequency Increment Test (FIT)
#'
#' @param v A [`numeric`] vector of frequencies.
#' @param t A [`numeric`] vector of time coordinates.
#' @param ... Extra parameter passed to [stats::t.test()].
#' @return A [`numeric`] vector containing the following components:
#'  \describe{
#'   \item{`t`}{the value of the test statistic.}
#'   \item{`p.value`}{the p-value for the test.}
#'  }
#' @author N. Frerebeau
#' @references
#'  Feder, A. F., Kryazhimskiy, S. & Plotkin, J. B. (2014). Identifying
#'  Signatures of Selection in Genetic Time Series. *Genetics*, 196(2),
#'  509-522. \doi{10.1534/genetics.113.158220}.
#' @keywords internal
#' @noRd
FIT <- function(v, t, ...) {
  ## Validation
  arkhe::assert_type(v, "numeric")
  arkhe::assert_type(t, "numeric")
  arkhe::assert_length(t, length(v))

  ## Remove zeros
  index <- v > 0
  v <- v[index]
  t <- t[index]
  if (length(t) < 3) {
    stop(tr_("A minimum of three time points is needed."), call. = FALSE)
  }

  ## Rescaled allele-frequency increments
  Y <- diff(v, lag = 1) /
    sqrt(2 * utils::head(v, -1) * (1 - utils::head(v, -1)) * diff(t, lag = 1))

  ## Statistics
  t_test <- stats::t.test(Y, ...)
  c(t = as.numeric(t_test$statistic), p.value = t_test$p.value)
}

# Plot =========================================================================
#' @export
#' @method plot IncrementTest
plot.IncrementTest <- function(x, calendar = get_calendar(),
                               col.neutral = "#004488", col.selection = "#BB5566",
                               col.roll = "grey", flip = FALSE, ncol = NULL,
                               xlab = NULL, ylab = NULL,
                               main = NULL, sub = NULL,
                               ann = graphics::par("ann"), axes = TRUE,
                               frame.plot = axes, ...) {
  ## Panel
  panel_fit <- function(x, y, ...) {
    graphics::lines(x, y, ...)
    graphics::points(x, y, ...)
  }

  ## Method for TimeSeries
  methods::callNextMethod(x = x, calendar = calendar, panel = panel_fit,
                          flip = flip, ncol = ncol,
                          xlab = xlab, ylab = ylab,
                          main = main, sub = sub,
                          ann = ann, axes = axes,
                          frame.plot = frame.plot,
                          col = c(col.roll, col.selection, col.neutral),
                          lwd = c(10, 1, 1), pch = 16, cex = 1, ...)
}

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

#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname seriate_average
#' @aliases seriate_average,data.frame-method
setMethod(
  f = "seriate_average",
  signature = c(object = "data.frame"),
  definition = function(object, margin = c(1, 2), axes = 1,
                        sup_row = NULL, sup_col = NULL, ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, margin = margin, axes = axes,
                         sup_row = sup_row, sup_col = sup_col, ...)
  }
)

#' @export
#' @rdname seriate_average
#' @aliases seriate_average,matrix-method
setMethod(
  f = "seriate_average",
  signature = c(object = "matrix"),
  definition = function(object, margin = c(1, 2), axes = 1,
                        sup_row = NULL, sup_col = NULL, ...) {
    ## Correspondence analysis
    corresp <- dimensio::ca(object, sup_row = sup_row, sup_col = sup_col)

    ## New PermutationOrder object
    as_seriation(corresp, margin = margin, axes = axes)
  }
)

# EVENT MODEL
#' @include AllGenerics.R
NULL

#' @export
#' @method summary EventDate
summary.EventDate <- function(object, ...) {
  summary(object@model, ...)
}

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

#' @method coef EventDate
#' @export
coef.EventDate <- function(object, calendar = NULL, ...) {
  z <- stats::coef(object@model, ...)
  if (is.null(calendar)) return(z)
  unclass(z) / aion::calendar_year(calendar) # Approximate
}

#' @export
#' @rdname model_event
#' @aliases coef,EventDate-method
setMethod("coef", "EventDate", coef.EventDate)

#' @method fitted EventDate
#' @export
fitted.EventDate <- function(object, calendar = NULL, ...) {
  z <- stats::fitted(object@model, ...)
  if (is.null(calendar)) return(z)
  aion::as_year(z, calendar = calendar)
}

#' @export
#' @rdname model_event
#' @aliases fitted,EventDate-method
setMethod("fitted", "EventDate", fitted.EventDate)

#' @method residuals EventDate
#' @export
residuals.EventDate <- function(object, calendar = NULL, ...) {
  z <- stats::residuals(object@model, ...)
  if (is.null(calendar)) return(z)
  unclass(z) / aion::calendar_year(calendar) # Approximate
}

#' @export
#' @rdname model_event
#' @aliases residuals,EventDate-method
setMethod("residuals", "EventDate", residuals.EventDate)

#' @method sigma EventDate
#' @export
sigma.EventDate <- function(object, calendar = NULL, ...) {
  z <- stats::sigma(object@model, ...)
  if (is.null(calendar)) return(z)
  unclass(z) / aion::calendar_year(calendar) # Approximate
}

#' @export
#' @rdname model_event
#' @aliases sigma,EventDate-method
setMethod("sigma", "EventDate", sigma.EventDate)

#' @method terms EventDate
#' @export
terms.EventDate <- function(x, ...) stats::terms(x@model, ...)

#' @export
#' @rdname model_event
#' @aliases terms,EventDate-method
setMethod("terms", "EventDate", terms.EventDate)

# SHOW METHODS
#' @include AllClasses.R
NULL

# EventDate ====================================================================
setMethod(
  f = "show",
  signature = "EventDate",
  definition = function(object) {
    show_ca <- utils::capture.output(methods::callNextMethod(object))
    summary_lm <- summary(object)
    sig <- sigma(object, calendar = get_calendar())

    cat(
      show_ca,
      "",
      paste0(tr_("Multiple Linear Regression"), ":"),
      sprintf(tr_("* R-squared: %s"), round(summary_lm$r.squared, 3)),
      sprintf(tr_("* Adjusted R-squared: %s"), round(summary_lm$adj.r.squared, 3)),
      sprintf(tr_("* Residual standard error: %s"), round(sig, 2)),
      sep = "\n"
    )
  }
)

# PermutationOrder =============================================================
setMethod(
  f = "show",
  signature = "PermutationOrder",
  definition = function(object) {
    k <- 50
    rows <- strtrim(paste0(object@rows_order, collapse = " "), k)
    columns <- strtrim(paste0(object@columns_order, collapse = " "), k)
    cat(
      paste0(tr_("Permutation order for matrix seriation"), ":"),
      sprintf(tr_("* Row order: %s"), paste0(rows, "...")),
      sprintf(tr_("* Column order: %s"), paste0(columns, "...")),
      sep = "\n"
    )
  }
)

# RefineCA =====================================================================
setMethod(
  f = "show",
  signature = "RefinePermutationOrder",
  definition = function(object) {
    value <- switch (object@margin, `1` = "Rows", `2` = "Columns")
    keep <- length(object@keep)
    total <- length(object@length)
    pc <- round(keep * 100 / total)
    methods::callNextMethod(object)
    cat(
      paste0("Partial bootstrap refinement", ":"),
      sprintf("* Cutoff value: %s", round(object@cutoff, digits = 2)),
      sprintf("* %s to keep: %d of %d (%g%%)", value, keep, total, pc),
      sep = "\n")
  }
)

# DATE MODEL
#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname event
#' @aliases event,data.frame,numeric-method
setMethod(
  f = "event",
  signature = c(object = "data.frame", dates = "numeric"),
  definition = function(object, dates, rank = NULL, sup_row = NULL,
                        calendar = CE(), ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, dates, rank = rank, sup_row = sup_row,
                         calendar = calendar, ...)
  }
)

#' @export
#' @rdname event
#' @aliases event,matrix,numeric-method
setMethod(
  f = "event",
  signature = c(object = "matrix", dates = "numeric"),
  definition = function(object, dates, calendar = CE(),
                        rank = NULL, sup_row = NULL, total = 5,
                        verbose = getOption("kairos.verbose"), ...) {
    ## Sample
    n <- nrow(object)

    ## Check dates
    if (arkhe::has_names(dates)) {
      i <- match(names(dates), rownames(object))
      old_dates <- dates
      dates <- rep(NA_real_, n)
      dates[i] <- old_dates
    }
    arkhe::assert_length(dates, nrow(object))

    ## Supplementary rows
    sup <- logical(n)
    sup[sup_row] <- TRUE

    data_ref <- object[!sup, , drop = FALSE]
    data_sup <- object[sup, , drop = FALSE]
    dates_ref <- dates[!sup]
    dates_sup <- dates[sup]

    ## Cleansing
    data_clean <- clean_total(data_ref, data_sup, dates_ref, dates_sup,
                              total = total, verbose = verbose)

    data <- rbind(data_clean$data_ref, data_clean$data_sup)
    dates <- c(data_clean$dates_ref, data_clean$dates_sup)
    sup <- seq_along(data_clean$dates_sup) + length(data_clean$dates_ref)

    ## Correspondance analysis
    if (is.null(rank)) {
      tmp <- dimensio::ca(data, rank = NULL, sup_row = sup, ...)
      eig <- dimensio::get_eigenvalues(tmp)
      rank <- which.max(eig$cumulative >= 60)
    }
    rank <- min(rank, dim(data) - 1)
    results_CA <- dimensio::ca(data, rank = rank, sup_row = sup, ...)

    ## Get row coordinates
    row_coord <- dimensio::get_coordinates(results_CA, margin = 1)

    ## Rata die
    if (is.null(calendar)) {
      rd <- aion::as_fixed(data_clean$dates_ref)
    } else {
      rd <- aion::fixed(data_clean$dates_ref, calendar = calendar)
    }

    ## Gaussian multiple linear regression model
    contexts <- data.frame(
      date = rd,
      row_coord[!row_coord$.sup, -ncol(row_coord), drop = FALSE]
    )
    fit <- stats::lm(date ~ ., data = contexts)

    ## FIXME: temporary workaround until aion::fixed() can deal with NA
    dates_not_na <- which(!is.na(dates))
    rt <- rep(NA_real_, length(dates))
    if (is.null(calendar)) {
      rt[dates_not_na] <- aion::as_fixed(dates[dates_not_na])
    } else {
      rt[dates_not_na] <- aion::fixed(dates[dates_not_na], calendar = calendar)
    }
    rt <- aion::as_fixed(rt)

    .EventDate(
      results_CA,
      dates = rt,
      model = fit,
      keep = seq_len(rank)
    )
  }
)

clean_total <- function(data_ref, data_sup, dates_ref, dates_sup,
                        total = 5, verbose = TRUE) {

  clean <- TRUE
  n_rm_col <- n_rm_ref <- n_rm_sup <- 0

  while (clean) {
    rm_col <- colSums(data_ref) < total
    n_rm_col <- n_rm_col + sum(rm_col)
    if (any(rm_col)) {
      data_ref <- data_ref[, !rm_col, drop = FALSE]
      data_sup <- data_sup[, !rm_col, drop = FALSE]
    }

    rm_row_ref <- rowSums(data_ref) < total
    n_rm_ref <- n_rm_ref + sum(rm_row_ref)
    if (any(rm_row_ref)) {
      data_ref <- data_ref[!rm_row_ref, , drop = FALSE]
      dates_ref <- dates_ref[!rm_row_ref]
    }

    rm_row_sup <- rowSums(data_sup) < total
    n_rm_sup <- n_rm_sup + sum(rm_row_sup)
    if (any(data_sup)) {
      data_sup <- data_sup[!rm_row_sup, , drop = FALSE]
      dates_sup <- dates_sup[!rm_row_sup]
    }

    if (!any(rm_col) & !any(rm_row_ref) & !any(rm_row_sup)) clean <- FALSE
  }

  if (isTRUE(verbose)) {
    if (n_rm_col > 0) {
      msg <- ngettext(n_rm_col,
                      "%d column has a grand total of less than %d.",
                      "%d columns have a grand total of less than %d.")
      message(sprintf(msg, n_rm_col, total))
    }
    if (n_rm_ref > 0) {
      msg <- ngettext(n_rm_ref,
                      "%d row has a grand total of less than %d.",
                      "%d rows have a grand total of less than %d.")
      message(sprintf(msg, n_rm_ref, total))
    }
    if (n_rm_sup > 0) {
      msg <- ngettext(n_rm_sup,
                      "%d supplementary row has a grand total of less than %d.",
                      "%d supplementary rows have a grand total of less than %d.")
      message(sprintf(msg, n_rm_sup, total))
    }
    n_rm <- n_rm_col + n_rm_ref + n_rm_sup
    if (n_rm > 0) {
      msg <- ngettext(n_rm, "It was omitted from the analysis.",
                      "They were omitted from the analysis.")
      message(msg)
    }
  }

  if (all(is.na(dates_ref))) {
    stop(tr_("All dates are missing!"), call. = FALSE)
  }

  list(
    data_ref = data_ref,
    data_sup = data_sup,
    dates_ref = dates_ref,
    dates_sup = dates_sup
  )
}

# Event ========================================================================
#' @export
#' @rdname predict_event
#' @aliases predict_event,EventDate,missing-method
setMethod(
  f = "predict_event",
  signature = c(object = "EventDate", data = "missing"),
  definition = function(object, margin = 1, level = 0.95,
                        calendar = get_calendar()) {
    data <- dimensio::get_coordinates(object, margin = margin)
    data <- as.matrix(data[, -ncol(data), drop = FALSE])

    .predict_event(object, data, level = level, calendar = calendar)
  }
)

#' @export
#' @rdname predict_event
#' @aliases predict_event,EventDate,matrix-method
setMethod(
  f = "predict_event",
  signature = c(object = "EventDate", data = "matrix"),
  definition = function(object, data, margin = 1, level = 0.95,
                        calendar = get_calendar()) {
    ## Correspondence analysis
    data <- dimensio::predict(object, data, margin = margin)

    .predict_event(object, data, level = level, calendar = calendar)
  }
)

.predict_event <- function(object, data, level = 0.95,
                           calendar = get_calendar()) {
  ## Predict event date
  fit_model <- object@model
  ca_event <- compute_event(fit_model, data, level = level)

  # TODO: check predicted dates consistency

  ## Convert
  ca_event <- as.data.frame(ca_event)
  if (is.null(calendar)) return(ca_event)

  ca_event[] <- lapply(
    X = ca_event,
    FUN = aion::as_year,
    calendar = calendar,
    decimal = TRUE
  )
  ca_event
}

#' @export
#' @rdname density_event
#' @aliases density_event,EventDate-method
setMethod(
  f = "density_event",
  signature = c(object = "EventDate"),
  definition = function(object, dates = NULL, calendar = NULL, n = 500, ...) {
    ## Get data
    rows <- predict_event(object, margin = 1, calendar = NULL)
    row_dates <- rows$date
    row_errors <- rows$error

    if (is.null(dates)) {
      dates <- seq(
        from = min(rows$lower, na.rm = TRUE) * 0.96,
        to = max(rows$upper, na.rm = TRUE) * 1.04,
        length.out = n
      )
    }
    if (is.null(calendar)) {
      dates <- aion::as_fixed(dates)
    } else {
      dates <- aion::fixed(dates, calendar = calendar)
    }

    ## Density estimation
    date_event <- mapply(
      FUN = stats::dnorm,
      mean = row_dates,
      sd = row_errors,
      MoreArgs = list(x = dates),
      SIMPLIFY = TRUE
    )

    colnames(date_event) <- rownames(rows)
    aion::series(object = date_event, time = dates)
  }
)

# Accumulation =================================================================
#' @export
#' @rdname predict_event
#' @aliases predict_accumulation,EventDate,missing-method
setMethod(
  f = "predict_accumulation",
  signature = c(object = "EventDate", data = "missing"),
  definition = function(object, level = 0.95,
                        calendar = get_calendar()) {
    data <- object@data
    methods::callGeneric(object, data = data, level = level, calendar = calendar)
  }
)

#' @export
#' @rdname predict_event
#' @aliases predict_accumulation,EventDate,matrix-method
setMethod(
  f = "predict_accumulation",
  signature = c(object = "EventDate", data = "matrix"),
  definition = function(object, data, level = 0.95,
                        calendar = get_calendar()) {
    ## Predict event date
    col_event <- predict_event(object, margin = 2, calendar = NULL)

    ## Accumulation time point estimate
    date_range <- seq(
      from = min(col_event$lower, na.rm = TRUE),
      to = max(col_event$upper, na.rm = TRUE),
      length.out = 1000
    )
    dens <- mapply(
      FUN = stats::dnorm,
      mean = col_event$date,
      sd = col_event$error,
      MoreArgs = list(date_range),
      SIMPLIFY = TRUE
    )
    acc <- apply(
      X = data,
      MARGIN = 1,
      FUN = function(x, density) {
        colSums(t(density) * as.numeric(x))
      },
      density = dens
    )

    alpha <- (1 - level) / 2
    probs <- c(0.5, alpha, 1 - alpha)
    quant <- quantile_density(x = date_range, y = acc, probs = probs)

    ## Convert
    quant <- as.data.frame(quant)
    colnames(quant) <- c("date", "lower", "upper")
    if (is.null(calendar)) return(quant)

    quant[] <- lapply(
      X = quant,
      FUN = aion::as_year,
      calendar = calendar,
      decimal = TRUE
    )
    quant
  }
)

#' @export
#' @rdname density_event
#' @aliases density_accumulation,EventDate-method
setMethod(
  f = "density_accumulation",
  signature = c(object = "EventDate"),
  definition = function(object, dates = NULL, calendar = NULL,
                        type = c("activity", "tempo"), n = 500, ...) {
    ## Validation
    type <- match.arg(type, several.ok = FALSE)

    ## Get data
    columns <- predict_event(object, margin = 2, calendar = NULL)
    col_dates <- columns$date
    col_errors <- columns$error

    if (is.null(dates)) {
      dates <- seq(
        from = min(columns$lower, na.rm = TRUE) * 0.96,
        to = max(columns$upper, na.rm = TRUE) * 1.04,
        length.out = n
      )
    }
    if (is.null(calendar)) {
      dates <- aion::as_fixed(dates)
    } else {
      dates <- aion::fixed(dates, calendar = calendar)
    }

    ## Weighted sum of the fabric dates
    counts <- dimensio::get_data(object)
    freq <- counts / rowSums(counts)

    ## Tempo vs activity plot
    fun <- switch(
      type,
      activity = stats::dnorm,
      tempo = stats::pnorm
    )

    ## Density estimation
    col_density <- mapply(
      FUN = fun,
      mean = col_dates,
      sd = col_errors,
      MoreArgs = list(dates),
      SIMPLIFY = TRUE
    )

    date_acc <- apply(
      X = freq,
      MARGIN = 1,
      FUN = function(x, density) {
        colSums(t(density) * as.numeric(x))
      },
      density = col_density
    )

    aion::series(object = date_acc, time = dates)
  }
)

# Helpers ======================================================================
#' Predict event dates
#'
#' @param fit A \code{\link[stats:lm]{multiple linear model}}.
#' @param data A [`numeric`] matrix giving the coordinates in CA space.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @return
#'  A four columns [`numeric`] matrix giving the predicted event dates,
#'  the corresponding confidence interval boundaries and the standard error of
#'  the predicted dates.
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
compute_event <- function(fit, data, level, error = TRUE) {
  data <- as.data.frame(data)

  date_predict <- stats::predict.lm(
    object = fit,
    newdata = data,
    se.fit = TRUE,
    interval = "confidence",
    level = level
  )

  # Three columns matrix: predicted value + CI boundaries
  results <- as.data.frame(date_predict$fit)
  rownames(results) <- rownames(data)
  colnames(results) <- c("date", "lower", "upper")
  if (error) results$error <- date_predict$se.fit

  results
}

# REFINE DATE MODEL
#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname resample_event
#' @aliases jackknife,EventDate-method
setMethod(
  f = "jackknife",
  signature = c(object = "EventDate"),
  definition = function(object, level = 0.95,
                        calendar = get_calendar(),
                        progress = getOption("kairos.progress"),
                        verbose = getOption("kairos.verbose"), ...) {
    ## Get data
    fit_model <- object@model
    fit_dates <- time(object)
    fit_data <- dimensio::get_data(object)
    fit_dim <- object@keep

    ## Compute event dates
    event <- predict_event(object, level = level, calendar = calendar)

    ## Jackknife model coefficients
    jack_values <- compute_date_jack(
      x = fit_data,
      dates = fit_dates,
      rank = length(fit_dim),
      progress = progress,
      verbose = verbose
    )
    jack_coef <- colMeans(jack_values)

    ## Change lm coefficients (UGLY)
    fit_model$coefficients <- jack_coef[c(1, fit_dim + 1)]

    ## Predict event date for each context
    row_coord <- dimensio::get_coordinates(object, margin = 1)
    jack_event <- compute_event(fit_model, row_coord, level = level, error = FALSE)

    results <- as.data.frame(jack_event)
    if (!is.null(calendar)) {
      results[] <- lapply(
        X = results,
        FUN = aion::as_year,
        calendar = calendar,
        decimal = TRUE
      )
    }
    results
  }
)

#' Jackknife Fabrics
#'
#' Compute date event jackknifed statistics for each replicated sample.
#' @param x A [`numeric`] matrix of count data.
#' @param dates A [`numeric`] vector of known dates.
#' @param rank A [`numeric`] value.
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @param ... Currently not used].
#' @return A `numeric` [`matrix`] of linear model coefficients.
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
compute_date_jack <- function(x, dates, rank = 10,
                              progress = FALSE, verbose = FALSE, ...) {
  m <- ncol(x)
  k <- seq_len(m)
  jack <- vector(mode = "list", length = m)

  progress_bar <- interactive() && isTRUE(progress)
  if (progress_bar) pbar <- utils::txtProgressBar(max = m, style = 3)

  for (j in k) {
    counts <- x[, -j, drop = FALSE]
    model <- event(counts, dates = dates, calendar = NULL, rank = rank, verbose = verbose)
    jack[[j]] <- coef(model, calendar = NULL) # Get model coefficients
    if (progress_bar) utils::setTxtProgressBar(pbar, j)
  }

  if (progress_bar) close(pbar)
  do.call(rbind, jack)
}

#' @export
#' @rdname resample_event
#' @aliases bootstrap,EventDate-method
setMethod(
  f = "bootstrap",
  signature = c(object = "EventDate"),
  definition = function(object, level = 0.95, probs = c(0.05, 0.95), n = 1000,
                        calendar = get_calendar(),
                        progress = getOption("kairos.progress"), ...) {
    ## Get data
    fit_model <- object@model
    fit_dim <- object@keep

    ## Partial bootstrap CA
    ## /!\ Be careful: EventDate inherits from CA
    ca_boot <- methods::callNextMethod(object, n = n)
    ca_rep_row <- dimensio::get_replications(ca_boot, margin = 1)

    ## Bootstrap assemblages
    event_rows <- apply(
      X = ca_rep_row,
      MARGIN = 1,
      FUN = function(x, axes, model, level, probs) {
        compute_date_boot(t(x), axes, model, level, probs)
      },
      axes = fit_dim,
      model = fit_model,
      level = level,
      probs = probs
    )

    results <- as.data.frame(t(event_rows))
    if (is.null(calendar)) return(results)

    results[] <- lapply(
      X = results,
      FUN = aion::as_year,
      calendar = calendar,
      decimal = TRUE
    )
    results
  }
)

#' Bootstrap Resampling of Assemblages
#'
#' Computes date event bootstraped statistics for each replicated sample.
#' @param x A [`numeric`] matrix of bootstrap replicates.
#' @param axes A [`numeric`] vector giving the subscripts of the CA components
#'  to keep.
#' @param model An object of class [`lm`][stats::lm()].
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @param probs A [`numeric`] vector of probabilities with values in \eqn{[0,1]}
#'  (see [stats:quantile()]).
#' @return A [`numeric`] vector with the following elements:
#'  \describe{
#'   \item{`min`}{Minimum value.}
#'   \item{`mean`}{Mean value (event date).}
#'   \item{`max`}{Maximum value.}
#'   \item{`Q05`}{Sample quantile to 0.05 probability.}
#'   \item{`Q95`}{Sample quantile to 0.95 probability.}
#'  }
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
compute_date_boot <- function(x, axes, model, level, probs = c(0.05, 0.95)) {
  ## Remove missing values
  coords <- stats::na.omit(x[, axes])

  ## Gaussian multiple linear regression model
  event <- compute_event(model, coords, level)[, 1]
  Q <- stats::quantile(event, probs = probs, names = FALSE)

  distrib <- c(min(event), mean(event), max(event), Q)
  quant <- paste0("Q", round(probs * 100, 0))
  names(distrib) <- c("min", "mean", "max", quant)
  return(distrib)
}

# PLOT DATES
#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @method plot EventDate
plot.EventDate <- function(x, type = c("activity", "tempo"), event = FALSE,
                           calendar = get_calendar(),
                           select = 1, n = 500, eps = 1e-09,
                           col.accumulation = "black", col.event = "red",
                           flip = FALSE, ncol = NULL,
                           xlab = NULL, ylab = NULL,
                           main = NULL, sub = NULL,
                           ann = graphics::par("ann"), axes = TRUE,
                           frame.plot = axes, ...) {
  ## Validation
  type <- match.arg(type, several.ok = FALSE)

  ## Event date
  date_event <- density_event(x, n = n)
  date_range <- aion::time(date_event)
  cases <- colnames(date_event)

  ## Accumulation time
  date_acc <- density_accumulation(x, dates = date_range, type = type, n = n)

  ## Selection
  if (is.null(select)) index <- seq_along(cases)
  else if (is.character(select)) index <- which(cases %in% select)
  else index <- as.numeric(select)

  k <- length(index)
  if (k == 0) stop(tr_("Wrong selection."), call. = FALSE)

  if (type != "activity" || !event) {
    date_event <- array(data = NA, dim = list(n, k, 1))
  }

  date_event <- date_event[, index, , drop = FALSE]
  date_acc <- date_acc[, index, , drop = FALSE]

  ## Cleaning
  # date_acc[date_acc < eps] <- NA
  date_event[date_event < eps] <- NA
  na_event <- apply(date_event, 1, function(x) all(is.na(x)))
  na_acc <- apply(date_acc, 1, function(x) all(is.na(x)))
  date_drop <- na_event & na_acc

  ## Time series
  ts <- array(data = c(date_acc, date_event), dim = c(n, k, 2))
  dimnames(ts) <- list(NULL, cases[index], c("accumulation", "event"))
  ts <- aion::series(object = ts[!date_drop, , , drop = FALSE],
                     time = aion::as_fixed(date_range[!date_drop]))

  panel <- switch(
    type,
    activity = 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", lty = list(...)$lty)
    },
    tempo = function(x, y, ...) graphics::lines(x, y, col = "black", lty = 1)
  )

  col <- c(col.accumulation, grDevices::adjustcolor(col.event, alpha = 0.5))
  aion::plot(ts, panel = panel, calendar = calendar,
             flip = flip, ncol = ncol, xlab = xlab, ylab = ylab,
             main = main, sub = sub, ann = ann, axes = axes,
             frame.plot = frame.plot,
             col = col, lty = c(0, 0))

  invisible(x)
}

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

#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname as_seriation
#' @aliases as_seriation,CA-method
setMethod(
  f = "as_seriation",
  signature = c(object = "CA"),
  definition = function(object, margin = c(1, 2), axes = 1) {
    ## Validation
    arkhe::assert_length(axes, 1)

    margin <- as.integer(margin)
    axes <- as.integer(axes)

    ## Original sequences
    data <- dimensio::get_data(object)
    i <- seq_len(nrow(data))
    j <- seq_len(ncol(data))

    ## Correspondence analysis
    rows <- dimensio::get_coordinates(object, margin = 1)
    cols <- dimensio::get_coordinates(object, margin = 2)

    ## Reorder in case if supplementary observations
    rows <- rows[order(object@rows@order), ]
    cols <- cols[order(object@columns@order), ]

    ## Seriation order
    row_coords <- if (any(margin == 1)) order(rows[, axes]) else i
    col_coords <- if (any(margin == 2)) order(cols[, axes]) else j

    ## New PermutationOrder object
    .AveragePermutationOrder(
      object,
      rows_order = as.integer(row_coords),
      columns_order = as.integer(col_coords)
    )
  }
)

# MUTATORS
#' @include AllClasses.R
NULL

# Getters ======================================================================
#' @export
#' @rdname series
#' @aliases time,EventDate-method
setMethod(
  f = "time",
  signature = "EventDate",
  definition = function(x, calendar = NULL) {
    z <- x@dates
    if (is.null(calendar)) return(z)
    aion::as_year(z, calendar = calendar, decimal = TRUE)
  }
)

#' @export
#' @rdname series
#' @aliases span,AoristicSum-method
setMethod(
  f = "span",
  signature = "AoristicSum",
  definition = function(x, calendar = NULL) {
    z <- x@span
    if (is.null(calendar)) return(unclass(z))
    aion::as_year(z, calendar = calendar, shift = FALSE)
  }
)

#' @export
#' @rdname mutators
#' @aliases weights,AoristicSum-method
setMethod("weights", "AoristicSum", function(object, ...) object@p)

#' @export
#' @rdname mutators
#' @aliases weights,CountApportion-method
setMethod("weights", "CountApportion", function(object, ...) object@p)

# Setters ======================================================================

# GENERIC METHODS
#' @include AllClasses.R
NULL

# Import S4 generics ===========================================================
#' @importMethodsFrom arkhe jackknife
#' @importMethodsFrom arkhe bootstrap
NULL

# Mutators =====================================================================
## 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

## Extract ---------------------------------------------------------------------
#' Get or Set Parts of an Object
#'
#' Getters and setters to retrieve or set parts of an object.
#' @param object An object from which to get or set element(s).
# @param value A possible value for the element(s) of `x`.
#' @param ... Currently not used.
# @return
#  * `set_*()` returns an object of the same sort as `x` with the new values
#    assigned.
#  * `get_*()` returns the part of `x`.
# @example inst/examples/ex-mutators.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name mutators
#' @rdname mutators
#' @aliases get set
NULL

## Time series -----------------------------------------------------------------
#' Sampling Times
#'
#' Get the times at which a time series was sampled.
#' @param x An \R object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL` (the default), *rata die* are
#'  returned.
#' @return
#'  A [`numeric`] vector.
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name series
#' @rdname series
NULL

## Coerce ----------------------------------------------------------------------
#' Coerce to a Data Frame
#'
#' @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 row.names,optional Currently not used.
#' @param ... Further parameters to be passed to [data.frame()].
#' @return
#'  A [`data.frame`] with an extra `time` column giving the (decimal) years at
#'  which the time series was sampled.
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name data.frame
#' @rdname data.frame
NULL

# Mean Ceramic Dates ===========================================================
#' Mean Ceramic Date
#'
#' Estimates the Mean Ceramic Date of an assemblage.
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data (absolute frequencies giving the number of
#'  individuals for each category, i.e. a contingency table). A [`data.frame`]
#'  will be coerced to a `numeric` `matrix` via [data.matrix()].
#' @param dates A length-\eqn{p} [`numeric`] vector of dates expressed in years.
#' @param calendar An [`aion::TimeScale-class`] object specifying the calendar
#'  of `dates` (see [aion::calendar()]). Defaults to Gregorian Common Era.
#' @param ... Currently not used.
#' @details
#'  The Mean Ceramic Date (MCD) is a point estimate of the occupation of an
#'  archaeological site (South 1977). The MCD is estimated as the weighted mean
#'  of the date midpoints of the ceramic types (based on absolute dates or the
#'  known production interval) found in a given assemblage. The weights are the
#'  relative frequencies of the respective types in the assemblage.
#'
#'  A bootstrapping procedure is used to estimate the confidence interval of a
#'  given MCD. For each assemblage, a large number of new bootstrap replicates
#'  is created, with the same sample size, by resampling the original
#'  assemblage with replacement. MCDs are calculated for each replicates and
#'  upper and lower boundaries of the confidence interval associated with each
#'  MCD are then returned.
#' @return
#'  A [`MeanDate-class`] object.
#' @references
#'  South, S. A. (1977). *Method and Theory in Historical Archaeology*.
#'  New York: Academic Press.
#' @example inst/examples/ex-mcd.R
#' @author N. Frerebeau
#' @family mean ceramic date tools
#' @family dating methods
#' @docType methods
#' @rdname mcd
#' @aliases mcd-method
setGeneric(
  name = "mcd",
  def = function(object, dates, ...) standardGeneric("mcd")
)

#' Resample Mean Ceramic Dates
#'
#' @description
#'  * `bootstrap()` generate bootstrap estimations of an [MCD][mcd()].
#'  * `jackknife()` generate jackknife estimations of an [MCD][mcd()].
#' @param object A [`MeanDate-class`] object (typically returned by [mcd()]).
#' @param n A non-negative [`integer`] specifying the number of bootstrap
#'  replications.
#' @param nsim A non-negative [`integer`] specifying the number of simulations.
#' @param f A [`function`] that takes a single numeric vector (the result of
#'  the resampling procedure) as argument.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]).
#' @return
#'  If `f` is `NULL`, `bootstrap()` and `jackknife()` return a [`data.frame`]
#'  with the following elements (else, returns the result of `f` applied to the
#'  `n` resampled values) :
#'  \describe{
#'   \item{original}{The observed value.}
#'   \item{mean}{The bootstrap/jackknife estimate of mean.}
#'   \item{bias}{The bootstrap/jackknife estimate of bias.}
#'   \item{error}{The boostrap/jackknife estimate of standard erro.}
#'  }
#' @author N. Frerebeau
#' @docType methods
#' @family mean ceramic date tools
#' @name resample_mcd
#' @rdname resample_mcd
NULL

#' MCD Plot
#'
#' @param x A [`MeanDate-class`] object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]).
#' @param interval A [`character`] string giving the type of confidence
#'  interval to be returned. It must be one "`student`" (the default),
#'  "`normal`", "`percentiles`" or "`range`" (min-max).
#'  Any unambiguous substring can be given.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#'  Only used if `interval` is not "`range`".
#' @param decreasing A [`logical`] scalar: should the sort be increasing or
#'  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,
#'  y and z axis 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 ... Further [graphical parameters][graphics::par].
#' @return
#'  `plot()` is called it for its side-effects: it results in a graphic being
#'  displayed (invisibly returns `x`).
#' @example inst/examples/ex-mcd.R
#' @author N. Frerebeau
#' @family mean ceramic date tools
#' @docType methods
#' @name plot_mcd
#' @rdname plot_mcd
NULL

# Event Dates ==================================================================
#' Event and Accumulation Dates
#'
#' Fits a date event model.
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data (absolute frequencies giving the number of
#'  individuals for each category, i.e. a contingency table). A [`data.frame`]
#'  will be coerced to a `numeric` `matrix` via [data.matrix()].
#' @param dates A [`numeric`] vector of dates. If named, the names must match
#'  the row names of `object`.
#' @param calendar An [`aion::TimeScale-class`] object specifying the calendar
#'  of `dates` (see [aion::calendar()]). Defaults to Gregorian Common Era.
#' @param rank An [`integer`] specifying the number of CA factorial components
#'  to be use for linear model fitting (see details). If `NULL` (the default),
#'  axes corresponding to at least 60% of the inertia will be used.
#' @param sup_row A [`numeric`] or [`logical`] vector specifying the indices of
#'  the supplementary rows.
#' @param total A length-one [`numeric`] vector specifying the minimum total of
#'  a row/column. Rows/columns whose total is less than this value will be
#'  omitted from the analysis.
#' @param verbose A [`logical`] scalar: should \R report extra information
#'  on progress?
#' @param ... Further arguments to be passed to internal methods.
#' @details
#'  This is an implementation of the chronological modeling method proposed by
#'  Bellanger and Husi (2012, 2013).
#'
#'  Event and accumulation dates are density estimates of the occupation and
#'  duration of an archaeological site (Bellanger and Husi 2012, 2013).
#'  The event date is an estimation of the *terminus post-quem* of an
#'  archaeological assemblage. The accumulation date represents the
#'  "chronological profile" of the assemblage. According to Bellanger and Husi
#'  (2012), accumulation date can be interpreted "at best \[...\] as a formation
#'  process reflecting the duration or succession of events on the scale of
#'  archaeological time, and at worst, as imprecise dating due to contamination
#'  of the context by residual or intrusive material." In other words,
#'  accumulation dates estimate occurrence of archaeological events and rhythms
#'  of the long term.
#'
#'  Dates are converted to *[rata die][aion::RataDie-class]* before any
#'  computation.
#'
#'  This method relies on strong archaeological and statistical assumptions
#'  (see `vignette("event")`).
#' @return
#'  An [`EventDate-class`] object.
#' @seealso [`plot()`][plot_event], [`predict_event()`][predict_event],
#'  [`predict_accumulation()`][predict_event], [`jackknife()`][resample_event],
#'  [`bootstrap()`][resample_event]
#' @references
#'  Bellanger, L. & Husi, P. (2013). Mesurer et modéliser le temps inscrit dans
#'  la matière à partir d'une source matérielle : la céramique médiévale.
#'  In *Mesure et Histoire Médiévale*. Histoire ancienne et médiévale.
#'  Paris: Publication de la Sorbonne, p. 119-134.
#'
#'  Bellanger, L. & Husi, P. (2012). Statistical Tool for Dating and
#'  Interpreting Archaeological Contexts Using Pottery. *Journal of
#'  Archaeological Science*, 39(4), 777-790. \doi{10.1016/j.jas.2011.06.031}.
#'
#'  Bellanger, L., Tomassone, R. & Husi, P. (2008). A Statistical Approach for
#'  Dating Archaeological Contexts. *Journal of Data Science*, 6, 135-154.
#'
#'  Bellanger, L., Husi, P. & Tomassone, R. (2006). Une approche statistique
#'  pour la datation de contextes archéologiques. *Revue de Statistique
#'  Appliquée*, 54(2), 65-81.
#'
#'  Bellanger, L., Husi, P. & Tomassone, R. (2006). Statistical Aspects of
#'  Pottery Quantification for the Dating of Some Archaeological Contexts.
#'  *Archaeometry*, 48(1), 169-183. \doi{10.1111/j.1475-4754.2006.00249.x}.
#'
#'  Poblome, J. & Groenen, P. J. F. (2003). Constrained Correspondence Analysis
#'  for Seriation of Sagalassos Tablewares. In Doerr, M. & Apostolis, S. (eds.),
#'  *The Digital Heritage of Archaeology*. Athens: Hellenic Ministry of Culture.
#' @example inst/examples/ex-event.R
#' @author N. Frerebeau
#' @family event date tools
#' @family dating methods
#' @docType methods
#' @aliases event-method
setGeneric(
  name = "event",
  def = function(object, dates, ...) standardGeneric("event"),
  valueClass = "EventDate"
)

#' Predict Event and Accumulation Dates
#'
#' Estimates the event and accumulation dates of an assemblage.
#' @param object An [`EventDate-class`] object.
#' @param data A `numeric` [`matrix`] or a [`data.frame`] of count data
#' (absolute frequencies) for which to predict event and accumulation dates.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @param margin A [`numeric`] vector giving the subscripts which the prediction
#'  will be applied over: `1` indicates rows, `2` indicates columns.
#' @param ... Further arguments to be passed to internal methods.
#' @return
#'  A [`data.frame`].
#' @references
#'  Bellanger, L. & Husi, P. (2013). Mesurer et modéliser le temps inscrit dans
#'  la matière à partir d'une source matérielle : la céramique médiévale.
#'  In *Mesure et Histoire Médiévale*. Histoire ancienne et médiévale.
#'  Paris: Publication de la Sorbonne, p. 119-134.
#'
#'  Bellanger, L. & Husi, P. (2012). Statistical Tool for Dating and
#'  Interpreting Archaeological Contexts Using Pottery. *Journal of
#'  Archaeological Science*, 39(4), 777-790. \doi{10.1016/j.jas.2011.06.031}.
#'
#'  Bellanger, L., Tomassone, R. & Husi, P. (2008). A Statistical Approach for
#'  Dating Archaeological Contexts. *Journal of Data Science*, 6, 135-154.
#'
#'  Bellanger, L., Husi, P. & Tomassone, R. (2006). Une approche statistique
#'  pour la datation de contextes archéologiques. *Revue de Statistique
#'  Appliquée*, 54(2), 65-81.
#'
#'  Bellanger, L., Husi, P. & Tomassone, R. (2006). Statistical Aspects of
#'  Pottery Quantification for the Dating of Some Archaeological Contexts.
#'  *Archaeometry*, 48(1), 169-183. \doi{10.1111/j.1475-4754.2006.00249.x}.
#' @example inst/examples/ex-event.R
#' @author N. Frerebeau
#' @family event date tools
#' @docType methods
#' @aliases predict_event-method
setGeneric(
  name = "predict_event",
  def = function(object, data, ...) standardGeneric("predict_event"),
  valueClass = "data.frame"
)

#' @rdname predict_event
#' @aliases predict_accumulation-method
setGeneric(
  name = "predict_accumulation",
  def = function(object, data, ...) standardGeneric("predict_accumulation"),
  valueClass = "data.frame"
)

#' Density of Event and Accumulation Dates
#'
#' Estimates the event and accumulation density.
#' @param object An [`EventDate-class`] object.
#' @param dates A [`numeric`] vector of dates expressed as `calendar` years or
#'  *rata die* (if `calendar` is `NULL`).
#' @param calendar An [`aion::TimeScale-class`] object specifying the calendar
#'  of `dates` (see [aion::calendar()]). If `NULL` (the default), *rata die* are
#'  expected.
#' @param type A [`character`] string indicating the type of plot.
#'  It must be one of "`activity`" (default) or "`tempo`" (see details).
#'  Any unambiguous substring can be given.
#' @param n A length-one non-negative [`numeric`] vector giving the desired
#'  length of the vector of quantiles for density computation.
#' @param ... Currently not used.
#' @return
#'  An [`aion::TimeSeries-class`] object.
#' @example inst/examples/ex-event.R
#' @author N. Frerebeau
#' @family event date tools
#' @docType methods
#' @aliases density_event-method
setGeneric(
  name = "density_event",
  def = function(object, ...) standardGeneric("density_event")
)

#' @rdname density_event
#' @aliases density_accumulation-method
setGeneric(
  name = "density_accumulation",
  def = function(object, ...) standardGeneric("density_accumulation")
)

#' Resample Event Dates
#'
#' @description
#'  * `bootstrap()` generate bootstrap estimations of an [event][event()].
#'  * `jackknife()` generate jackknife estimations of an [event][event()].
#' @param object An [`EventDate-class`] object (typically returned by [event()]).
#' @param n A non-negative [`integer`] specifying the number of bootstrap
#'  replications.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @param probs A [`numeric`] vector of probabilities with values in
#'  \eqn{[0,1]}.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL`, *rata die* are returned.
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @param verbose A [`logical`] scalar: should \R report extra information
#'  on progress?
#' @param ... Further arguments to be passed to internal methods.
#' @details
#'  If `jackknife()` is used, one type/fabric is removed at a
#'  time and all statistics are recalculated. In this way, one can assess
#'  whether certain type/fabric has a substantial influence on the date
#'  estimate.
#'  A three columns `data.frame` is returned, giving the results of
#'  the resampling procedure (jackknifing fabrics) for each assemblage (in rows)
#'  with the following columns:
#'  \describe{
#'   \item{`mean`}{The jackknife mean (event date).}
#'   \item{`lower`}{The lower boundary of the confidence interval.}
#'   \item{`upper`}{The upper boundary of the confidence interval.}
#'  }
#'
#'  If `bootstrap()` is used, a large number of new bootstrap assemblages is
#'  created, with the same sample size, by resampling each of the original
#'  assemblage with replacement. Then, examination of the bootstrap statistics
#'  makes it possible to pinpoint assemblages that require further
#'  investigation.
#'
#'  A five columns `data.frame` is returned, giving the bootstrap
#'  distribution statistics for each replicated assemblage (in rows)
#'  with the following columns:
#'  \describe{
#'   \item{`min`}{Minimum value.}
#'   \item{`mean`}{Mean value (event date).}
#'   \item{`max`}{Maximum value.}
#'   \item{`Q5`}{Sample quantile to 0.05 probability.}
#'   \item{`Q95`}{Sample quantile to 0.95 probability.}
#'  }
#' @return
#'  A [`data.frame`].
#' @author N. Frerebeau
#' @docType methods
#' @family event date tools
#' @name resample_event
#' @rdname resample_event
NULL

#' Extract Event Date Model Results
#'
#' @description
#'  * `summary()` summarizes linear model fit.
#'  * `coef()` extracts model coefficients (see [stats::coef()]).
#'  * `fitted()` extracts model fitted values (see [stats::fitted()]).
#'  * `residuals()` extracts model residuals (see [stats::residuals()]).
#'  * `sigma()` extracts the residual standard deviation (see [stats::sigma()]).
#'  * `terms()` extracts model terms (see [stats::terms()]).
#' @param x,object An [`EventDate-class`] object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]). If `NULL` (the default), *rata die* are
#'  returned.
#' @param ... Currently not used.
#' @example inst/examples/ex-event.R
#' @author N. Frerebeau
#' @family event date tools
#' @docType methods
#' @name model_event
#' @rdname model_event
NULL

#' Plot Event and Accumulation Dates
#'
#' Produces an activity or a tempo plot.
#' @param x An [`EventDate-class`] object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]).
#' @param type A [`character`] string indicating the type of plot.
#'  It must be one of "`activity`" (default) or "`tempo`" (see details).
#'  Any unambiguous substring can be given.
#' @param event A [`logical`] scalar: should the distribution of the event date
#'  be displayed? Only used if type is "`activity`".
#' @param select A [`numeric`] or [`character`] vector giving the selection of
#'  the assemblage that are drawn.
#' @param n A length-one non-negative [`numeric`] vector giving the desired
#'  length of the vector of quantiles for density computation.
#' @param eps A length-one [`numeric`] value giving the cutoff below which
#'  values will be removed.
#' @param col.accumulation A color specification for the accumulation density
#'  curve.
#' @param col.event A color specification for the event density curve.
#' @inheritParams aion::plot
#' @section Event and Acccumulation Dates:
#'  `plot()` displays the probability estimate density curves of archaeological
#'  assemblage dates (*event* and *accumulation* dates; Bellanger and Husi
#'  2012). The *event* date is plotted as a line, while the *accumulation* date
#'  is shown as a grey filled area.
#'
#'  The accumulation date can be displayed as a tempo plot (Dye 2016) or an
#'  activity plot (Philippe and Vibet 2020):
#'  \describe{
#'   \item{`tempo`}{A tempo plot estimates the cumulative occurrence of
#'   archaeological events, such as the slope of the plot directly reflects the
#'   pace of change.}
#'   \item{`activity`}{An activity plot displays the first derivative of the
#'   tempo plot.}
#'  }
#' @return
#'  `plot()` is called it for its side-effects: it results in a graphic being
#'  displayed (invisibly returns `x`).
#' @references
#'  Bellanger, L. & Husi, P. (2012). Statistical Tool for Dating and
#'  Interpreting Archaeological Contexts Using Pottery. *Journal of
#'  Archaeological Science*, 39(4), 777-790. \doi{10.1016/j.jas.2011.06.031}.
#'
#'  Dye, T. S. (2016). Long-Term Rhythms in the Development of Hawaiian
#'  Social Stratification. *Journal of Archaeological Science*, 71, 1-9.
#'  \doi{10.1016/j.jas.2016.05.006}.
#'
#'  Philippe, A. & Vibet, M.-A. (2020). Analysis of Archaeological Phases Using
#'  the R Package ArchaeoPhases. *Journal of Statistical Software, Code
#'  Snippets*, 93(1), 1-25. \doi{10.18637/jss.v093.c01}.
#' @seealso [event()]
#' @example inst/examples/ex-event.R
#' @author N. Frerebeau
#' @family event date tools
#' @seealso [event()]
#' @docType methods
#' @name plot_event
#' @rdname plot_event
NULL

# Aoristic Analysis ============================================================
#' Aoristic Analysis
#'
#' Computes the aoristic sum.
#' @param x,y A [`numeric`] vector giving the lower and upper boundaries of the
#'  time intervals, respectively. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param step A length-one [`integer`] vector giving the step size, i.e. the
#'  width of each time step in the time series (defaults to \eqn{1},
#'  i.e. annual level).
#' @param start A length-one [`numeric`] vector giving the beginning of the time
#'  window.
#' @param end A length-one [`numeric`] vector giving the end of the time
#'  window.
#' @param calendar An [`aion::TimeScale-class`] object specifying the calendar
#'  of `x` and `y` (see [aion::calendar()]). Defaults to Gregorian Common Era.
#' @param weight A [`logical`] scalar: should the aoristic sum be weighted by
#'  the length of periods (default). If `FALSE` the aoristic sum is the number
#'  of elements within a time block.
#' @param groups A [`factor`] vector in the sense that `as.factor(groups)`
#'  defines the grouping. If `x` is a `list` (or a `data.frame`), `groups` can
#'  be a length-one vector giving the index of the grouping component (column)
#'  of `x`.
#' @param ... Currently not used.
#' @details
#'  Aoristic analysis is used to determine the probability of contemporaneity of
#'  archaeological sites or assemblages. The aoristic analysis distributes the
#'  probability of an event uniformly over each temporal fraction of the period
#'  considered. The aoristic sum is then the distribution of the total number of
#'  events to be assumed within this period.
#'
#'  Muller and Hinz (2018) pointed out that the overlapping of temporal
#'  intervals related to period categorization and dating accuracy is likely to
#'  bias the analysis. They proposed a weighting method to overcome this
#'  problem. This method is not implemented here (for the moment), see the
#'  [\pkg{aoristAAR} package](https://github.com/ISAAKiel/aoristAAR).
#' @return
#'  An [`AoristicSum-class`] object.
#' @references
#'  Crema, E. R. (2012). Modelling Temporal Uncertainty in Archaeological
#'  Analysis. *Journal of Archaeological Method and Theory*, 19(3): 440-61.
#'  \doi{10.1007/s10816-011-9122-3}.
#'
#'  Johnson, I. (2004). Aoristic Analysis: Seeds of a New Approach to Mapping
#'  Archaeological Distributions through Time. *In* Ausserer, K. F., Börner, W.,
#'  Goriany, M. & Karlhuber-Vöckl, L. (ed.), *Enter the Past - The E-Way into
#'  the Four Dimensions of Cultural Heritage*, Oxford: Archaeopress, p. 448-52.
#'  BAR International Series 1227.
#'  \doi{10.15496/publikation-2085}
#'
#'  Müller-Scheeßel, N. & Hinz, M. (2018). *Aoristic Research in R: Correcting
#'  Temporal Categorizations in Archaeology*. Presented at the Human History and
#'  Digital Future (CAA 2018), Tubingen, March 21.
#'  <https://www.youtube.com/watch?v=bUBukex30QI>.
#'
#'  Palmisano, A., Bevan, A. & Shennan, S. (2017). Comparing Archaeological
#'  Proxies for Long-Term Population Patterns: An Example from Central Italy.
#'  *Journal of Archaeological Science*, 87: 59-72.
#'  \doi{10.1016/j.jas.2017.10.001}.
#'
#'  Ratcliffe, J. H. (2000). Aoristic Analysis: The Spatial Interpretation of
#'  Unspecific Temporal Events. *International Journal of Geographical
#'  Information Science*, 14(7): 669-79. \doi{10.1080/136588100424963}.
#'
#'  Ratcliffe, J. H. (2002). Aoristic Signatures and the Spatio-Temporal
#'  Analysis of High Volume Crime Patterns. *Journal of Quantitative
#'  Criminology*, 18(1): 23-43. \doi{10.1023/A:1013240828824}.
#' @example inst/examples/ex-aoristic.R
#' @author N. Frerebeau
#' @family aoristic analysis
#' @docType methods
#' @aliases aoristic-method
setGeneric(
  name = "aoristic",
  def = function(x, y, ...) standardGeneric("aoristic"),
  valueClass = "AoristicSum"
)

#' Rate of Change
#'
#' Computes the rate of change from an aoristic analysis.
#' @param object An [`AoristicSum-class`] object.
#' @param n A non-negative [`integer`] giving the number of replications (see
#'  details).
#' @param ... Currently not used.
#' @return
#'  A [`RateOfChange-class`] object.
#' @references
#'  Baxter, M. J. & Cool, H. E. M. (2016). Reinventing the Wheel? Modelling
#'  Temporal Uncertainty with Applications to Brooch Distributions in Roman
#'  Britain. *Journal of Archaeological Science*, 66: 120-27.
#'  \doi{10.1016/j.jas.2015.12.007}.
#'
#'  Crema, E. R. (2012). Modelling Temporal Uncertainty in Archaeological
#'  Analysis. *Journal of Archaeological Method and Theory*, 19(3): 440-61.
#'  \doi{10.1007/s10816-011-9122-3}.
#' @example inst/examples/ex-aoristic.R
#' @author N. Frerebeau
#' @family aoristic analysis
#' @docType methods
#' @aliases roc-method
setGeneric(
  name = "roc",
  def = function(object, ...) standardGeneric("roc"),
  valueClass = "RateOfChange"
)

#' Plot Aoristic Analysis
#'
#' @param x An [`AoristicSum-class`] object.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]).
#' @param type A [`character`] string specifying whether bar or density should
#'  be plotted? It must be one of "`bar`" or "`density`". Any unambiguous
#'  substring can be given.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @inheritParams 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-aoristic.R
#' @author N. Frerebeau
#' @family aoristic analysis
#' @docType methods
#' @name plot_aoristic
#' @rdname plot_aoristic
NULL

# Chronological Apportioning ===================================================
#' Chronological Apportioning
#'
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data (absolute frequencies giving the number of
#'  individuals for each category, i.e. a contingency table). A [`data.frame`]
#'  will be coerced to a `numeric` `matrix` via [data.matrix()].
#' @param s0 A length-\eqn{m} [`numeric`] vector giving the site beginning dates
#'  expressed in CE years (BCE years must be given as negative numbers).
#' @param s1 A length-\eqn{m} [`numeric`] vector giving the site end dates
#'  expressed in CE years (BCE years must be given as negative numbers).
#' @param t0 A length-\eqn{p} [`numeric`] vector giving the type beginning dates
#'  expressed in CE years (BCE years must be given as negative numbers).
#' @param t1 A length-\eqn{p} [`numeric`] vector giving the type end dates
#'  expressed in CE years (BCE years must be given as negative numbers).
#' @param from A length-one [`numeric`] vector giving the beginning of the
#'  period of interest (in years CE).
#' @param to A length-one [`numeric`] vector giving the end of the period of
#'  interest (in years CE).
#' @param step A length-one [`integer`] vector giving the step size, i.e. the
#'  width of each time step for apportioning (in years CE; defaults to
#'  \eqn{25}).
#' @param method A [`character`] string specifying the distribution to be used
#'  (type popularity curve). It must be one of "`uniform`" (uniform
#'  distribution) or "`truncated`" (truncated standard normal distribution).
#'  Any unambiguous substring can be given.
#' @param z An [`integer`] value giving the lower and upper truncation points
#'  (defaults to \eqn{2}). Only used if `method` is "`truncated`".
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @param ... Currently not used.
#' @return
#'  A [`CountApportion-class`] object.
#' @references
#'  Roberts, J. M., Mills, B. J., Clark, J. J., Haas, W. R., Huntley, D. L. &
#'  Trowbridge, M. A. (2012). A Method for Chronological Apportioning of Ceramic
#'  Assemblages. *Journal of Archaeological Science*, 39(5): 1513-20.
#'  \doi{10.1016/j.jas.2011.12.022}.
#' @example inst/examples/ex-apportion.R
#' @author N. Frerebeau
#' @family chronological apportioning methods
#' @family chronological analysis
#' @docType methods
#' @aliases apportion-method
setGeneric(
  name = "apportion",
  def = function(object, ...) standardGeneric("apportion"),
  valueClass = "CountApportion"
)

# Frequency Increment Test =====================================================
#' Frequency Increment Test
#'
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data (absolute frequencies giving the number of
#'  individuals for each category, i.e. a contingency table). A [`data.frame`]
#'  will be coerced to a `numeric` `matrix` via [data.matrix()].
#' @param dates A length-\eqn{m} [`numeric`] vector of dates.
#' @param calendar An [`aion::TimeScale-class`] object specifying the calendar
#'  of `dates` (see [aion::calendar()]). Defaults to Gregorian Common Era.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @param roll A [`logical`] scalar: should each time series be subsetted to
#'  look for episodes of selection?
#' @param window An odd [`integer`] giving the size of the rolling
#'  window. Only used if `roll` is `TRUE`.
#' @param ... Currently not used.
#' @details
#'  The Frequency Increment Test (FIT) rejects neutrality if the distribution
#'  of normalized variant frequency increments exhibits a mean that deviates
#'  significantly from zero.
#'
#'  If `roll` is `TRUE`, each time series is subsetted according to `window` to
#'  see if episodes of selection can be identified among variables that might
#'  not show overall selection.
#' @return
#'  An [`IncrementTest-class`] object.
#' @example inst/examples/ex-fit.R
#' @references
#'  Feder, A. F., Kryazhimskiy, S. & Plotkin, J. B. (2014). Identifying
#'  Signatures of Selection in Genetic Time Series. *Genetics*, 196(2):
#'  509-522. \doi{10.1534/genetics.113.158220}.
#' @author N. Frerebeau
#' @family chronological analysis
#' @docType methods
#' @aliases fit-method
setGeneric(
  name = "fit",
  def = function(object, dates, ...) standardGeneric("fit"),
  valueClass = "IncrementTest"
)

#' Detection of Selective Processes
#'
#' Produces an abundance *vs* time diagram.
#' @param x An [`IncrementTest-class`] object to be plotted.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]).
#' @param col.neutral,col.selection,col.roll A vector of colors.
#' @inheritParams aion::plot
#' @details
#'  Results of the frequency increment test can be displayed on an abundance
#'  *vs* time diagram aid in the detection and quantification of selective
#'  processes in the archaeological record. If `roll` is `TRUE`, each time
#'  series is subsetted according to `window` to see if episodes of selection
#'  can be identified among decoration types that might not show overall
#'  selection. If so, shading highlights the data points where
#'  [fit()] identifies selection.
#' @return
#'  `plot()` is called it for its side-effects: it results in a graphic being
#'  displayed (invisibly returns `x`).
#' @note
#'  Displaying FIT results on an abundance *vs* time diagram is adapted from Ben
#'  Marwick's [original idea](https://github.com/benmarwick/signatselect/).
#' @example inst/examples/ex-fit.R
#' @author N. Frerebeau
#' @family chronological analysis
#' @docType methods
#' @name plot_fit
#' @rdname plot_fit
NULL

# Plot =========================================================================
#' Abundance vs Time Plot
#'
#' Produces an abundance *vs* time diagram.
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data (absolute frequencies giving the number of
#'  individuals for each category, i.e. a contingency table). A [`data.frame`]
#'  will be coerced to a `numeric` `matrix` via [data.matrix()].
#' @param dates A [`numeric`] vector of dates.
#' @param calendar An [`aion::TimeScale-class`] object specifying the target
#'  calendar (see [aion::calendar()]).
#' @param ... Further parameters to be passed to [aion::plot()].
#' @return
#'   `plot_time()` is called it for its side-effects: it results in a graphic
#'   being displayed (invisibly returns `object`).
#' @example inst/examples/ex-plot.R
#' @author N. Frerebeau
#' @family plotting methods
#' @docType methods
#' @aliases plot_time-method
setGeneric(
  name = "plot_time",
  def = function(object, dates, ...) standardGeneric("plot_time")
)

# Seriation Methods ============================================================
## Reciprocal ranking ----------------------------------------------------------
#' Reciprocal Ranking Seriation
#'
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data (absolute frequencies giving the number of
#'  individuals for each category, i.e. a contingency table). A [`data.frame`]
#'  will be coerced to a `numeric` `matrix` via [data.matrix()].
#' @param EPPM A [`logical`] scalar: should the seriation be computed on EPPM
#'  instead of raw data?
#' @param margin A [`numeric`] vector giving the subscripts which the
#'  rearrangement will be applied over: `1` indicates rows, `2` indicates
#'  columns, `c(1, 2)` indicates rows then columns, `c(2, 1)` indicates columns
#'  then rows.
#' @param stop An [`integer`] giving the stopping rule (i.e. maximum number of
#'  iterations) to avoid infinite loop.
#' @param ... Currently not used.
#' @details
#'  This procedure iteratively rearrange rows and/or columns according to their
#'  weighted rank in the data matrix until convergence.
#'
#'  Note that this procedure could enter into an infinite loop. If no
#'  convergence is reached before the maximum number of iterations, it stops
#'  with a warning.
#' @return
#'  A [RankPermutationOrder-class] object.
#' @references
#'  Desachy, B. (2004). Le sériographe EPPM: un outil informatisé de sériation
#'  graphique pour tableaux de comptages. *Revue archéologique de Picardie*,
#'  3(1), 39-56. \doi{10.3406/pica.2004.2396}.
#'
#'  Dunnell, R. C. (1970). Seriation Method and Its Evaluation. *American
#'  Antiquity*, 35(03), 305-319. \doi{10.2307/278341}.
#'
#'  Ihm, P. (2005). A Contribution to the History of Seriation in Archaeology.
#'  In C. Weihs & W. Gaul (Eds.), *Classification: The Ubiquitous
#'  Challenge*. Berlin Heidelberg: Springer, p. 307-316.
#'  \doi{10.1007/3-540-28084-7_34}.
#' @example inst/examples/ex-seriation.R
#' @author N. Frerebeau
#' @family seriation methods
#' @docType methods
#' @aliases seriate_rank-method
setGeneric(
  name = "seriate_rank",
  def = function(object, ...) standardGeneric("seriate_rank"),
  valueClass = "PermutationOrder"
)

## Average Ranking -------------------------------------------------------------
#' Correspondence Analysis-Based Seriation
#'
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data (absolute frequencies giving the number of
#'  individuals for each category, i.e. a contingency table). A [`data.frame`]
#'  will be coerced to a `numeric` `matrix` via [data.matrix()].
#' @param margin A [`numeric`] vector giving the subscripts which the
#'  rearrangement will be applied over: `1` indicates rows, `2` indicates
#'  columns, `c(1, 2)` indicates rows then columns, `c(2, 1)` indicates columns
#'  then rows.
#' @param axes An [`integer`] vector giving the subscripts of the CA axes to be
#'  used.
#' @param sup_row A `vector` specifying the indices of the supplementary rows
#'  (see [dimensio::ca()]).
#' @param sup_col A `vector` specifying the indices of the supplementary columns
#'  (see [dimensio::ca()]).
#' @param ... Currently not used.
#' @details
#'  Correspondence analysis (CA) is an effective method for the seriation of
#'  archaeological assemblages. The order of the rows and columns is given by
#'  the coordinates along one dimension of the CA space, assumed to account for
#'  temporal variation. The direction of temporal change within the
#'  correspondence analysis space is arbitrary: additional information is needed
#'  to determine the actual order in time.
#' @return
#'  An [`AveragePermutationOrder-class`] object.
#' @references
#'  Ihm, P. (2005). A Contribution to the History of Seriation in Archaeology.
#'  In C. Weihs & W. Gaul (Eds.), *Classification: The Ubiquitous
#'  Challenge*. Berlin Heidelberg: Springer, p. 307-316.
#'  \doi{10.1007/3-540-28084-7_34}.
#' @seealso [dimensio::ca()]
#' @example inst/examples/ex-seriation.R
#' @author N. Frerebeau
#' @family seriation methods
#' @docType methods
#' @aliases seriate_average-method
setGeneric(
  name = "seriate_average",
  def = function(object, ...) standardGeneric("seriate_average"),
  valueClass = "PermutationOrder"
)

# @rdname seriate_constrain
# @aliases seriate_constrain-method
# setGeneric(
#   name = "seriate_constrain",
#   def = function(object, constrain, ...) standardGeneric("seriate_constrain"),
#   valueClass = "PermutationOrder"
# )

# @rdname seriate_idds
# @aliases seriate_idds-method
# setGeneric(
#   name = "seriate_idds",
#   def = function(object, ...) standardGeneric("seriate_idds"),
#   valueClass = "PermutationOrder"
# )

#' Coerce an \R Object to a Seriation Order
#'
#' @param object An \R object.
#' @param margin A [`numeric`] vector giving the subscripts which the
#'  rearrangement will be applied over: `1` indicates rows, `2` indicates
#'  columns, `c(1, 2)` indicates rows then columns, `c(2, 1)` indicates columns
#'  then rows.
#' @param axes An [`integer`] vector giving the subscripts of the CA axes to be
#'  used.
#' @param ... Currently not used.
#' @return
#'  A [`PermutationOrder-class`] object.
#' @author N. Frerebeau
#' @family seriation methods
#' @docType methods
#' @aliases as_seriation-method
setGeneric(
  name = "as_seriation",
  def = function(object, ...) standardGeneric("as_seriation")
)

## Refine ----------------------------------------------------------------------
#' Refine CA-based Seriation
#'
#' @param object A [`PermutationOrder-class`] object (typically returned by
#'  [seriate_average()]) or a [`dimensio::BootstrapCA-class`] object (typically
#'  returned by [dimensio::bootstrap()]).
#' @param cutoff A function that takes a numeric vector as argument and returns
#'  a single numeric value (see below).
#' @param n A non-negative [`integer`] giving the number of bootstrap
#'  replications.
#' @param margin A length-one [`numeric`] vector giving the subscripts which the
#'  refinement will be applied over: `1` indicates rows, `2` indicates columns.
#' @param axes An [`integer`] vector giving the subscripts of the CA axes to be
#'  used.
#' @param ... Currently not used.
#' @details
#'  `refine()` allows to identify samples that are subject to sampling error or
#'  samples that have underlying structural relationships and might be
#'  influencing the ordering along the CA space.
#'
#'  This relies on a partial bootstrap approach to CA-based seriation where each
#'  sample is replicated `n` times. The maximum dimension length of the convex
#'  hull around the sample point cloud allows to remove samples for a given
#'  `cutoff` value.
#'
#'  According to Peebles and Schachner (2012), "\[this\] point removal procedure
#'  \[results in\] a reduced dataset where the position of individuals within
#'  the CA are highly stable and which produces an ordering consistent with the
#'  assumptions of frequency seriation."
#'
#'  See `vignette("seriation")`.
#' @return
#'  A [`list`] with the following elements:
#'  \describe{
#'   \item{`length`}{A [`numeric`] vector giving the convex hull maximum
#'   dimension length.}
#'   \item{`cutoff`}{A [`numeric`] value giving the cutoff value for samples
#'   selection.}
#'   \item{`exclude`}{An [`integer`] vector giving the subscript of the
#'   observations to be removed.}
#'   \item{`margin`}{A [`numeric`] value specifying the dimension along which
#'   the refinement procedure has been applied: `1` indicates rows,
#'   `2` indicates columns.}
#'  }
#' @references
#'  Peeples, M. A., & Schachner, G. (2012). Refining correspondence
#'  analysis-based ceramic seriation of regional data sets. *Journal of
#'  Archaeological Science*, 39(8), 2818-2827.
#'  \doi{10.1016/j.jas.2012.04.040}.
#' @seealso [dimensio::bootstrap()]
#' @author N. Frerebeau
#' @docType methods
#' @family seriation methods
#' @aliases refine-method
setGeneric(
  name = "refine",
  def = function(object, ...) standardGeneric("refine")
)

## Assess ----------------------------------------------------------------------
#' Statistical Significance of Seriation Solutions
#'
#' Tests the significance of seriation solutions.
#' @param object A [`PermutationOrder-class`] object giving the permutation
#'  order for rows and columns (typically returned by [seriate_average()]).
#' @param axes An [`integer`] vector giving the subscripts of the CA axes to be
#'  used.
#' @param n A non-negative [`integer`] giving the number of bootstrap
#'  replications.
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @param ... Currently not used.
#' @return
#'  A [`list`] with the following elements:
#'  \describe{
#'   \item{`random`}{A [`numeric`] vector giving the randomized total number of
#'   modes values.}
#'   \item{`observed`}{A [`numeric`] value giving the observed total number of
#'   modes.}
#'   \item{`expected`}{A [`numeric`] value giving the expected total number of
#'   modes if all types had unimodal distributions.}
#'   \item{`maximum`}{A [`numeric`] value giving the maximum total number of
#'   modes.}
#'   \item{`coef`}{A [`numeric`] value giving the seriation coefficient (a value
#'   close to 1 indicates a strong fit to the seriation model, while a value
#'   close to 0 indicates a poor fit).}
#'  }
#' @references
#'  Porčić, M. (2013). The Goodness of Fit and Statistical Significance of
#'  Seriation Solutions. *Journal of Archaeological Science*, 40(12): 4552-4559.
#'  \doi{10.1016/j.jas.2013.07.013}.
#' @example inst/examples/ex-assess.R
#' @author N. Frerebeau
#' @docType methods
#' @family seriation methods
#' @aliases assess-method
setGeneric(
  name = "assess",
  def = function(object, ...) standardGeneric("assess")
)

## Permute ---------------------------------------------------------------------
#' Rearrange a Data Matrix
#'
#' Rearranges a data matrix according to a permutation order.
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data (absolute frequencies giving the number of
#'  individuals for each category, i.e. a contingency table).
#' @param order A [`PermutationOrder-class`] object giving the permutation
#'  order for rows and columns.
#' @param ... Currently not used.
#' @return
#'  A permuted `matrix` or a permuted `data.frame` (the same as `object`).
#' @example inst/examples/ex-seriation.R
#' @author N. Frerebeau
#' @family seriation methods
#' @docType methods
#' @aliases permute-method
setGeneric(
  name = "permute",
  def = function(object, order, ...) standardGeneric("permute")
)

#' Permutation Order
#'
#' Returns the seriation order for rows and/or columns.
#' @param object A [`PermutationOrder-class`] object giving the permutation
#'  order for rows and columns.
#' @param ... Currently not used.
#' @return
#'  An [`integer`] vector.
#' @example inst/examples/ex-seriation.R
#' @author N. Frerebeau
#' @family seriation methods
#' @docType methods
#' @name order
#' @rdname order
NULL

#' @rdname order
#' @aliases order_rows-method
setGeneric(
  name = "order_rows",
  def = function(object, ...) standardGeneric("order_rows")
)

#' @rdname order
#' @aliases order_columns-method
setGeneric(
  name = "order_columns",
  def = function(object, ...) standardGeneric("order_columns")
)

# PLOT LINE
#' @include AllClasses.R AllGenerics.R
NULL

# data.frame ===================================================================
#' @export
#' @rdname plot_time
#' @aliases plot_time,data.frame,numeric-method
setMethod(
  f = "plot_time",
  signature = c(object = "data.frame", dates = "numeric"),
  definition = function(object, dates,
                        calendar = get_calendar(), ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, dates = dates, calendar = calendar, ...)
    invisible(object)
  }
)

# matrix =======================================================================
#' @export
#' @rdname plot_time
#' @aliases plot_time,matrix,numeric-method
setMethod(
  f = "plot_time",
  signature = c(object = "matrix", dates = "numeric"),
  definition = function(object, dates,
                        calendar = get_calendar(), ...) {
    ## Validation
    arkhe::assert_type(dates, "numeric")
    arkhe::assert_length(dates, nrow(object))

    ## Convert to rata die
    if (is.null(calendar)) {
      dates <- aion::as_fixed(dates)
    } else {
      dates <- aion::fixed(dates, calendar = calendar)
    }

    ## Prepare data
    ts <- aion::series(object = object, time = dates)

    ## Plot
    aion::plot(x = ts, calendar = calendar, ...)

    invisible(object)
  }
)

# COERCE
#' @include AllGenerics.R
NULL

# To data.frame ================================================================
#' @export
#' @method as.data.frame MeanDate
as.data.frame.MeanDate <- function(x, ..., calendar = get_calendar()) {
  data.frame(
    sample = rownames(x) %||% paste0("S", seq_len(nrow(x))),
    time = aion::time(x, calendar = calendar)
  )
}

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

#' @export
#' @method as.data.frame AoristicSum
as.data.frame.AoristicSum <- function(x, ..., calendar = get_calendar()) {
  block_start <- utils::head(x@breaks, -1)
  block_end <- utils::tail(x@breaks, -1)

  aorist <- x[, , 1, drop = TRUE]
  dim(aorist) <- c(nrow(x), ncol(x))
  colnames(aorist) <- colnames(x) %||% paste0("S", seq_len(ncol(x)))

  data.frame(
    start = aion::as_year(block_start, calendar = calendar),
    end = aion::as_year(block_end, calendar = calendar),
    aorist
  )
}

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

#' @method as.data.frame IncrementTest
#' @export
as.data.frame.IncrementTest <- function(x, row.names = NULL, optional = FALSE, ...) {
  data.frame(
    t = x@statistic,
    p.value = x@p_value,
    row.names = colnames(x),
    ...
  )
}

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

# REFINE MATRIX SERIATION
#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname refine
#' @aliases refine,AveragePermutationOrder-method
setMethod(
  f = "refine",
  signature = c(object = "AveragePermutationOrder"),
  definition = function(object, cutoff, margin = 1, axes = 1, n = 30, ...) {
    ## Partial bootstrap CA
    ## /!\ Be careful: AveragePermutationOrder inherits from CA
    object <- dimensio::bootstrap(object, n = n)
    methods::callGeneric(object, cutoff = cutoff, margin = margin, axes = axes, ...)
  }
)

#' @export
#' @rdname refine
#' @aliases refine,BootstrapCA-method
setMethod(
  f = "refine",
  signature = c(object = "BootstrapCA"),
  definition = function(object, cutoff, margin = 1, axes = 1, ...) {
    ## Validation
    arkhe::assert_function(cutoff)
    arkhe::assert_length(margin, 1)
    arkhe::assert_length(axes, 1)

    ## Get data
    ca_rep <- dimensio::get_replications(object, margin = margin)

    ## Compute convex hull
    hull <- apply(
      X = ca_rep,
      MARGIN = 1,
      FUN = function(x, axes) compute_chull(t(x), c(1, 2)),
      axes = axes
    )

    ## Get convex hull maximal dimension length for each sample
    len <- vapply(
      X = hull,
      FUN = function(x) max(stats::dist(x, method = "euclidean")),
      FUN.VALUE = double(1)
    )

    ## Get cutoff values
    limit <- cutoff(len)

    ## Samples to be excluded
    sup <- which(len >= limit)

    list(
      length = len,
      cutoff = limit,
      exclude = which(len >= limit),
      margin = as.integer(margin)
    )
  }
)

#' Convex Hull of CA Coordinates
#'
#' Compute convex hull area for each replicated sample
#' @param x A [`numeric`] matrix of bootstrap replicates.
#' @param axes A length-two [`numeric`] vector giving the subscripts
#'  of the CA components to use.
#' @return A [`matrix`].
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
compute_chull <- function(x, axes) {
  # Remove missing values
  clean <- stats::na.omit(x[, axes])
  # Get convex hull coordinates
  points <- grDevices::chull(clean)
  # Repeat first point
  hull <- clean[c(points, points[1]), axes]
  colnames(hull) <- c("x", "y")

  return(hull)
}

#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname permute
#' @aliases permute,data.frame,PermutationOrder-method
setMethod(
  f = "permute",
  signature = c(object = "data.frame", order = "PermutationOrder"),
  definition = function(object, order) {
    # Rearrange data.frame
    object[order@rows_order, order@columns_order]
  }
)

#' @export
#' @rdname permute
#' @aliases permute,matrix,PermutationOrder-method
setMethod(
  f = "permute",
  signature = c(object = "matrix", order = "PermutationOrder"),
  definition = function(object, order) {
    # Rearrange matrix
    object[order@rows_order, order@columns_order]
  }
)

#' @export
#' @rdname order
#' @aliases order_rows,PermutationOrder-method
setMethod(
  f = "order_rows",
  signature = c("PermutationOrder"),
  definition = function(object) object@rows_order
)

#' @export
#' @rdname order
#' @aliases order_columns,PermutationOrder-method
setMethod(
  f = "order_columns",
  signature = c("PermutationOrder"),
  definition = function(object) object@columns_order
)

1		#' @include AllClasses.R AllGenerics.R
2		NULL
3
4		#' @export
5		#' @rdname seriate_rank
6		#' @aliases seriate_rank,data.frame-method
7		setMethod(
8		f = "seriate_rank",
9		signature = c(object = "data.frame"),
10		definition = function(object, EPPM = FALSE, margin = c(1, 2), stop = 100) {
11	4x	object <- data.matrix(object)
12	4x	methods::callGeneric(object, EPPM = EPPM, margin = margin, stop = stop)
13		}
14		)
15
16		#' @export
17		#' @rdname seriate_rank
18		#' @aliases seriate_rank,matrix-method
19		setMethod(
20		f = "seriate_rank",
21		signature = c(object = "matrix"),
22		definition = function(object, EPPM = FALSE, margin = c(1, 2), stop = 100) {
23		## Validation
24	5x	margin <- as.integer(margin)
25	5x	stop <- as.integer(stop)
26
27	5x	data <- object
28	1x	if (EPPM) data <- tabula::eppm(object)
29
30		## Compute ranks
31		## margin = 1 : on rows
32		## margin = 2 : on columns
33	5x	reorder <- function(x, margin) {
34	11x	i <- seq_len(nrow(x))
35	11x	j <- seq_len(ncol(x))
36	4x	if (margin == 1) k <- colSums(t(x) * j) / rowSums(x)
37	7x	if (margin == 2) k <- colSums(x * i) / colSums(x)
38	11x	order(k)
39		}
40
41	5x	start <- 0
42	5x	index <- list(rows = seq_len(nrow(data)), columns = seq_len(ncol(data)))
43	5x	convergence <- FALSE
44	5x	while (!convergence) {
45	9x	old_index <- index
46		# Rearrange along margins
47	9x	for (k in margin) {
48	11x	index[[k]] <- index[[k]][reorder(data[index[[1]], index[[2]]], margin = k)]
49		}
50		# Loop counter
51	9x	convergence <- identical(index, old_index)
52	9x	start <- start + 1
53	9x	if (start >= stop) {
54	1x	msg <- tr_("Convergence not reached (possible infinite loop).")
55	1x	warning(msg, call. = FALSE)
56	1x	break
57		}
58		}
59
60		# New PermutationOrder object
61	5x	.RankPermutationOrder(
62	5x	rows_order = as.integer(index[[1]]),
63	5x	columns_order = as.integer(index[[2]])
64		)
65		}
66		)

1		# SUBSET
2		#' @include AllGenerics.R
3		NULL
4
5		# Extract ======================================================================
6		## [ ---------------------------------------------------------------------------
7		#' @export
8		#' @rdname subset
9		#' @aliases [,MeanDate-method
10		setMethod(
11		f = "[",
12		signature = c(x = "MeanDate"),
13		function(x, i, j, k, drop = FALSE) {
14	422x	z <- x@.Data
15	422x	time <- x@.Time
16	422x	dates <- x@dates
17
18	422x	z <- z[i, j, k, drop = drop]
19	422x	if (!missing(i)) {
20	!	if (is.character(i)) i <- match(i, dimnames(x)[1L])
21	422x	time <- time[i]
22		}
23	422x	if (!missing(j)) {
24	!	if (is.character(j)) j <- match(j, dimnames(x)[2L])
25	!	dates <- dates[j]
26		}
27
28	420x	if (isTRUE(drop)) return(z)
29	2x	methods::initialize(x, z, .Time = time, dates = dates)
30		}
31		)
32
33		#' @export
34		#' @rdname subset
35		#' @aliases [,IncrementTest-method
36		setMethod(
37		f = "[",
38		signature = c(x = "IncrementTest"),
39		function(x, i, j, k, drop = FALSE) {
40	72x	z <- x@.Data
41	72x	time <- x@.Time
42	72x	statistic <- x@statistic
43	72x	p_value <- x@p_value
44
45	72x	z <- z[i, j, k, drop = drop]
46	72x	if (!missing(i)) {
47	!	if (is.character(i)) i <- match(i, dimnames(x)[1L])
48	!	time <- time[i]
49		}
50	72x	if (!missing(j)) {
51	!	if (is.character(j)) j <- match(j, dimnames(x)[2L])
52	72x	statistic <- statistic[j]
53	72x	p_value <- p_value[j]
54		}
55
56	54x	if (isTRUE(drop)) return(z)
57	18x	methods::initialize(x, z, .Time = time, statistic = statistic, p_value = p_value)
58		}
59		)
60
61		## [[ --------------------------------------------------------------------------
62		#' Extract Parts of an Object
63		#'
64		#' @inheritParams [[
65		#' @author N. Frerebeau
66		#' @keywords internal
67		#' @noRd
68		extract_slot <- function(x, i) {
69	!	class_name <- class(x)
70	!	i <- match.arg(i, choices = methods::slotNames(class_name),
71	!	several.ok = FALSE)
72	!	data <- methods::slot(x, i)
73	!	data
74		}
75
76		#' @export
77		#' @rdname subset
78		#' @aliases [[,PermutationOrder,ANY,missing-method
79		setMethod(
80		f = "[[",
81		signature = c(x = "PermutationOrder", i = "ANY", j = "missing"),
82		definition = extract_slot
83		)

1		# MEAN CERAMIC DATE
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		# MCD ==========================================================================
6		#' @export
7		#' @rdname mcd
8		#' @aliases mcd,numeric,numeric-method
9		setMethod(
10		f = "mcd",
11		signature = c(object = "numeric", dates = "numeric"),
12		definition = function(object, dates, calendar = CE()) {
13	!	object <- matrix(object, nrow = 1)
14	!	methods::callGeneric(object, dates = dates, calendar = calendar)
15		}
16		)
17
18		#' @export
19		#' @rdname mcd
20		#' @aliases mcd,data.frame,numeric-method
21		setMethod(
22		f = "mcd",
23		signature = c(object = "data.frame", dates = "numeric"),
24		definition = function(object, dates, calendar = CE()) {
25	2x	object <- data.matrix(object)
26	2x	methods::callGeneric(object, dates = dates, calendar = calendar)
27		}
28		)
29
30		#' @export
31		#' @rdname mcd
32		#' @aliases mcd,matrix,numeric-method
33		setMethod(
34		f = "mcd",
35		signature = c(object = "matrix", dates = "numeric"),
36		definition = function(object, dates, calendar = CE()) {
37		## Validation
38	2x	arkhe::assert_length(dates, NCOL(object))
39
40	1x	if (is.null(calendar)) {
41	!	dates <- aion::as_fixed(dates) # Assume rata die
42		} else {
43	1x	dates <- aion::fixed(dates, calendar = calendar)
44		}
45
46		## Calculate MCD
47	1x	mcd_dates <- apply(X = object, MARGIN = 1, FUN = .mcd, dates = dates)
48	1x	mcd_dates <- aion::as_fixed(mcd_dates)
49
50	1x	ts <- aion::series(object, time = mcd_dates)
51	1x	rownames(ts) <- rownames(object)
52	1x	.MeanDate(ts, dates = dates)
53		}
54		)
55
56		.mcd <- function(x, dates) {
57	8400x	stats::weighted.mean(x = dates, w = x)
58		}
59
60		# Resample =====================================================================
61		#' @export
62		#' @rdname resample_mcd
63		#' @aliases bootstrap,MeanDate-method
64		setMethod(
65		f = "bootstrap",
66		signature = c(object = "MeanDate"),
67		definition = function(object, n = 1000, f = NULL,
68		calendar = get_calendar()) {
69
70	!	w <- object
71	!	m <- nrow(w)
72	!	p <- ncol(w)
73	!	seq_col <- seq_len(p)
74
75	!	dates <- aion::as_year(object@dates, calendar = calendar)
76	!	theta <- function(x, counts, dates) {
77	!	.mcd(counts[x], dates[x])
78		}
79
80	!	results <- vector(mode = "list", length = m)
81	!	for (i in seq_len(m)) {
82	!	results[[i]] <- arkhe::bootstrap(
83	!	object = seq_col,
84	!	do = theta,
85	!	n = n,
86	!	counts = w[i, , 1, drop = TRUE],
87	!	dates = dates,
88	!	f = f
89		)
90		}
91
92	!	results <- do.call(rbind, results)
93	!	rownames(results) <- rownames(w)
94	!	as.data.frame(results)
95		}
96		)
97
98		#' @export
99		#' @rdname resample_mcd
100		#' @aliases jackknife,MeanDate-method
101		setMethod(
102		f = "jackknife",
103		signature = c(object = "MeanDate"),
104		definition = function(object, f = NULL,
105		calendar = get_calendar()) {
106
107	1x	w <- object
108	1x	m <- nrow(w)
109	1x	p <- ncol(w)
110
111	1x	dates <- aion::as_year(object@dates, calendar = calendar)
112	1x	theta <- function(x, counts, dates) {
113	7980x	.mcd(counts[x], dates[x])
114		}
115
116	1x	results <- vector(mode = "list", length = m)
117	1x	for (i in seq_len(m)) {
118	420x	results[[i]] <- arkhe::jackknife(
119	420x	object = seq_len(p),
120	420x	do = theta,
121	420x	counts = w[i, , 1, drop = TRUE],
122	420x	dates = dates,
123	420x	f = f
124		)
125		}
126	1x	results <- do.call(rbind, results)
127	1x	rownames(results) <- rownames(w)
128	1x	results
129		}
130		)
131
132		#' @export
133		#' @rdname resample_mcd
134		#' @aliases simulate,MeanDate-method
135		setMethod(
136		f = "simulate",
137		signature = c(object = "MeanDate"),
138		definition = function(object, nsim = 1000) {
139
140	!	results <- apply(
141	!	X = object[, , 1, drop = TRUE],
142	!	MARGIN = 1,
143	!	FUN = resample,
144	!	do = .mcd,
145	!	n = nsim,
146	!	dates = object@dates,
147	!	f = function(x) x
148		)
149
150	!	.SimulationMeanDate(object, replications = t(results))
151		}
152		)
153
154		# Plot =========================================================================
155		#' @export
156		#' @method plot MeanDate
157		plot.MeanDate <- function(x, calendar = get_calendar(),
158		decreasing = TRUE,
159		main = NULL, sub = NULL,
160		ann = graphics::par("ann"), axes = TRUE,
161		frame.plot = axes,
162		panel.first = NULL, panel.last = NULL, ...) {
163		## Get data
164	2x	n <- NROW(x)
165	2x	n_seq <- seq_len(n)
166	2x	sites <- rownames(x) %\|\|% paste0("S", n_seq)
167
168		## Graphical parameters
169	2x	col <- list(...)$col %\|\|% graphics::par("col")
170	2x	bg <- list(...)$bg %\|\|% graphics::par("bg")
171	2x	pch <- list(...)$pch %\|\|% c(16)
172	2x	cex <- list(...)$cex %\|\|% graphics::par("cex")
173	2x	lwd <- list(...)$lwd %\|\|% graphics::par("lwd")
174	2x	cex.axis <- list(...)$cex.axis %\|\|% graphics::par("cex.axis")
175	2x	col.axis <- list(...)$col.axis %\|\|% graphics::par("col.axis")
176	2x	font.axis <- list(...)$font.axis %\|\|% graphics::par("font.axis")
177	2x	if (length(col) != n) col <- rep(col, length.out = n)
178	2x	if (length(bg) != n) bg <- rep(bg, length.out = n)
179	2x	if (length(pch) != n) pch <- rep(pch, length.out = n)
180	2x	if (length(cex) != n) cex <- rep(cex, length.out = n)
181	2x	if (length(lwd) != n) lwd <- rep(lwd, length.out = n)
182
183		## Save and restore
184	2x	mar <- graphics::par("mar")
185	2x	mar[2] <- inch2line(sites, cex = cex.axis) + 0.5
186	2x	old_par <- graphics::par(mar = mar)
187	2x	on.exit(graphics::par(old_par))
188
189		## Open new window
190	2x	grDevices::dev.hold()
191	2x	on.exit(grDevices::dev.flush(), add = TRUE)
192	2x	graphics::plot.new()
193
194		## Set plotting coordinates
195	2x	years <- time(x, calendar = calendar)
196	2x	xlim <- range(years)
197	2x	ylim <- c(1, n)
198	2x	graphics::plot.window(xlim = xlim, ylim = ylim)
199
200		## Evaluate pre-plot expressions
201	2x	panel.first
202
203		## Order data
204	2x	k <- order(years, decreasing = decreasing)
205
206		## Plot
207	2x	graphics::points(x = years[k], y = n_seq, col = col[k],
208	2x	bg = bg[k], pch = pch[k], cex = cex[k], lwd = lwd[k])
209
210		## Evaluate post-plot and pre-axis expressions
211	2x	panel.last
212
213		## Construct Axis
214	2x	if (axes) {
215	2x	aion::year_axis(side = 1, format = TRUE, calendar = calendar,
216	2x	current_calendar = calendar, cex.axis = cex.axis,
217	2x	col.axis = col.axis, font.axis = font.axis)
218	2x	graphics::axis(side = 2, at = n_seq, labels = sites[k],
219	2x	xpd = NA, cex.axis = cex.axis, las = 1,
220	2x	col.axis = col.axis, font.axis = font.axis)
221		}
222
223		## Plot frame
224	2x	if (frame.plot) {
225	2x	graphics::box()
226		}
227
228		## Add annotation
229	2x	if (ann) {
230	2x	xlab <- format(calendar)
231	2x	ylab <- NULL
232	2x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
233		}
234
235	2x	invisible(x)
236		}
237
238		#' @export
239		#' @rdname plot_mcd
240		#' @aliases plot,MeanDate,missing-method
241		setMethod("plot", c(x = "MeanDate", y = "missing"), plot.MeanDate)
242
243		#' @export
244		#' @method plot SimulationMeanDate
245		plot.SimulationMeanDate <- function(x, calendar = get_calendar(),
246		interval = "student", level = 0.80,
247		decreasing = TRUE,
248		main = NULL, sub = NULL,
249		ann = graphics::par("ann"), axes = TRUE,
250		frame.plot = axes,
251		panel.first = NULL, panel.last = NULL, ...) {
252		## Get data
253	!	n <- NROW(x)
254	!	n_seq <- seq_len(n)
255	!	sites <- rownames(x) %\|\|% paste0("S1", n_seq)
256
257		## Compute confidence interval
258	!	fun <- function(x) conf(x, level = level, type = interval)
259	!	inter <- apply(X = x@replications, MARGIN = 1, FUN = fun)
260	!	inter <- apply(X = inter, MARGIN = 1, FUN = as_year, calendar = calendar)
261
262		## Graphical parameters
263	!	col <- list(...)$col %\|\|% graphics::par("col")
264	!	bg <- list(...)$bg %\|\|% graphics::par("bg")
265	!	pch <- list(...)$pch %\|\|% c(16)
266	!	cex <- list(...)$cex %\|\|% graphics::par("cex")
267	!	lty <- list(...)$lty %\|\|% graphics::par("lty")
268	!	lwd <- list(...)$lwd %\|\|% graphics::par("lwd")
269	!	if (length(col) != n) col <- rep(col, length.out = n)
270	!	if (length(bg) != n) bg <- rep(bg, length.out = n)
271	!	if (length(pch) != n) pch <- rep(pch, length.out = n)
272	!	if (length(cex) != n) cex <- rep(cex, length.out = n)
273	!	if (length(lty) != n) lty <- rep(lty, length.out = n)
274	!	if (length(lwd) != n) lwd <- rep(lwd, length.out = n)
275	!	cex.axis <- list(...)$cex.axis %\|\|% graphics::par("cex.axis")
276	!	col.axis <- list(...)$col.axis %\|\|% graphics::par("col.axis")
277	!	font.axis <- list(...)$font.axis %\|\|% graphics::par("font.axis")
278
279		## Save and restore
280	!	mar <- graphics::par("mar")
281	!	mar[2] <- inch2line(sites, cex = cex.axis) + 0.5
282	!	old_par <- graphics::par(mar = mar)
283	!	on.exit(graphics::par(old_par))
284
285		## Open new window
286	!	grDevices::dev.hold()
287	!	on.exit(grDevices::dev.flush(), add = TRUE)
288	!	graphics::plot.new()
289
290		## Set plotting coordinates
291	!	years <- time(x, calendar = calendar)
292	!	xlim <- range(inter)
293	!	ylim <- c(1, n)
294	!	graphics::plot.window(xlim = xlim, ylim = ylim)
295
296		## Evaluate pre-plot expressions
297	!	panel.first
298
299		## Order data
300	!	k <- order(years, decreasing = decreasing)
301
302		## Plot
303	!	graphics::segments(x0 = inter[k, 2], y0 = n_seq,
304	!	x1 = inter[k, 3], y1 = n_seq,
305	!	col = col[k], lty = lty[k], lwd = lwd[k])
306	!	graphics::points(x = years[k], y = n_seq, col = col[k],
307	!	bg = bg[k], pch = pch[k], cex = cex[k], lwd = lwd[k])
308
309		## Evaluate post-plot and pre-axis expressions
310	!	panel.last
311
312		## Construct Axis
313	!	if (axes) {
314	!	aion::year_axis(side = 1, format = TRUE, calendar = calendar,
315	!	current_calendar = calendar, cex.axis = cex.axis,
316	!	col.axis = col.axis, font.axis = font.axis)
317	!	graphics::axis(side = 2, at = n_seq, labels = sites[k],
318	!	xpd = NA, cex.axis = cex.axis, las = 1,
319	!	col.axis = col.axis, font.axis = font.axis)
320		}
321
322		## Plot frame
323	!	if (frame.plot) {
324	!	graphics::box()
325		}
326
327		## Add annotation
328	!	if (ann) {
329		## Caption
330	!	cap <- sprintf(tr_("Simulated assemblages, %d replications."), NCOL(x@replications))
331	!	xlab <- format(calendar)
332	!	ylab <- NULL
333	!	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
334		}
335
336	!	invisible(x)
337		}
338
339		#' @export
340		#' @rdname plot_mcd
341		#' @aliases plot,SimulationMeanDate,missing-method
342		setMethod("plot", c(x = "SimulationMeanDate", y = "missing"),
343		plot.SimulationMeanDate)
344
345		conf <- function(x, type = c("percentiles", "student", "normal", "range"),
346		level = 0.80) {
347	!	type <- match.arg(type, several.ok = FALSE)
348	!	if (type == "range") {
349	!	conf <- range(x, na.rm = FALSE)
350	!	} else if (type == "percentiles") {
351		## Confidence interval as described in Kintigh 1989
352	!	k <- (1 - level) / 2
353	!	conf <- stats::quantile(x, probs = c(k, 1 - k), names = FALSE)
354		} else {
355		## Confidence interval
356	!	conf <- arkhe::confidence_mean(x, level = level, type = type)
357		}
358
359	!	result <- c(mean(x), conf)
360	!	names(result) <- c("mean", "lower", "upper")
361	!	result
362		}

1		# HELPERS
2
3		## https://michaelchirico.github.io/potools/articles/developers.html
4		tr_ <- function(...) {
5	3x	enc2utf8(gettext(paste0(...), domain = "R-kairos"))
6		}
7
8		make_par <- function(params, x, n = 0) {
9	!	p <- params[[x]] %\|\|% graphics::par()[[x]]
10	!	if (n > 0) p <- rep(p, length.out = n)
11	!	p
12		}
13
14		#' Plotting Dimensions of Character Strings
15		#'
16		#' Convert string length in inch to number of (margin) lines.
17		#' @param x A [`character`] vector of string whose length is to be calculated.
18		#' @param ... Further parameter to be passed to [graphics::strwidth()]`, such as
19		#' `cex`.
20		#' @return
21		#' A [`numeric`] vector (maximum string width in units of margin lines).
22		#' @note For internal use only.
23		#' @family graphic tools
24		#' @keywords internal
25		#' @noRd
26		inch2line <- function(x, ...) {
27	2x	(max(graphics::strwidth(x, units = "inch", ...)) /
28	2x	graphics::par("cin")[2] + graphics::par("mgp")[2]) * graphics::par("cex")
29		}
30
31		#' Compute x Limits
32		#'
33		#' Computes x limits for a time series according to a given calendar.
34		#' This ensures that the x axis is always in chronological order.
35		#' @param x A [`TimeSeries-class`] object.
36		#' @param calendar A [`TimeScale-class`] object.
37		#' @return A length-two [`numeric`] vector.
38		#' @keywords internal
39		#' @noRd
40		xlim <- function(x, calendar) {
41	!	x <- range(time(x, calendar = NULL))
42	!	if (is.null(calendar)) return(x)
43	!	as_year(x, calendar = calendar)
44		}
45
46		#' Indices of a rolling window
47		#'
48		#' @param x An object.
49		#' @param window An [`integer`] scalar giving the window size.
50		#' @return A [`list`] with the following components:
51		#' \describe{
52		#' \item{i}{An [`integer`] vector of indices.}
53		#' \item{w}{An [`integer`] vector of indices giving the indice of the window
54		#' mid-point.}
55		#' }
56		#' @keywords internal
57		#' @noRd
58		roll <- function(x, window = 3) {
59		## Validation
60	1x	arkhe::assert_odd(window)
61
62	1x	n <- if (is.matrix(x) \|\| is.data.frame(x)) nrow(x) else length(x)
63	1x	i <- seq_len(n) # Indices of the rows
64
65		## Matrix of rolling-window indices of length w
66	1x	w <- stats::embed(i, window)[, window:1]
67	1x	inds <- as.vector(t(w)) # Flatten indices
68
69		## Window mid-point
70	1x	m <- w[, ceiling(window / 2)]
71
72	1x	list(i = inds, w = rep(m, each = window))
73		}
74
75		#' Resample
76		#'
77		#' Simulates observations from a multinomial distribution.
78		#' @param x A [`numeric`] vector of count data (absolute frequencies).
79		#' @param do A [`function`] that takes `x` as an argument
80		#' and returns a single numeric value.
81		#' @param n A non-negative [`integer`] specifying the number of bootstrap
82		#' replications.
83		#' @param size A non-negative [`integer`] specifying the sample size.
84		#' @param f A [`function`] that takes a single numeric vector (the result of
85		#' `do`) as argument.
86		#' @param ... Extra arguments passed to `do`.
87		#' @return
88		#' If `f` is `NULL`, `resample()` returns the `n` values of `do`. Else,
89		#' returns the result of `f` applied to the `n` values of `do`.
90		#' @seealso [stats::rmultinom()]
91		#' @keywords internal
92		#' @noRd
93		resample <- function(x, do, n, size = sum(x), ..., f = NULL) {
94		## Validation
95	!	arkhe::assert_count(x)
96
97	!	prob <- x / sum(x)
98	!	replicates <- stats::rmultinom(n, size = size, prob = prob)
99	!	values <- apply(X = replicates, MARGIN = 2, FUN = do, ...)
100	!	if (is.function(f)) values <- f(values)
101	!	values
102		}
103
104		#' Quantiles of a Density Estimate
105		#'
106		#' @param x A [`numeric`] vector giving the n coordinates of the points where
107		#' the density is estimated.
108		#' @param y A `numeric` [`matrix`] of density estimates (each column is a
109		#' density estimate).
110		#' @param probs A [`numeric`] vector of probabilities with values in
111		#' \eqn{[0,1]}.
112		#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
113		#' computation proceeds?
114		#' @param ... Currently not used.
115		#' @return
116		#' A `numeric` [`matrix`] containing the quantiles.
117		#' @keywords internal
118		#' @noRd
119		quantile_density <- function(x, y, probs = seq(0, 1, 0.25), na.rm = FALSE, ...) {
120	1x	eps <- 100 * .Machine$double.eps
121	1x	if (anyNA(probs) \| any(probs < -eps \| probs > 1 + eps))
122	!	stop(sprintf("%s outside [0,1]", sQuote("probs")))
123
124	1x	q <- apply(
125	1x	X = y,
126	1x	MARGIN = 2,
127	1x	FUN = function(y, x, probs, na.rm) {
128	420x	np <- length(probs)
129	420x	qs <- rep(NA_real_, np)
130	420x	if (na.rm) {
131	!	i <- !is.na(x) & !is.na(y)
132	!	x <- x[i]
133	!	y <- y[i]
134		}
135	420x	if (np > 0) {
136	420x	nn <- length(x)
137	420x	Fx <- cumsum(y * c(0, diff(x)))
138	420x	Fx <- Fx / Fx[nn]
139	420x	for (j in seq_len(np)) {
140	1260x	ii <- min(which(Fx >= probs[j]))
141	1260x	if (!is.na(ii) && ii >= 1 && ii <= nn) qs[j] <- x[ii]
142		}
143	420x	qs
144		}
145		},
146	1x	x = x,
147	1x	probs = probs,
148	1x	na.rm = na.rm
149		)
150
151	1x	if (!is.null(dim(q))) {
152	1x	q <- t(q)
153	1x	colnames(q) <- paste0(round(probs * 100, digits = 0), "%")
154		}
155	1x	q
156		}

1		# CHRONOLOGICAL APPORTIONING
2		#' @include AllGenerics.R AllClasses.R
3		NULL
4
5		#' @export
6		#' @rdname apportion
7		#' @aliases apportion,data.frame-method
8		setMethod(
9		f = "apportion",
10		signature = c(object = "data.frame"),
11		definition = function(object, s0, s1, t0, t1, from = min(s0), to = max(s1),
12		step = 25, method = c("uniform", "truncated"), z = 2,
13		progress = getOption("kairos.progress")) {
14	!	object <- data.matrix(object)
15	!	methods::callGeneric(object, s0, s1, t0, t1, from = from, to = to,
16	!	step = step, method = method, z = z,
17	!	progress = progress)
18		}
19		)
20
21		#' @export
22		#' @rdname apportion
23		#' @aliases apportion,matrix-method
24		setMethod(
25		f = "apportion",
26		signature = c(object = "matrix"),
27		definition = function(object, s0, s1, t0, t1, from = min(s0), to = max(s1),
28		step = 25, method = c("uniform", "truncated"), z = 2,
29		progress = getOption("kairos.progress")) {
30		## Validation
31	2x	method <- match.arg(method, several.ok = FALSE)
32	2x	arkhe::assert_lower(s0, s1, strict = FALSE)
33	2x	arkhe::assert_lower(t0, t1, strict = FALSE)
34	2x	arkhe::assert_lower(from, to, strict = TRUE)
35
36		## Get data
37	2x	n_site <- nrow(object)
38	2x	n_type <- ncol(object)
39	2x	span <- to - from
40
41		## Number of periods (rounded toward the smallest integer)
42	2x	n_periode <- ceiling(span / step)
43	2x	t <- seq_len(n_periode)
44	2x	periode <- paste(from + (t - 1) * step, from + t * step, sep = "_")
45
46		## Empty array to store apportioning probabilities
47	2x	p <- array(data = NA_real_, dim = c(nrow(object), ncol(object), n_periode))
48	2x	dimnames(p) <- c(dimnames(object), list(periode))
49	2x	a <- p
50
51		## Round the site dates to be multiples of the step size
52		## (rounded toward the smallest/largest multiple)
53	2x	s0[s0 %% step != 0] <- floor(s0[s0 %% step != 0] / step) * step
54	2x	s1[s1 %% step != 0] <- ceiling(s1[s1 %% step != 0] / step) * step
55
56		## Type midpoints and life spans
57	2x	m <- (t1 + t0) / 2
58	2x	g <- t1 - t0
59
60		## Distribution function
61	2x	fun <- switch (
62	2x	method,
63	2x	uniform = dist_uniform,
64	2x	truncated = dist_truncated
65		)
66
67		## Apportion
68	2x	k_site <- seq_len(n_site)
69	2x	k_type <- seq_len(n_type)
70
71	2x	progress_bar <- interactive() && isTRUE(progress)
72	!	if (progress_bar) pbar <- utils::txtProgressBar(max = n_site, style = 3)
73
74	2x	for (i in k_site) {
75	2x	for (j in k_type) {
76		## If the type lies outside the known site occupation or study interval
77	!	if (t1[j] <= s0[i] \| t0[j] >= s1[i]) next # Do not apportion
78
79		## Earliest and latest overlaps between the ware and site
80	18x	vj0 <- max(s0[i], t0[j])
81	18x	vj1 <- min(s1[i], t1[j])
82
83		## Get the apportioning probability for any period
84	18x	qjt0 <- vapply(X = s0[i] + (t - 1) * step, FUN = max,
85	18x	FUN.VALUE = numeric(1), t0[j])
86	18x	qjt1 <- vapply(X = s0[i] + t * step, FUN = min,
87	18x	FUN.VALUE = numeric(1), t1[j])
88
89	18x	p[i, j, ] <- fun(vj0, vj1, qjt0, qjt1, g[j], m[j], z)
90		}
91
92	!	if (progress_bar) utils::setTxtProgressBar(pbar, i)
93		}
94
95	!	if (progress_bar) close(pbar)
96
97		## Remove negative values
98	2x	p[p < 0] <- 0
99
100		## Apportion
101	2x	a[] <- apply(X = p, MARGIN = 3, FUN = function(x, counts) x * counts,
102	2x	counts = object)
103
104	2x	.CountApportion(
105	2x	a,
106	2x	p = p,
107	2x	method = method,
108	2x	from = from,
109	2x	to = to,
110	2x	step = step
111		)
112		}
113		)
114
115		dist_uniform <- function(v_j0, v_j1, q_jt0, q_jt1, ...) {
116	9x	(q_jt1 - q_jt0) / (v_j1 - v_j0)
117		}
118		dist_truncated <- function(v_j0, v_j1, q_jt0, q_jt1, g, m, z) {
119	9x	z_j0 <- (v_j0 - m) / (g / (2 * z))
120	9x	z_j1 <- (v_j1 - m) / (g / (2 * z))
121
122	9x	z_jt0 <- (q_jt0 - m) / (g / (2 * z))
123	9x	z_jt1 <- (q_jt1 - m) / (g / (2 * z))
124
125	9x	(phi(z_jt1, z) - phi(z_jt0, z)) / (phi(z_j1, z) - phi(z_j0, z))
126		}
127		phi <- function(x, z) {
128	36x	extraDistr::ptnorm(x, mean = 0, sd = 1, a = -z, b = z)
129		}

1		# AORISTIC ANALYSIS
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		# Aoristic sum =================================================================
6		#' @export
7		#' @rdname aoristic
8		#' @aliases aoristic,numeric,numeric-method
9		setMethod(
10		f = "aoristic",
11		signature = c(x = "numeric", y = "numeric"),
12		definition = function(x, y, step = 1, start = min(x), end = max(y),
13		calendar = CE(), weight = TRUE, groups = NULL) {
14		## Validation
15	6x	n <- length(x)
16	6x	arkhe::assert_length(y, n)
17	4x	if (is.null(groups)) groups <- rep(NA, n)
18	6x	groups <- as.character(groups)
19	6x	arkhe::assert_length(groups, n)
20
21		## Grouping
22	5x	lvl <- unique(groups)
23	5x	grp <- factor(groups, levels = lvl, exclude = NULL)
24	5x	spl <- split(seq_along(x), f = grp)
25	5x	m <- length(spl)
26
27		## Convert to rata die
28		## NB: argument are not evaluated until they are used.
29		## ('start' and 'end' must go first)
30	5x	breaks <- seq(from = start, to = end, by = step)
31	5x	breaks <- aion::fixed(breaks, calendar = calendar)
32	5x	x <- aion::fixed(x, calendar = calendar)
33	5x	y <- aion::fixed(y, calendar = calendar)
34
35		## Set up breaks
36	5x	block_start <- utils::head(breaks, -1)
37	5x	block_end <- utils::tail(breaks, -1)
38	5x	block_mid <- (block_start + block_end) / 2
39
40		## Number of time blocks
41	5x	n_blocks <- length(breaks) - 1
42	5x	ao_probs <- array(data = 0, dim = c(n, n_blocks, m))
43
44		## Aoristic probability
45	5x	span <- abs(y - x)
46	5x	for (k in seq_len(m)) {
47	6x	s <- spl[[k]]
48	6x	for (j in seq_len(n_blocks)) {
49	24x	a <- block_start[[j]]
50	24x	b <- block_end[[j]]
51
52	24x	max_start <- pmax(x[s], a)
53	24x	min_end <- pmin(y[s], b)
54	24x	overlap <- (min_end - max_start) / span[s]
55	24x	overlap[overlap <= 0] <- 0
56	4x	if (!weight) overlap[overlap > 0] <- 1
57	24x	ao_probs[s, j, k] <- overlap
58		}
59		}
60
61		## Aoristic sum
62	5x	ao_sum <- apply(X = ao_probs, MARGIN = c(2, 3), FUN = sum)
63
64	5x	if (!anyNA(groups)) {
65	1x	dimnames(ao_probs)[[3]] <- colnames(ao_sum) <- lvl
66		}
67
68	5x	ts <- aion::series(
69	5x	ao_sum,
70	5x	time = block_mid
71		)
72	5x	.AoristicSum(
73	5x	ts,
74	5x	breaks = breaks,
75	5x	span = span,
76	5x	groups = groups,
77	5x	p = ao_probs
78		)
79		}
80		)
81
82		#' @export
83		#' @rdname aoristic
84		#' @aliases aoristic,ANY,missing-method
85		setMethod(
86		f = "aoristic",
87		signature = c(x = "ANY", y = "missing"),
88		definition = function(x, step = 1, start = NULL, end = NULL,
89		calendar = CE(), weight = TRUE, groups = NULL) {
90	7x	xy <- grDevices::xy.coords(x)
91	7x	g <- groups
92	7x	if (is.list(x) && !is.null(g) && length(g) == 1) {
93	1x	g <- x[[g]]
94	1x	if (is.null(g)) {
95	1x	txt <- tr_("%s is a list, but does not have component %s.")
96	1x	msg <- sprintf(txt, sQuote("x"), sQuote(groups))
97	1x	stop(msg, call. = FALSE)
98		}
99		}
100
101		## Start/end
102	6x	if (is.null(start)) start <- min(xy$x, na.rm = TRUE)
103	6x	if (is.null(end)) end <- max(xy$y, na.rm = TRUE)
104
105	6x	methods::callGeneric(x = xy$x, y = xy$y, step = step, start = start,
106	6x	end = end, calendar = calendar, weight = weight,
107	6x	groups = g)
108		}
109		)
110
111		# Rate of change ===============================================================
112		#' @export
113		#' @rdname roc
114		#' @aliases roc,AoristicSum-method
115		setMethod(
116		f = "roc",
117		signature = c(object = "AoristicSum"),
118		definition = function(object, n = 100) {
119		## Get time span of the blocks
120	1x	mid <- aion::time(object, calendar = NULL)
121	1x	step <- mean(diff(mid))
122	1x	end <- utils::tail(mid, -1) + step / 2
123
124		## Get aoristic weights
125	1x	w <- object@p
126
127	1x	n <- as.integer(n)
128	1x	m <- dim(w)[[1]]
129	1x	p <- dim(w)[[2]]
130	1x	q <- dim(w)[[3]]
131
132		## Monte Carlo simulation
133	1x	roc <- array(data = 0, dim = c(p-1, n, q))
134	1x	for (j in seq_len(q)) {
135	2x	sim <- array(data = 0, dim = c(m, p, n))
136	2x	for (i in seq_len(m)) {
137	20x	probs <- w[i, , j]
138	10x	if (sum(probs) == 0) next
139	10x	counts <- replicate(n, {
140	1000x	tmp <- numeric(p)
141	1000x	k <- sample(seq_len(p), size = 1, prob = probs)
142	1000x	tmp[k] <- 1
143	1000x	tmp
144		})
145	10x	sim[i, , ] <- counts
146		}
147	2x	total <- apply(X = sim, MARGIN = c(2, 3), FUN = sum)
148	2x	roc[, , j] <- apply(X = total, MARGIN = 2, FUN = diff, lag = 1)
149		}
150
151	1x	dimnames(roc)[[3]] <- dimnames(w)[[3]]
152
153	1x	ts <- aion::series(
154	1x	object = roc / step,
155	1x	time = end
156		)
157	1x	.RateOfChange(
158	1x	ts,
159	1x	replicates = n,
160	1x	groups = unique(object@groups)
161		)
162		}
163		)
164
165		# Plot =========================================================================
166		## AoristicSum -----------------------------------------------------------------
167		#' @export
168		#' @method plot AoristicSum
169		plot.AoristicSum <- function(x, calendar = get_calendar(),
170		type = c("bar"),
171		flip = FALSE, ncol = NULL,
172		main = NULL, sub = NULL,
173		ann = graphics::par("ann"), axes = TRUE,
174		frame.plot = axes,
175		panel.first = NULL, panel.last = NULL, ...) {
176		## Validation
177	3x	type <- match.arg(type, several.ok = FALSE)
178
179		## Plot
180	3x	panel_bar <- function(x, y, ...) {
181	4x	step <- mean(diff(x))
182	4x	start <- x - step / 2
183	4x	end <- x + step / 2
184	4x	graphics::rect(xleft = start, xright = end, ybottom = 0, ytop = y, ...)
185		}
186	3x	panel_density <- function(x, y, ...) {
187	!	step <- mean(diff(x))
188	!	s <- sample(x = x, size = 10^5, replace = TRUE, prob = y)
189	!	d <- stats::density(x = s, bw = step)
190	!	x <- d$x
191	!	y <- d$y
192	!	graphics::polygon(x = c(x, rev(x)), y = c(y, numeric(length(y))),
193	!	border = NA, ...)
194	!	graphics::lines(x = x, y = y)
195		}
196
197	3x	panel_fun <- switch (
198	3x	type,
199	3x	bar = panel_bar,
200	3x	density = panel_density
201		)
202
203		## Method for TimeSeries
204	3x	methods::callNextMethod(x, facet = "multiple", calendar = calendar,
205	3x	panel = panel_fun, flip = flip, ncol = ncol,
206	3x	xlab = NULL, ylab = NULL,
207	3x	main = main, sub = sub,
208	3x	ann = ann, axes = axes,
209	3x	frame.plot = frame.plot,
210	3x	panel.first = panel.first,
211	3x	panel.last = panel.last, ...)
212
213
214	3x	invisible(x)
215		}
216
217		#' @export
218		#' @rdname plot_aoristic
219		#' @aliases plot,AoristicSum,missing-method
220		setMethod("plot", c(x = "AoristicSum", y = "missing"), plot.AoristicSum)
221
222		#' @export
223		#' @method image AoristicSum
224		image.AoristicSum <- function(x, calendar = get_calendar(), ...) {
225	!	methods::callNextMethod()
226	!	invisible(x)
227		}
228
229		#' @export
230		#' @rdname plot_aoristic
231		#' @aliases image,AoristicSum-method
232		setMethod("image", c(x = "AoristicSum"), image.AoristicSum)
233
234		## RateOfChange ----------------------------------------------------------------
235		#' @export
236		#' @method plot RateOfChange
237		plot.RateOfChange <- function(x, calendar = get_calendar(),
238		level = 0.95, flip = FALSE, ncol = NULL,
239		main = NULL, sub = NULL,
240		ann = graphics::par("ann"), axes = TRUE,
241		frame.plot = axes,
242		panel.first = NULL, panel.last = NULL, ...) {
243		## Validation
244	1x	n <- dim(x)[[3L]]
245	1x	n_seq <- seq_len(n)
246	1x	if (is.null(ncol)) ncol <- if (n > 4) 2 else 1
247	1x	n_row <- ceiling(n / ncol)
248	1x	ylabs <- dimnames(x)[[3L]] %\|\|% c(tr_("Rate of Change"))
249
250		## Compute mean and confidence interval
251	1x	ci <- apply(
252	1x	X = x,
253	1x	MARGIN = c(1, 3),
254	1x	FUN = function(x, level) {
255	6x	ci <- arkhe::confidence_mean(x, level = level)
256	6x	c(mean = mean(x), ci)
257		},
258	1x	level = level
259		)
260
261		## Graphical parameters
262		## Save and restore
263	1x	old_par <- graphics::par(
264	1x	mar = c(0, 5.1, 0, if (flip) 5.1 else 2.1),
265	1x	oma = c(6, 0, 5, 0),
266	1x	mfcol = c(n_row, ncol)
267		)
268	1x	on.exit(graphics::par(old_par))
269
270	1x	cex.lab <- list(...)$cex.lab %\|\|% graphics::par("cex.lab")
271	1x	col.lab <- list(...)$col.lab %\|\|% graphics::par("col.lab")
272	1x	font.lab <- list(...)$font.lab %\|\|% graphics::par("font.lab")
273	1x	cex.axis <- list(...)$cex.axis %\|\|% graphics::par("cex.axis")
274	1x	col.axis <- list(...)$col.axis %\|\|% graphics::par("col.axis")
275	1x	font.axis <- list(...)$font.axis %\|\|% graphics::par("font.axis")
276	1x	cex.main <- list(...)$cex.main %\|\|% graphics::par("cex.main")
277	1x	font.main <- list(...)$font.main %\|\|% graphics::par("font.main")
278	1x	col.main <- list(...)$col.main %\|\|% graphics::par("col.main")
279	1x	col <- list(...)$col %\|\|% c("grey")
280	1x	if (length(col) != n) col <- rep(col, n)
281
282	1x	years <- aion::time(x, calendar = calendar)
283	1x	xlim <- range(years)
284	1x	for (i in n_seq) {
285	2x	xi <- ci[, , i, drop = TRUE]
286
287		## Open new window
288	2x	grDevices::dev.hold()
289	2x	on.exit(grDevices::dev.flush(), add = TRUE)
290	2x	graphics::plot.new()
291
292		## Set plotting coordinates
293	2x	ylim <- range(xi, na.rm = TRUE)
294	2x	graphics::plot.window(xlim = xlim, ylim = ylim)
295
296		## Evaluate pre-plot expressions
297	2x	panel.first
298
299		## Plot
300	2x	graphics::abline(h = 0, lty = "dotted", lwd = 1)
301	2x	graphics::polygon(x = c(years, rev(years)),
302	2x	y = c(xi[2, ], rev(xi[3, ])),
303	2x	col = col[i], border = NA)
304	2x	graphics::lines(x = years, y = xi[1, ])
305	2x	graphics::points(x = years, y = xi[1, ], pch = 16)
306
307		## Evaluate post-plot and pre-axis expressions
308	2x	panel.last
309
310		## Construct Axis
311	2x	do_x <- i %% n_row == 0 \|\| i == n
312	2x	y_side <- if (i %% 2 \|\| !flip) 2 else 4
313	2x	if (axes) {
314	2x	if (do_x) {
315	1x	year_axis(side = 1, format = TRUE, calendar = calendar,
316	1x	current_calendar = calendar, cex.axis = cex.axis,
317	1x	col.axis = col.axis, font.axis = font.axis)
318		}
319	2x	graphics::axis(side = y_side, xpd = NA, cex.axis = cex.axis,
320	2x	col.axis = col.axis, font.axis = font.axis, las = 1)
321		}
322
323		## Plot frame
324	2x	if (frame.plot) {
325	2x	graphics::box()
326		}
327
328		## Add annotation
329	2x	if (ann) {
330	2x	if (do_x) {
331	1x	xlab <- format(calendar)
332	1x	graphics::mtext(xlab, side = 1, line = 3, cex = cex.lab, col = col.lab,
333	1x	font = font.lab)
334		}
335	2x	graphics::mtext(ylabs[[i]], side = y_side, line = 4, cex = cex.lab,
336	2x	col = col.lab, font = font.lab)
337		}
338		}
339
340		## Add annotation
341	1x	if (ann) {
342	1x	graphics::par(mfcol = c(1, 1))
343	1x	graphics::mtext(main, side = 3, line = 3, cex = cex.main, font = font.main,
344	1x	col = col.main)
345		}
346
347	1x	invisible(x)
348		}
349
350		#' @export
351		#' @rdname plot_aoristic
352		#' @aliases plot,RateOfChange,missing-method
353		setMethod("plot", c(x = "RateOfChange", y = "missing"), plot.RateOfChange)

1		# FREQUENCY INCREMENT TEST
2		#' @include AllGenerics.R AllClasses.R
3		NULL
4
5		# FIT ==========================================================================
6		#' @export
7		#' @rdname fit
8		#' @aliases fit,data.frame,numeric-method
9		setMethod(
10		f = "fit",
11		signature = c(object = "data.frame", dates = "numeric"),
12		definition = function(object, dates, calendar = CE(),
13		level = 0.95, roll = FALSE, window = 3) {
14	2x	object <- data.matrix(object)
15	2x	methods::callGeneric(object, dates = dates, calendar = calendar,
16	2x	level = level, roll = roll, window = window)
17		}
18		)
19
20		#' @export
21		#' @rdname fit
22		#' @aliases fit,matrix,numeric-method
23		setMethod(
24		f = "fit",
25		signature = c(object = "matrix", dates = "numeric"),
26		definition = function(object, dates, calendar = CE(),
27		level = 0.95, roll = FALSE, window = 3) {
28		## Validation
29	2x	alpha <- 1 - level
30	2x	half_window <- (window - 1) / 2
31	2x	arkhe::assert_length(dates, nrow(object))
32
33		## Convert to rata die
34	2x	if (is.null(calendar)) {
35	2x	dates <- aion::as_fixed(dates)
36		} else {
37	!	dates <- aion::fixed(dates, calendar = calendar)
38		}
39
40		## Compute test
41	2x	results <- test_FIT(object, dates, roll = FALSE)
42	2x	results <- results[[1L]]
43
44		## Check results
45	2x	failed <- is.na(results$p.value)
46	2x	if (any(failed)) {
47	!	txt <- ngettext(sum(failed), "%d element was skipped: %s.",
48	!	"%d elements were skipped: %s.")
49	!	msg <- sprintf(txt, sum(failed),
50	!	paste0(rownames(results)[failed], collapse = ", "))
51	!	warning(msg, call. = FALSE)
52		}
53
54		## Rolling window
55	2x	fit_roll <- matrix(FALSE, nrow = nrow(object), ncol = ncol(object))
56	2x	if (roll) {
57	1x	results_roll <- test_FIT(object, dates, roll = roll, window = window)
58	1x	results_roll <- lapply(
59	1x	X = results_roll,
60	4x	FUN = function(x, alpha) { x$p.value <= alpha },
61	1x	alpha = alpha
62		)
63	1x	fit_row <- seq(from = half_window + 1, to = nrow(fit_roll) - half_window)
64	1x	fit_roll[fit_row, ] <- do.call(rbind, results_roll)
65
66	1x	fit_roll <- apply(
67	1x	X = fit_roll,
68	1x	MARGIN = 2,
69	1x	FUN = function(x, half_window) {
70	9x	vapply(
71	9x	X = seq_along(x),
72	9x	FUN = function(x, y, k) {
73	72x	max <- length(y)
74	72x	lower <- x - k
75	72x	lower[lower < 1] <- 1
76	72x	upper <- x + k
77	72x	upper[upper > max] <- max
78	72x	any(y[lower:upper], na.rm = TRUE)
79		},
80	9x	FUN.VALUE = logical(1),
81	9x	y = x,
82	9x	k = half_window
83		)
84		},
85	1x	half_window = half_window
86		)
87		}
88
89		## Build array
90	2x	select <- results$p.value <= alpha
91	2x	x <- y <- z <- object
92	2x	x[, !select] <- NA
93	2x	y[, select] <- NA
94	2x	z[!fit_roll] <- NA
95
96	2x	res <- array(
97	2x	data = c(z, x, y),
98	2x	dim = c(dim(object), 3),
99	2x	dimnames = c(dimnames(object), list(c("roll", "selection", "neutral")))
100		)
101
102	2x	ts <- aion::series(object = res, time = dates)
103	2x	.IncrementTest(
104	2x	ts,
105	2x	statistic = results$t,
106	2x	parameter = 1L,
107	2x	p_value = results$p.value
108		)
109		}
110		)
111
112		#' @param x A \eqn{m x p} [`numeric`] matrix.
113		#' @param time A length-\eqn{m} [`numeric`] vector.
114		#' @param roll A [`logical`] scalar.
115		#' @param window An [`integer`].
116		#' @return A list of \eqn{p x 3} [`data.frame`].
117		#' @author N. Frerebeau
118		#' @keywords internal
119		#' @noRd
120		test_FIT <- function(x, time, roll = FALSE, window = 3, ...) {
121		## Prepare data
122	3x	time <- as.numeric(time)
123
124		## Compute frequency as ratio of count of type of interest to all other types
125	3x	count_others <- lapply(
126	3x	X = seq_len(ncol(x)),
127	27x	FUN = function(i, data) { rowSums(data[, -i]) },
128	3x	data = x
129		)
130	3x	freq <- x / do.call(cbind, count_others)
131
132		## Compute test
133	3x	if (!roll) {
134	2x	roll_list <- list(seq_len(nrow(freq)))
135		} else {
136	1x	roll_index <- roll(freq, window = window)
137	1x	roll_list <- split(x = roll_index[["i"]], f = roll_index[["w"]])
138		}
139
140	3x	k <- 1
141	3x	results <- vector(mode = "list", length = length(roll_list))
142	3x	for (i in roll_list) {
143	6x	roll_freq <- freq[i, ]
144	6x	roll_time <- time[i]
145	6x	roll_date <- roll_time[ceiling(window / 2)]
146
147	6x	fit <- apply(
148	6x	X = roll_freq,
149	6x	MARGIN = 2,
150	6x	FUN = function(v, t, ...) {
151	54x	tryCatch(
152	54x	error = function(cnd) c(t = NA, p.value = NA),
153	54x	{ FIT(v, t, ...) }
154		)
155		},
156	6x	t = roll_time,
157		...
158		)
159
160	6x	results[[k]] <- data.frame(
161	6x	column = colnames(fit),
162	6x	time = roll_date,
163	6x	t(fit),
164	6x	row.names = NULL
165		)
166	6x	k <- k + 1
167		}
168
169	3x	return(results)
170		}
171
172		#' Frequency Increment Test (FIT)
173		#'
174		#' @param v A [`numeric`] vector of frequencies.
175		#' @param t A [`numeric`] vector of time coordinates.
176		#' @param ... Extra parameter passed to [stats::t.test()].
177		#' @return A [`numeric`] vector containing the following components:
178		#' \describe{
179		#' \item{`t`}{the value of the test statistic.}
180		#' \item{`p.value`}{the p-value for the test.}
181		#' }
182		#' @author N. Frerebeau
183		#' @references
184		#' Feder, A. F., Kryazhimskiy, S. & Plotkin, J. B. (2014). Identifying
185		#' Signatures of Selection in Genetic Time Series. Genetics, 196(2),
186		#' 509-522. \doi{10.1534/genetics.113.158220}.
187		#' @keywords internal
188		#' @noRd
189		FIT <- function(v, t, ...) {
190		## Validation
191	54x	arkhe::assert_type(v, "numeric")
192	54x	arkhe::assert_type(t, "numeric")
193	54x	arkhe::assert_length(t, length(v))
194
195		## Remove zeros
196	54x	index <- v > 0
197	54x	v <- v[index]
198	54x	t <- t[index]
199	54x	if (length(t) < 3) {
200	!	stop(tr_("A minimum of three time points is needed."), call. = FALSE)
201		}
202
203		## Rescaled allele-frequency increments
204	54x	Y <- diff(v, lag = 1) /
205	54x	sqrt(2 * utils::head(v, -1) * (1 - utils::head(v, -1)) * diff(t, lag = 1))
206
207		## Statistics
208	54x	t_test <- stats::t.test(Y, ...)
209	54x	c(t = as.numeric(t_test$statistic), p.value = t_test$p.value)
210		}
211
212		# Plot =========================================================================
213		#' @export
214		#' @method plot IncrementTest
215		plot.IncrementTest <- function(x, calendar = get_calendar(),
216		col.neutral = "#004488", col.selection = "#BB5566",
217		col.roll = "grey", flip = FALSE, ncol = NULL,
218		xlab = NULL, ylab = NULL,
219		main = NULL, sub = NULL,
220		ann = graphics::par("ann"), axes = TRUE,
221		frame.plot = axes, ...) {
222		## Panel
223	2x	panel_fit <- function(x, y, ...) {
224	54x	graphics::lines(x, y, ...)
225	54x	graphics::points(x, y, ...)
226		}
227
228		## Method for TimeSeries
229	2x	methods::callNextMethod(x = x, calendar = calendar, panel = panel_fit,
230	2x	flip = flip, ncol = ncol,
231	2x	xlab = xlab, ylab = ylab,
232	2x	main = main, sub = sub,
233	2x	ann = ann, axes = axes,
234	2x	frame.plot = frame.plot,
235	2x	col = c(col.roll, col.selection, col.neutral),
236	2x	lwd = c(10, 1, 1), pch = 16, cex = 1, ...)
237		}
238
239		#' @export
240		#' @rdname plot_fit
241		#' @aliases plot,IncrementTest,missing-method
242		setMethod("plot", c(x = "IncrementTest", y = "missing"), plot.IncrementTest)

1		# EVENT MODEL
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @method summary EventDate
7		summary.EventDate <- function(object, ...) {
8	1x	summary(object@model, ...)
9		}
10
11		#' @export
12		#' @rdname model_event
13		#' @aliases summary,EventDate,missing-method
14		setMethod("summary", c(object = "EventDate"), summary.EventDate)
15
16		#' @method coef EventDate
17		#' @export
18		coef.EventDate <- function(object, calendar = NULL, ...) {
19	18x	z <- stats::coef(object@model, ...)
20	18x	if (is.null(calendar)) return(z)
21	!	unclass(z) / aion::calendar_year(calendar) # Approximate
22		}
23
24		#' @export
25		#' @rdname model_event
26		#' @aliases coef,EventDate-method
27		setMethod("coef", "EventDate", coef.EventDate)
28
29		#' @method fitted EventDate
30		#' @export
31		fitted.EventDate <- function(object, calendar = NULL, ...) {
32	!	z <- stats::fitted(object@model, ...)
33	!	if (is.null(calendar)) return(z)
34	!	aion::as_year(z, calendar = calendar)
35		}
36
37		#' @export
38		#' @rdname model_event
39		#' @aliases fitted,EventDate-method
40		setMethod("fitted", "EventDate", fitted.EventDate)
41
42		#' @method residuals EventDate
43		#' @export
44		residuals.EventDate <- function(object, calendar = NULL, ...) {
45	!	z <- stats::residuals(object@model, ...)
46	!	if (is.null(calendar)) return(z)
47	!	unclass(z) / aion::calendar_year(calendar) # Approximate
48		}
49
50		#' @export
51		#' @rdname model_event
52		#' @aliases residuals,EventDate-method
53		setMethod("residuals", "EventDate", residuals.EventDate)
54
55		#' @method sigma EventDate
56		#' @export
57		sigma.EventDate <- function(object, calendar = NULL, ...) {
58	!	z <- stats::sigma(object@model, ...)
59	!	if (is.null(calendar)) return(z)
60	!	unclass(z) / aion::calendar_year(calendar) # Approximate
61		}
62
63		#' @export
64		#' @rdname model_event
65		#' @aliases sigma,EventDate-method
66		setMethod("sigma", "EventDate", sigma.EventDate)
67
68		#' @method terms EventDate
69		#' @export
70	!	terms.EventDate <- function(x, ...) stats::terms(x@model, ...)
71
72		#' @export
73		#' @rdname model_event
74		#' @aliases terms,EventDate-method
75		setMethod("terms", "EventDate", terms.EventDate)

1		# SHOW METHODS
2		#' @include AllClasses.R
3		NULL
4
5		# EventDate ====================================================================
6		setMethod(
7		f = "show",
8		signature = "EventDate",
9		definition = function(object) {
10	!	show_ca <- utils::capture.output(methods::callNextMethod(object))
11	!	summary_lm <- summary(object)
12	!	sig <- sigma(object, calendar = get_calendar())
13
14	!	cat(
15	!	show_ca,
16		"",
17	!	paste0(tr_("Multiple Linear Regression"), ":"),
18	!	sprintf(tr_("* R-squared: %s"), round(summary_lm$r.squared, 3)),
19	!	sprintf(tr_("* Adjusted R-squared: %s"), round(summary_lm$adj.r.squared, 3)),
20	!	sprintf(tr_("* Residual standard error: %s"), round(sig, 2)),
21	!	sep = "\n"
22		)
23		}
24		)
25
26		# PermutationOrder =============================================================
27		setMethod(
28		f = "show",
29		signature = "PermutationOrder",
30		definition = function(object) {
31	!	k <- 50
32	!	rows <- strtrim(paste0(object@rows_order, collapse = " "), k)
33	!	columns <- strtrim(paste0(object@columns_order, collapse = " "), k)
34	!	cat(
35	!	paste0(tr_("Permutation order for matrix seriation"), ":"),
36	!	sprintf(tr_("* Row order: %s"), paste0(rows, "...")),
37	!	sprintf(tr_("* Column order: %s"), paste0(columns, "...")),
38	!	sep = "\n"
39		)
40		}
41		)
42
43		# RefineCA =====================================================================
44		setMethod(
45		f = "show",
46		signature = "RefinePermutationOrder",
47		definition = function(object) {
48	!	value <- switch (object@margin, `1` = "Rows", `2` = "Columns")
49	!	keep <- length(object@keep)
50	!	total <- length(object@length)
51	!	pc <- round(keep * 100 / total)
52	!	methods::callNextMethod(object)
53	!	cat(
54	!	paste0("Partial bootstrap refinement", ":"),
55	!	sprintf("* Cutoff value: %s", round(object@cutoff, digits = 2)),
56	!	sprintf("* %s to keep: %d of %d (%g%%)", value, keep, total, pc),
57	!	sep = "\n")
58		}
59		)

1		# DATE MODEL
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname event
7		#' @aliases event,data.frame,numeric-method
8		setMethod(
9		f = "event",
10		signature = c(object = "data.frame", dates = "numeric"),
11		definition = function(object, dates, rank = NULL, sup_row = NULL,
12		calendar = CE(), ...) {
13	20x	object <- data.matrix(object)
14	20x	methods::callGeneric(object, dates, rank = rank, sup_row = sup_row,
15	20x	calendar = calendar, ...)
16		}
17		)
18
19		#' @export
20		#' @rdname event
21		#' @aliases event,matrix,numeric-method
22		setMethod(
23		f = "event",
24		signature = c(object = "matrix", dates = "numeric"),
25		definition = function(object, dates, calendar = CE(),
26		rank = NULL, sup_row = NULL, total = 5,
27		verbose = getOption("kairos.verbose"), ...) {
28		## Sample
29	20x	n <- nrow(object)
30
31		## Check dates
32	20x	if (arkhe::has_names(dates)) {
33	1x	i <- match(names(dates), rownames(object))
34	1x	old_dates <- dates
35	1x	dates <- rep(NA_real_, n)
36	1x	dates[i] <- old_dates
37		}
38	20x	arkhe::assert_length(dates, nrow(object))
39
40		## Supplementary rows
41	20x	sup <- logical(n)
42	20x	sup[sup_row] <- TRUE
43
44	20x	data_ref <- object[!sup, , drop = FALSE]
45	20x	data_sup <- object[sup, , drop = FALSE]
46	20x	dates_ref <- dates[!sup]
47	20x	dates_sup <- dates[sup]
48
49		## Cleansing
50	20x	data_clean <- clean_total(data_ref, data_sup, dates_ref, dates_sup,
51	20x	total = total, verbose = verbose)
52
53	20x	data <- rbind(data_clean$data_ref, data_clean$data_sup)
54	20x	dates <- c(data_clean$dates_ref, data_clean$dates_sup)
55	20x	sup <- seq_along(data_clean$dates_sup) + length(data_clean$dates_ref)
56
57		## Correspondance analysis
58	20x	if (is.null(rank)) {
59	!	tmp <- dimensio::ca(data, rank = NULL, sup_row = sup, ...)
60	!	eig <- dimensio::get_eigenvalues(tmp)
61	!	rank <- which.max(eig$cumulative >= 60)
62		}
63	20x	rank <- min(rank, dim(data) - 1)
64	20x	results_CA <- dimensio::ca(data, rank = rank, sup_row = sup, ...)
65
66		## Get row coordinates
67	20x	row_coord <- dimensio::get_coordinates(results_CA, margin = 1)
68
69		## Rata die
70	20x	if (is.null(calendar)) {
71	18x	rd <- aion::as_fixed(data_clean$dates_ref)
72		} else {
73	2x	rd <- aion::fixed(data_clean$dates_ref, calendar = calendar)
74		}
75
76		## Gaussian multiple linear regression model
77	20x	contexts <- data.frame(
78	20x	date = rd,
79	20x	row_coord[!row_coord$.sup, -ncol(row_coord), drop = FALSE]
80		)
81	20x	fit <- stats::lm(date ~ ., data = contexts)
82
83		## FIXME: temporary workaround until aion::fixed() can deal with NA
84	20x	dates_not_na <- which(!is.na(dates))
85	20x	rt <- rep(NA_real_, length(dates))
86	20x	if (is.null(calendar)) {
87	18x	rt[dates_not_na] <- aion::as_fixed(dates[dates_not_na])
88		} else {
89	2x	rt[dates_not_na] <- aion::fixed(dates[dates_not_na], calendar = calendar)
90		}
91	20x	rt <- aion::as_fixed(rt)
92
93	20x	.EventDate(
94	20x	results_CA,
95	20x	dates = rt,
96	20x	model = fit,
97	20x	keep = seq_len(rank)
98		)
99		}
100		)
101
102		clean_total <- function(data_ref, data_sup, dates_ref, dates_sup,
103		total = 5, verbose = TRUE) {
104
105	20x	clean <- TRUE
106	20x	n_rm_col <- n_rm_ref <- n_rm_sup <- 0
107
108	20x	while (clean) {
109	34x	rm_col <- colSums(data_ref) < total
110	34x	n_rm_col <- n_rm_col + sum(rm_col)
111	34x	if (any(rm_col)) {
112	!	data_ref <- data_ref[, !rm_col, drop = FALSE]
113	!	data_sup <- data_sup[, !rm_col, drop = FALSE]
114		}
115
116	34x	rm_row_ref <- rowSums(data_ref) < total
117	34x	n_rm_ref <- n_rm_ref + sum(rm_row_ref)
118	34x	if (any(rm_row_ref)) {
119	14x	data_ref <- data_ref[!rm_row_ref, , drop = FALSE]
120	14x	dates_ref <- dates_ref[!rm_row_ref]
121		}
122
123	34x	rm_row_sup <- rowSums(data_sup) < total
124	34x	n_rm_sup <- n_rm_sup + sum(rm_row_sup)
125	34x	if (any(data_sup)) {
126	!	data_sup <- data_sup[!rm_row_sup, , drop = FALSE]
127	!	dates_sup <- dates_sup[!rm_row_sup]
128		}
129
130	20x	if (!any(rm_col) & !any(rm_row_ref) & !any(rm_row_sup)) clean <- FALSE
131		}
132
133	20x	if (isTRUE(verbose)) {
134	!	if (n_rm_col > 0) {
135	!	msg <- ngettext(n_rm_col,
136	!	"%d column has a grand total of less than %d.",
137	!	"%d columns have a grand total of less than %d.")
138	!	message(sprintf(msg, n_rm_col, total))
139		}
140	!	if (n_rm_ref > 0) {
141	!	msg <- ngettext(n_rm_ref,
142	!	"%d row has a grand total of less than %d.",
143	!	"%d rows have a grand total of less than %d.")
144	!	message(sprintf(msg, n_rm_ref, total))
145		}
146	!	if (n_rm_sup > 0) {
147	!	msg <- ngettext(n_rm_sup,
148	!	"%d supplementary row has a grand total of less than %d.",
149	!	"%d supplementary rows have a grand total of less than %d.")
150	!	message(sprintf(msg, n_rm_sup, total))
151		}
152	!	n_rm <- n_rm_col + n_rm_ref + n_rm_sup
153	!	if (n_rm > 0) {
154	!	msg <- ngettext(n_rm, "It was omitted from the analysis.",
155	!	"They were omitted from the analysis.")
156	!	message(msg)
157		}
158		}
159
160	20x	if (all(is.na(dates_ref))) {
161	!	stop(tr_("All dates are missing!"), call. = FALSE)
162		}
163
164	20x	list(
165	20x	data_ref = data_ref,
166	20x	data_sup = data_sup,
167	20x	dates_ref = dates_ref,
168	20x	dates_sup = dates_sup
169		)
170		}
171
172		# Event ========================================================================
173		#' @export
174		#' @rdname predict_event
175		#' @aliases predict_event,EventDate,missing-method
176		setMethod(
177		f = "predict_event",
178		signature = c(object = "EventDate", data = "missing"),
179		definition = function(object, margin = 1, level = 0.95,
180		calendar = get_calendar()) {
181	12x	data <- dimensio::get_coordinates(object, margin = margin)
182	12x	data <- as.matrix(data[, -ncol(data), drop = FALSE])
183
184	12x	.predict_event(object, data, level = level, calendar = calendar)
185		}
186		)
187
188		#' @export
189		#' @rdname predict_event
190		#' @aliases predict_event,EventDate,matrix-method
191		setMethod(
192		f = "predict_event",
193		signature = c(object = "EventDate", data = "matrix"),
194		definition = function(object, data, margin = 1, level = 0.95,
195		calendar = get_calendar()) {
196		## Correspondence analysis
197	!	data <- dimensio::predict(object, data, margin = margin)
198
199	!	.predict_event(object, data, level = level, calendar = calendar)
200		}
201		)
202
203		.predict_event <- function(object, data, level = 0.95,
204		calendar = get_calendar()) {
205		## Predict event date
206	12x	fit_model <- object@model
207	12x	ca_event <- compute_event(fit_model, data, level = level)
208
209		# TODO: check predicted dates consistency
210
211		## Convert
212	12x	ca_event <- as.data.frame(ca_event)
213	11x	if (is.null(calendar)) return(ca_event)
214
215	1x	ca_event[] <- lapply(
216	1x	X = ca_event,
217	1x	FUN = aion::as_year,
218	1x	calendar = calendar,
219	1x	decimal = TRUE
220		)
221	1x	ca_event
222		}
223
224		#' @export
225		#' @rdname density_event
226		#' @aliases density_event,EventDate-method
227		setMethod(
228		f = "density_event",
229		signature = c(object = "EventDate"),
230		definition = function(object, dates = NULL, calendar = NULL, n = 500, ...) {
231		## Get data
232	4x	rows <- predict_event(object, margin = 1, calendar = NULL)
233	4x	row_dates <- rows$date
234	4x	row_errors <- rows$error
235
236	4x	if (is.null(dates)) {
237	4x	dates <- seq(
238	4x	from = min(rows$lower, na.rm = TRUE) * 0.96,
239	4x	to = max(rows$upper, na.rm = TRUE) * 1.04,
240	4x	length.out = n
241		)
242		}
243	4x	if (is.null(calendar)) {
244	4x	dates <- aion::as_fixed(dates)
245		} else {
246	!	dates <- aion::fixed(dates, calendar = calendar)
247		}
248
249		## Density estimation
250	4x	date_event <- mapply(
251	4x	FUN = stats::dnorm,
252	4x	mean = row_dates,
253	4x	sd = row_errors,
254	4x	MoreArgs = list(x = dates),
255	4x	SIMPLIFY = TRUE
256		)
257
258	4x	colnames(date_event) <- rownames(rows)
259	4x	aion::series(object = date_event, time = dates)
260		}
261		)
262
263		# Accumulation =================================================================
264		#' @export
265		#' @rdname predict_event
266		#' @aliases predict_accumulation,EventDate,missing-method
267		setMethod(
268		f = "predict_accumulation",
269		signature = c(object = "EventDate", data = "missing"),
270		definition = function(object, level = 0.95,
271		calendar = get_calendar()) {
272	1x	data <- object@data
273	1x	methods::callGeneric(object, data = data, level = level, calendar = calendar)
274		}
275		)
276
277		#' @export
278		#' @rdname predict_event
279		#' @aliases predict_accumulation,EventDate,matrix-method
280		setMethod(
281		f = "predict_accumulation",
282		signature = c(object = "EventDate", data = "matrix"),
283		definition = function(object, data, level = 0.95,
284		calendar = get_calendar()) {
285		## Predict event date
286	1x	col_event <- predict_event(object, margin = 2, calendar = NULL)
287
288		## Accumulation time point estimate
289	1x	date_range <- seq(
290	1x	from = min(col_event$lower, na.rm = TRUE),
291	1x	to = max(col_event$upper, na.rm = TRUE),
292	1x	length.out = 1000
293		)
294	1x	dens <- mapply(
295	1x	FUN = stats::dnorm,
296	1x	mean = col_event$date,
297	1x	sd = col_event$error,
298	1x	MoreArgs = list(date_range),
299	1x	SIMPLIFY = TRUE
300		)
301	1x	acc <- apply(
302	1x	X = data,
303	1x	MARGIN = 1,
304	1x	FUN = function(x, density) {
305	420x	colSums(t(density) * as.numeric(x))
306		},
307	1x	density = dens
308		)
309
310	1x	alpha <- (1 - level) / 2
311	1x	probs <- c(0.5, alpha, 1 - alpha)
312	1x	quant <- quantile_density(x = date_range, y = acc, probs = probs)
313
314		## Convert
315	1x	quant <- as.data.frame(quant)
316	1x	colnames(quant) <- c("date", "lower", "upper")
317	1x	if (is.null(calendar)) return(quant)
318
319	!	quant[] <- lapply(
320	!	X = quant,
321	!	FUN = aion::as_year,
322	!	calendar = calendar,
323	!	decimal = TRUE
324		)
325	!	quant
326		}
327		)
328
329		#' @export
330		#' @rdname density_event
331		#' @aliases density_accumulation,EventDate-method
332		setMethod(
333		f = "density_accumulation",
334		signature = c(object = "EventDate"),
335		definition = function(object, dates = NULL, calendar = NULL,
336		type = c("activity", "tempo"), n = 500, ...) {
337		## Validation
338	4x	type <- match.arg(type, several.ok = FALSE)
339
340		## Get data
341	4x	columns <- predict_event(object, margin = 2, calendar = NULL)
342	4x	col_dates <- columns$date
343	4x	col_errors <- columns$error
344
345	4x	if (is.null(dates)) {
346	!	dates <- seq(
347	!	from = min(columns$lower, na.rm = TRUE) * 0.96,
348	!	to = max(columns$upper, na.rm = TRUE) * 1.04,
349	!	length.out = n
350		)
351		}
352	4x	if (is.null(calendar)) {
353	4x	dates <- aion::as_fixed(dates)
354		} else {
355	!	dates <- aion::fixed(dates, calendar = calendar)
356		}
357
358		## Weighted sum of the fabric dates
359	4x	counts <- dimensio::get_data(object)
360	4x	freq <- counts / rowSums(counts)
361
362		## Tempo vs activity plot
363	4x	fun <- switch(
364	4x	type,
365	4x	activity = stats::dnorm,
366	4x	tempo = stats::pnorm
367		)
368
369		## Density estimation
370	4x	col_density <- mapply(
371	4x	FUN = fun,
372	4x	mean = col_dates,
373	4x	sd = col_errors,
374	4x	MoreArgs = list(dates),
375	4x	SIMPLIFY = TRUE
376		)
377
378	4x	date_acc <- apply(
379	4x	X = freq,
380	4x	MARGIN = 1,
381	4x	FUN = function(x, density) {
382	1680x	colSums(t(density) * as.numeric(x))
383		},
384	4x	density = col_density
385		)
386
387	4x	aion::series(object = date_acc, time = dates)
388		}
389		)
390
391		# Helpers ======================================================================
392		#' Predict event dates
393		#'
394		#' @param fit A \code{\link[stats:lm]{multiple linear model}}.
395		#' @param data A [`numeric`] matrix giving the coordinates in CA space.
396		#' @param level A length-one [`numeric`] vector giving the confidence level.
397		#' @return
398		#' A four columns [`numeric`] matrix giving the predicted event dates,
399		#' the corresponding confidence interval boundaries and the standard error of
400		#' the predicted dates.
401		#' @author N. Frerebeau
402		#' @keywords internal
403		#' @noRd
404		compute_event <- function(fit, data, level, error = TRUE) {
405	13x	data <- as.data.frame(data)
406
407	13x	date_predict <- stats::predict.lm(
408	13x	object = fit,
409	13x	newdata = data,
410	13x	se.fit = TRUE,
411	13x	interval = "confidence",
412	13x	level = level
413		)
414
415		# Three columns matrix: predicted value + CI boundaries
416	13x	results <- as.data.frame(date_predict$fit)
417	13x	rownames(results) <- rownames(data)
418	13x	colnames(results) <- c("date", "lower", "upper")
419	12x	if (error) results$error <- date_predict$se.fit
420
421	13x	results
422		}

1		# REFINE DATE MODEL
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname resample_event
7		#' @aliases jackknife,EventDate-method
8		setMethod(
9		f = "jackknife",
10		signature = c(object = "EventDate"),
11		definition = function(object, level = 0.95,
12		calendar = get_calendar(),
13		progress = getOption("kairos.progress"),
14		verbose = getOption("kairos.verbose"), ...) {
15		## Get data
16	1x	fit_model <- object@model
17	1x	fit_dates <- time(object)
18	1x	fit_data <- dimensio::get_data(object)
19	1x	fit_dim <- object@keep
20
21		## Compute event dates
22	1x	event <- predict_event(object, level = level, calendar = calendar)
23
24		## Jackknife model coefficients
25	1x	jack_values <- compute_date_jack(
26	1x	x = fit_data,
27	1x	dates = fit_dates,
28	1x	rank = length(fit_dim),
29	1x	progress = progress,
30	1x	verbose = verbose
31		)
32	1x	jack_coef <- colMeans(jack_values)
33
34		## Change lm coefficients (UGLY)
35	1x	fit_model$coefficients <- jack_coef[c(1, fit_dim + 1)]
36
37		## Predict event date for each context
38	1x	row_coord <- dimensio::get_coordinates(object, margin = 1)
39	1x	jack_event <- compute_event(fit_model, row_coord, level = level, error = FALSE)
40
41	1x	results <- as.data.frame(jack_event)
42	1x	if (!is.null(calendar)) {
43	1x	results[] <- lapply(
44	1x	X = results,
45	1x	FUN = aion::as_year,
46	1x	calendar = calendar,
47	1x	decimal = TRUE
48		)
49		}
50	1x	results
51		}
52		)
53
54		#' Jackknife Fabrics
55		#'
56		#' Compute date event jackknifed statistics for each replicated sample.
57		#' @param x A [`numeric`] matrix of count data.
58		#' @param dates A [`numeric`] vector of known dates.
59		#' @param rank A [`numeric`] value.
60		#' @param progress A [`logical`] scalar: should a progress bar be displayed?
61		#' @param ... Currently not used].
62		#' @return A `numeric` [`matrix`] of linear model coefficients.
63		#' @author N. Frerebeau
64		#' @keywords internal
65		#' @noRd
66		compute_date_jack <- function(x, dates, rank = 10,
67		progress = FALSE, verbose = FALSE, ...) {
68	1x	m <- ncol(x)
69	1x	k <- seq_len(m)
70	1x	jack <- vector(mode = "list", length = m)
71
72	1x	progress_bar <- interactive() && isTRUE(progress)
73	!	if (progress_bar) pbar <- utils::txtProgressBar(max = m, style = 3)
74
75	1x	for (j in k) {
76	18x	counts <- x[, -j, drop = FALSE]
77	18x	model <- event(counts, dates = dates, calendar = NULL, rank = rank, verbose = verbose)
78	18x	jack[[j]] <- coef(model, calendar = NULL) # Get model coefficients
79	!	if (progress_bar) utils::setTxtProgressBar(pbar, j)
80		}
81
82	!	if (progress_bar) close(pbar)
83	1x	do.call(rbind, jack)
84		}
85
86		#' @export
87		#' @rdname resample_event
88		#' @aliases bootstrap,EventDate-method
89		setMethod(
90		f = "bootstrap",
91		signature = c(object = "EventDate"),
92		definition = function(object, level = 0.95, probs = c(0.05, 0.95), n = 1000,
93		calendar = get_calendar(),
94		progress = getOption("kairos.progress"), ...) {
95		## Get data
96	!	fit_model <- object@model
97	!	fit_dim <- object@keep
98
99		## Partial bootstrap CA
100		## /!\ Be careful: EventDate inherits from CA
101	!	ca_boot <- methods::callNextMethod(object, n = n)
102	!	ca_rep_row <- dimensio::get_replications(ca_boot, margin = 1)
103
104		## Bootstrap assemblages
105	!	event_rows <- apply(
106	!	X = ca_rep_row,
107	!	MARGIN = 1,
108	!	FUN = function(x, axes, model, level, probs) {
109	!	compute_date_boot(t(x), axes, model, level, probs)
110		},
111	!	axes = fit_dim,
112	!	model = fit_model,
113	!	level = level,
114	!	probs = probs
115		)
116
117	!	results <- as.data.frame(t(event_rows))
118	!	if (is.null(calendar)) return(results)
119
120	!	results[] <- lapply(
121	!	X = results,
122	!	FUN = aion::as_year,
123	!	calendar = calendar,
124	!	decimal = TRUE
125		)
126	!	results
127		}
128		)
129
130		#' Bootstrap Resampling of Assemblages
131		#'
132		#' Computes date event bootstraped statistics for each replicated sample.
133		#' @param x A [`numeric`] matrix of bootstrap replicates.
134		#' @param axes A [`numeric`] vector giving the subscripts of the CA components
135		#' to keep.
136		#' @param model An object of class [`lm`][stats::lm()].
137		#' @param level A length-one [`numeric`] vector giving the confidence level.
138		#' @param probs A [`numeric`] vector of probabilities with values in \eqn{[0,1]}
139		#' (see [stats:quantile()]).
140		#' @return A [`numeric`] vector with the following elements:
141		#' \describe{
142		#' \item{`min`}{Minimum value.}
143		#' \item{`mean`}{Mean value (event date).}
144		#' \item{`max`}{Maximum value.}
145		#' \item{`Q05`}{Sample quantile to 0.05 probability.}
146		#' \item{`Q95`}{Sample quantile to 0.95 probability.}
147		#' }
148		#' @author N. Frerebeau
149		#' @keywords internal
150		#' @noRd
151		compute_date_boot <- function(x, axes, model, level, probs = c(0.05, 0.95)) {
152		## Remove missing values
153	!	coords <- stats::na.omit(x[, axes])
154
155		## Gaussian multiple linear regression model
156	!	event <- compute_event(model, coords, level)[, 1]
157	!	Q <- stats::quantile(event, probs = probs, names = FALSE)
158
159	!	distrib <- c(min(event), mean(event), max(event), Q)
160	!	quant <- paste0("Q", round(probs * 100, 0))
161	!	names(distrib) <- c("min", "mean", "max", quant)
162	!	return(distrib)
163		}

1		# PLOT DATES
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @method plot EventDate
7		plot.EventDate <- function(x, type = c("activity", "tempo"), event = FALSE,
8		calendar = get_calendar(),
9		select = 1, n = 500, eps = 1e-09,
10		col.accumulation = "black", col.event = "red",
11		flip = FALSE, ncol = NULL,
12		xlab = NULL, ylab = NULL,
13		main = NULL, sub = NULL,
14		ann = graphics::par("ann"), axes = TRUE,
15		frame.plot = axes, ...) {
16		## Validation
17	4x	type <- match.arg(type, several.ok = FALSE)
18
19		## Event date
20	4x	date_event <- density_event(x, n = n)
21	4x	date_range <- aion::time(date_event)
22	4x	cases <- colnames(date_event)
23
24		## Accumulation time
25	4x	date_acc <- density_accumulation(x, dates = date_range, type = type, n = n)
26
27		## Selection
28	!	if (is.null(select)) index <- seq_along(cases)
29	4x	else if (is.character(select)) index <- which(cases %in% select)
30	!	else index <- as.numeric(select)
31
32	4x	k <- length(index)
33	1x	if (k == 0) stop(tr_("Wrong selection."), call. = FALSE)
34
35	3x	if (type != "activity" \|\| !event) {
36	2x	date_event <- array(data = NA, dim = list(n, k, 1))
37		}
38
39	3x	date_event <- date_event[, index, , drop = FALSE]
40	3x	date_acc <- date_acc[, index, , drop = FALSE]
41
42		## Cleaning
43		# date_acc[date_acc < eps] <- NA
44	3x	date_event[date_event < eps] <- NA
45	3x	na_event <- apply(date_event, 1, function(x) all(is.na(x)))
46	3x	na_acc <- apply(date_acc, 1, function(x) all(is.na(x)))
47	3x	date_drop <- na_event & na_acc
48
49		## Time series
50	3x	ts <- array(data = c(date_acc, date_event), dim = c(n, k, 2))
51	3x	dimnames(ts) <- list(NULL, cases[index], c("accumulation", "event"))
52	3x	ts <- aion::series(object = ts[!date_drop, , , drop = FALSE],
53	3x	time = aion::as_fixed(date_range[!date_drop]))
54
55	3x	panel <- switch(
56	3x	type,
57	3x	activity = function(x, y, ...) {
58	8x	graphics::polygon(x = c(x, rev(x)), y = c(y, rep(0, length(y))),
59	8x	border = NA, ...)
60	8x	graphics::lines(x, y, col = "black", lty = list(...)$lty)
61		},
62	3x	tempo = function(x, y, ...) graphics::lines(x, y, col = "black", lty = 1)
63		)
64
65	3x	col <- c(col.accumulation, grDevices::adjustcolor(col.event, alpha = 0.5))
66	3x	aion::plot(ts, panel = panel, calendar = calendar,
67	3x	flip = flip, ncol = ncol, xlab = xlab, ylab = ylab,
68	3x	main = main, sub = sub, ann = ann, axes = axes,
69	3x	frame.plot = frame.plot,
70	3x	col = col, lty = c(0, 0))
71
72	3x	invisible(x)
73		}
74
75		#' @export
76		#' @rdname plot_event
77		#' @aliases plot,EventDate,missing-method
78		setMethod("plot", c(x = "EventDate", y = "missing"), plot.EventDate)

1		# MUTATORS
2		#' @include AllClasses.R
3		NULL
4
5		# Getters ======================================================================
6		#' @export
7		#' @rdname series
8		#' @aliases time,EventDate-method
9		setMethod(
10		f = "time",
11		signature = "EventDate",
12		definition = function(x, calendar = NULL) {
13	1x	z <- x@dates
14	1x	if (is.null(calendar)) return(z)
15	!	aion::as_year(z, calendar = calendar, decimal = TRUE)
16		}
17		)
18
19		#' @export
20		#' @rdname series
21		#' @aliases span,AoristicSum-method
22		setMethod(
23		f = "span",
24		signature = "AoristicSum",
25		definition = function(x, calendar = NULL) {
26	1x	z <- x@span
27	1x	if (is.null(calendar)) return(unclass(z))
28	!	aion::as_year(z, calendar = calendar, shift = FALSE)
29		}
30		)
31
32		#' @export
33		#' @rdname mutators
34		#' @aliases weights,AoristicSum-method
35	2x	setMethod("weights", "AoristicSum", function(object, ...) object@p)
36
37		#' @export
38		#' @rdname mutators
39		#' @aliases weights,CountApportion-method
40	!	setMethod("weights", "CountApportion", function(object, ...) object@p)
41
42		# Setters ======================================================================

1		# PLOT LINE
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		# data.frame ===================================================================
6		#' @export
7		#' @rdname plot_time
8		#' @aliases plot_time,data.frame,numeric-method
9		setMethod(
10		f = "plot_time",
11		signature = c(object = "data.frame", dates = "numeric"),
12		definition = function(object, dates,
13		calendar = get_calendar(), ...) {
14	1x	object <- data.matrix(object)
15	1x	methods::callGeneric(object, dates = dates, calendar = calendar, ...)
16	1x	invisible(object)
17		}
18		)
19
20		# matrix =======================================================================
21		#' @export
22		#' @rdname plot_time
23		#' @aliases plot_time,matrix,numeric-method
24		setMethod(
25		f = "plot_time",
26		signature = c(object = "matrix", dates = "numeric"),
27		definition = function(object, dates,
28		calendar = get_calendar(), ...) {
29		## Validation
30	1x	arkhe::assert_type(dates, "numeric")
31	1x	arkhe::assert_length(dates, nrow(object))
32
33		## Convert to rata die
34	1x	if (is.null(calendar)) {
35	!	dates <- aion::as_fixed(dates)
36		} else {
37	1x	dates <- aion::fixed(dates, calendar = calendar)
38		}
39
40		## Prepare data
41	1x	ts <- aion::series(object = object, time = dates)
42
43		## Plot
44	1x	aion::plot(x = ts, calendar = calendar, ...)
45
46	1x	invisible(object)
47		}
48		)

1		# COERCE
2		#' @include AllGenerics.R
3		NULL
4
5		# To data.frame ================================================================
6		#' @export
7		#' @method as.data.frame MeanDate
8		as.data.frame.MeanDate <- function(x, ..., calendar = get_calendar()) {
9	1x	data.frame(
10	1x	sample = rownames(x) %\|\|% paste0("S", seq_len(nrow(x))),
11	1x	time = aion::time(x, calendar = calendar)
12		)
13		}
14
15		#' @export
16		#' @rdname data.frame
17		#' @aliases as.data.frame,MeanDate-method
18		setMethod("as.data.frame", "MeanDate", as.data.frame.MeanDate)
19
20		#' @export
21		#' @method as.data.frame AoristicSum
22		as.data.frame.AoristicSum <- function(x, ..., calendar = get_calendar()) {
23	1x	block_start <- utils::head(x@breaks, -1)
24	1x	block_end <- utils::tail(x@breaks, -1)
25
26	1x	aorist <- x[, , 1, drop = TRUE]
27	1x	dim(aorist) <- c(nrow(x), ncol(x))
28	1x	colnames(aorist) <- colnames(x) %\|\|% paste0("S", seq_len(ncol(x)))
29
30	1x	data.frame(
31	1x	start = aion::as_year(block_start, calendar = calendar),
32	1x	end = aion::as_year(block_end, calendar = calendar),
33	1x	aorist
34		)
35		}
36
37		#' @export
38		#' @rdname data.frame
39		#' @aliases as.data.frame,AoristicSum-method
40		setMethod("as.data.frame", "AoristicSum", as.data.frame.AoristicSum)
41
42		#' @method as.data.frame IncrementTest
43		#' @export
44		as.data.frame.IncrementTest <- function(x, row.names = NULL, optional = FALSE, ...) {
45	1x	data.frame(
46	1x	t = x@statistic,
47	1x	p.value = x@p_value,
48	1x	row.names = colnames(x),
49		...
50		)
51		}
52
53		#' @export
54		#' @rdname data.frame
55		#' @aliases as.data.frame,IncrementTest-method
56		setMethod("as.data.frame", "IncrementTest", as.data.frame.IncrementTest)

1		# REFINE MATRIX SERIATION
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname refine
7		#' @aliases refine,AveragePermutationOrder-method
8		setMethod(
9		f = "refine",
10		signature = c(object = "AveragePermutationOrder"),
11		definition = function(object, cutoff, margin = 1, axes = 1, n = 30, ...) {
12		## Partial bootstrap CA
13		## /!\ Be careful: AveragePermutationOrder inherits from CA
14	!	object <- dimensio::bootstrap(object, n = n)
15	!	methods::callGeneric(object, cutoff = cutoff, margin = margin, axes = axes, ...)
16		}
17		)
18
19		#' @export
20		#' @rdname refine
21		#' @aliases refine,BootstrapCA-method
22		setMethod(
23		f = "refine",
24		signature = c(object = "BootstrapCA"),
25		definition = function(object, cutoff, margin = 1, axes = 1, ...) {
26		## Validation
27	!	arkhe::assert_function(cutoff)
28	!	arkhe::assert_length(margin, 1)
29	!	arkhe::assert_length(axes, 1)
30
31		## Get data
32	!	ca_rep <- dimensio::get_replications(object, margin = margin)
33
34		## Compute convex hull
35	!	hull <- apply(
36	!	X = ca_rep,
37	!	MARGIN = 1,
38	!	FUN = function(x, axes) compute_chull(t(x), c(1, 2)),
39	!	axes = axes
40		)
41
42		## Get convex hull maximal dimension length for each sample
43	!	len <- vapply(
44	!	X = hull,
45	!	FUN = function(x) max(stats::dist(x, method = "euclidean")),
46	!	FUN.VALUE = double(1)
47		)
48
49		## Get cutoff values
50	!	limit <- cutoff(len)
51
52		## Samples to be excluded
53	!	sup <- which(len >= limit)
54
55	!	list(
56	!	length = len,
57	!	cutoff = limit,
58	!	exclude = which(len >= limit),
59	!	margin = as.integer(margin)
60		)
61		}
62		)
63
64		#' Convex Hull of CA Coordinates
65		#'
66		#' Compute convex hull area for each replicated sample
67		#' @param x A [`numeric`] matrix of bootstrap replicates.
68		#' @param axes A length-two [`numeric`] vector giving the subscripts
69		#' of the CA components to use.
70		#' @return A [`matrix`].
71		#' @author N. Frerebeau
72		#' @keywords internal
73		#' @noRd
74		compute_chull <- function(x, axes) {
75		# Remove missing values
76	!	clean <- stats::na.omit(x[, axes])
77		# Get convex hull coordinates
78	!	points <- grDevices::chull(clean)
79		# Repeat first point
80	!	hull <- clean[c(points, points[1]), axes]
81	!	colnames(hull) <- c("x", "y")
82
83	!	return(hull)
84		}