tabula coverage - 93.32%

Files
Source

#' @include AllGenerics.R AllClasses.R
NULL

# Shannon ======================================================================
variance_shannon <- function(x, base = exp(1), na.rm = FALSE, ...) {
  ## Validation
  x <- x[x > 0] # Remove unobserved species
  if (na.rm) x <- stats::na.omit(x) # Remove NAs
  if (anyNA(x)) return(NA)

  N <- sum(x)
  S <- length(x)
  p <- x / N

  a <- sum(p * (log(p, base = base))^2)
  b <- sum(p * log(p, base = base))^2

  var <- ((a - b) / N) + ((S - 1) / (2 * N^2))
  var
}

#' @export
#' @rdname test
#' @aliases test_shannon,numeric,numeric-method
setMethod(
  f = "test_shannon",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, ...) {
    ## Validation
    arkhe::assert_length(y, length(x))

    ## Calculate the number of individuals
    Nx <- sum(x, na.rm = TRUE)
    Ny <- sum(y, na.rm = TRUE)

    ## Calculate Shannon diversity
    ## See PAST documentation, p. 208
    d <- (sum(x > 0) - 1) / (2 * sum(x))
    Hx <- index_shannon(x, ...)
    Hy <- index_shannon(y, ...)

    ## Calculate Shannon variance
    Vx <- variance_shannon(x, ...)
    Vy <- variance_shannon(y, ...)

    ## t test statistic
    tt <- (Hx - Hy) / sqrt(Vx + Vy)

    ## Degrees of freedom
    df <- (Vx + Vy)^2 / sum(c(Vx, Vy)^2 / c(Nx, Ny))

    ## p value
    p <- 2 * (1 - stats::pt(q = abs(tt), df = df))

    list(
      statistic = tt,
      parameter = df,
      p.value = p
    )
  }
)

#' @export
#' @describeIn test Produces two sided pairwise comparisons.
#' @aliases test_shannon,matrix,missing-method
setMethod(
  f = "test_shannon",
  signature = c(x = "matrix", y = "missing"),
  definition = function(x, adjust = "holm", ...) {
    ## Get the names of the assemblages
    row_names <- rownames(x)
    if (length(row_names) != 0) {
      row_names <- factor(row_names, levels = unique(row_names))
    } else {
      row_names <- factor(seq_len(nrow(x)))
    }

    ## Compute t test
    compare <- function(i, j) {
      test_shannon(x[i, ], x[j, ], ...)$p.value
    }

    result <- stats::pairwise.table(
      compare.levels = compare,
      level.names = row_names,
      p.adjust.method = adjust
    )
    result
  }
)

#' @export
#' @describeIn test Produces two sided pairwise comparisons.
#' @aliases test_shannon,data.frame,missing-method
setMethod(
  f = "test_shannon",
  signature = c(x = "data.frame", y = "missing"),
  definition = function(x, adjust = "holm", ...) {
    x <- data.matrix(x)
    methods::callGeneric(x, adjust = adjust, ...)
  }
)

# Simpson ======================================================================
variance_simpson <- function(x, na.rm = FALSE, ...) {
  ## Validation
  x <- x[x > 0] # Remove unobserved species
  if (na.rm) x <- stats::na.omit(x) # Remove NAs
  if (anyNA(x)) return(NA)

  N <- sum(x)
  S <- length(x)
  p <- x / N

  a <- 4 * N * (N - 1) * (N - 2) * sum(p^3)
  b <- 2 * N * (N - 1) * sum(p^2)
  c <- 2 * N * (N - 1) * (2 * N - 3) * sum(p^2)^2

  var <- (a + b - c) / (N^2 * (N - 1)^2)
  var
}

#' @export
#' @rdname test
#' @aliases test_simpson,numeric,numeric-method
setMethod(
  f = "test_simpson",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, adjust = "holm", ...) {
    ## Validation
    arkhe::assert_length(y, length(x))

    ## Calculate the number of individuals
    Nx <- sum(x, na.rm = TRUE)
    Ny <- sum(y, na.rm = TRUE)

    ## Calculate Shannon diversity
    Hx <- index_simpson(x, ...)
    Hy <- index_simpson(y, ...)

    ## Calculate Shannon variance
    Vx <- variance_simpson(x, ...)
    Vy <- variance_simpson(y, ...)

    ## t test statistic
    tt <- (Hx - Hy) / sqrt(Vx + Vy)

    ## Degrees of freedom
    df <- (Vx + Vy)^2 / sum(c(Vx, Vy)^2 / c(Nx, Ny))

    ## p value
    p <- 2 * (1 - stats::pt(q = abs(tt), df = df))

    list(
      statistic = tt,
      parameter = df,
      p.value = p
    )
  }
)

#' @export
#' @describeIn test Produces two sided pairwise comparisons.
#' @aliases test_simpson,matrix,missing-method
setMethod(
  f = "test_simpson",
  signature = c(x = "matrix", y = "missing"),
  definition = function(x, adjust = "holm", ...) {
    ## Get the names of the assemblages
    row_names <- rownames(x)
    if (length(row_names) != 0) {
      row_names <- factor(row_names, levels = unique(row_names))
    } else {
      row_names <- factor(seq_len(nrow(x)))
    }

    ## Compute t test
    compare <- function(i, j) {
      test_simpson(x[i, ], x[j, ])$p.value
    }

    result <- stats::pairwise.table(
      compare.levels = compare,
      level.names = row_names,
      p.adjust.method = adjust
    )
    result
  }
)

#' @export
#' @describeIn test Produces two sided pairwise comparisons.
#' @aliases test_simpson,data.frame,missing-method
setMethod(
  f = "test_simpson",
  signature = c(x = "data.frame", y = "missing"),
  definition = function(x, adjust = "holm", ...) {
    x <- data.matrix(x)
    methods::callGeneric(x, adjust = adjust)
  }
)

# MATRIX PLOT

# Plot =========================================================================
#' Matrix Plot
#'
#' @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 panel A [`function`] in the form `function(x, y, z, color, ...)`
#'  which gives the action to be carried out in each panel of the display.
#' @param diag A [`logical`] scalar indicating whether the diagonal of the
#'  matrix should be plotted. Only used if `object` is a symmetric matrix.
#' @param upper A [`logical`] scalar indicating whether the upper triangle of
#'  the matrix should be plotted. Only used if `object` is a symmetric matrix.
#' @param lower A [`logical`] scalar indicating whether the lower triangle of
#'  the matrix should be plotted. Only used if `object` is a symmetric matrix.
#' @param freq A [`logical`] scalar indicating whether conditional proportions
#'  given `margins` should be used (i.e. entries of `object`, divided by the
#'  appropriate marginal sums).
#' @param margin An [`integer`] vector giving the margins to split by:
#'  `1` indicates individuals/rows (the default), `2` indicates
#'  variables/columns. Only used if `freq` is `TRUE`.
#' @param scale A [`logical`] scalar indicating whether data should be rescaled
#'  to \eqn{[-1,1]}. Only used if `freq` if `FALSE`.
#' @param drop_zero A [`logical`] scalar: should zeros be discarded?
#' @param color A vector of colors or a [`function`] that when called with a
#'  single argument (an integer specifying the number of colors) returns a
#'  vector of colors.
#' @param midpoint A [`numeric`] value specifying the data midpoint.
#' @param axes A [`logical`] scalar: should axes be drawn on the plot? It will
#'  omit labels where they would abut or overlap previously drawn labels.
#' @param legend A [`logical`] scalar: should a legend be displayed?
#' @param asp A length-one [`numeric`] vector, giving the aspect ratio
#'  \eqn{y/x}.
#' @param ... Further arguments to be passed to `panel`.
#' @keywords internal
plot_matrix <- function(object, panel, diag = TRUE, upper = TRUE, lower = TRUE,
                        freq = FALSE, margin = 1, scale = TRUE, drop_zero = TRUE,
                        color = graphics::par("fg"), midpoint = NULL,
                        axes = TRUE, legend = TRUE, asp = 1, ...) {
  ## Validation
  if (is_incidence(object)) {
    object[] <- as.numeric(object)
    legend <- FALSE
  }

  ## Prepare data
  n <- nrow(object)
  m <- ncol(object)
  seq_row <- rev(seq_len(n))
  seq_col <- seq_len(m)
  lab_row <- rownames(object) %||% seq_row
  lab_col <- colnames(object) %||% seq_col

  data <- prepare(object, diag = diag, upper = upper, lower = lower,
                  freq = freq, margin = margin, scale = scale,
                  drop_zero = drop_zero,
                  palette = color, midpoint = midpoint)

  ## Graphical parameters
  cex.axis <- graphics::par("cex.axis")
  col.axis <- graphics::par("col.axis")
  font.axis <- graphics::par("font.axis")

  ## Save and restore
  d <- inch2line("M", cex = cex.axis)
  old_par <- graphics::par("mar", "plt")
  on.exit(graphics::par(old_par))

  mar_left <- inch2line(lab_row, cex = cex.axis)
  mar_top <- inch2line(lab_col, cex = cex.axis)
  mar_right <- if (legend) inch2line("999%", cex = cex.axis) else d
  graphics::par(mar = c(d, mar_left, mar_top, mar_right))

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

  ## Squish plotting area
  pin <- graphics::par("pin")
  plt <- graphics::par("plt")

  ## Add horizontal space for the legend (5%)
  legend_width <- if (legend) max(1, m / 20) else 0

  if (isTRUE(asp)) asp <- 1
  if (!isFALSE(asp) && !is.na(asp)) {
    aspect_ratio <- n / (m + legend_width)
    pin_y <- pin[1] * aspect_ratio * asp

    if (pin_y < pin[2]) {
      ## Squish vertically
      graphics::par(pin = c(pin[1], pin_y))
      graphics::par(plt = c(plt[1:2], graphics::par('plt')[3:4]))
    } else {
      ## Squish horizontally
      pin_x <- pin[2] / aspect_ratio / asp
      graphics::par(pin = c(pin_x, pin[2]))
      graphics::par(plt = c(graphics::par('plt')[1:2], plt[3:4]))
    }
  }

  ## Set plotting coordinates
  xlim <- c(0, m + legend_width) + 0.5
  ylim <- c(0, n) + 0.5
  graphics::plot.window(xlim = xlim, ylim = ylim, xaxs = "i", yaxs = "i", asp = asp)

  ## Plot
  panel(x = data$x, y = data$y, z = data$scaled, color = data$color, ...)

  ## Construct axis
  if (axes) {
    graphics::axis(side = 2, at = seq_row, labels = lab_row, las = 2,
                   lty = 0, cex.axis = cex.axis, col.axis = col.axis,
                   font.axis = font.axis)
    graphics::axis(side = 3, at = seq_col, labels = lab_col, las = 2,
                   lty = 0, cex.axis = cex.axis, col.axis = col.axis,
                   font.axis = font.axis)
  }

  ## Legend
  if (legend) {
    lgd <- attr(data, "legend")
    legend_gradient(x = m, y = n, labels = lgd$labels,
                    at = lgd$at, width = legend_width, col = lgd$colors)
  }
}

legend_gradient <- function(x, y, labels, at, width, col) {
  legend_image <- grDevices::as.raster(col)
  legend_y <- (at - min(at)) * y / diff(range(at)) + 0.5

  graphics::rasterImage(legend_image, xleft = x + 1, ybottom = max(legend_y),
                        xright = x + 1 + width, ytop = min(legend_y))
  graphics::segments(x0 = x + 1, y0 = legend_y,
                     x1 = x + 1 + width, y1 = legend_y,
                     col = "white")
  graphics::polygon(x = c(x, x + width, x + width, x) + 1,
                    y = c(0.5, 0.5, max(legend_y), max(legend_y)),
                    col = NA, border = "black")
  graphics::axis(side = 4, at = legend_y, labels = labels, las = 2)
}

# Panels =======================================================================
panel_bertin <- function(x, y, z, color, ..., space = 0.05) {
  y_bottom <- y - 0.5
  y_top <- y - 0.5 + z - space * 2
  y_top[y_top < y_bottom] <- y_bottom[y_top < y_bottom]
  graphics::rect(
    xleft = x - 0.5 + space,
    ybottom = y_bottom,
    xright = x + 0.5 - space,
    ytop = y_top,
    col = color,
    border = "black"
  )
}
panel_matrigraph <- function(x, y, z, ..., reverse = FALSE) {
  pvi <- data.frame(x = x, y = y, z = z)
  pvi$z <- pvi$z / 100
  pvi_plus <- pvi[pvi$z > 1, ]
  pvi_plus$z <- pvi_plus$z - 1
  pvi_plus$z[pvi_plus$z > 1] <- 1

  if (reverse) {
    pvi_minus <- pvi[pvi$z < 1, ]
    pvi_minus$z <- 1 - pvi_minus$z
  } else {
    pvi_minus <- pvi
    pvi_minus$z[pvi_minus$z > 1] <- 1
  }

  pvi_minus$z <- pvi_minus$z * 0.5
  pvi_plus$z <- pvi_plus$z * 0.5

  col_bkg <- if (reverse) "darkgrey" else "white"
  col_minus <- if (reverse) "white" else "darkgrey"

  graphics::rect(xleft = pvi$x - 0.5, ybottom = pvi$y - 0.5,
                 xright = pvi$x + 0.5, ytop = pvi$y + 0.5,
                 col = col_bkg, border = NA)
  graphics::rect(xleft = pvi_minus$x - pvi_minus$z,
                 ybottom = pvi_minus$y - pvi_minus$z,
                 xright = pvi_minus$x + pvi_minus$z,
                 ytop = pvi_minus$y + pvi_minus$z,
                 col = col_minus, border = NA)
  graphics::rect(xleft = pvi_plus$x - pvi_plus$z,
                 ybottom = pvi_plus$y - pvi_plus$z,
                 xright = pvi_plus$x + pvi_plus$z,
                 ytop = pvi_plus$y + pvi_plus$z,
                 col = "black", border = NA)
}
panel_tiles <- function(x, y, color, ...) {
  width <- 0.5
  graphics::rect(
    xleft = x - width,
    ybottom = y - width,
    xright = x + width,
    ytop = y + width,
    col = color,
    border = NA
  )
}
panel_spot <- function(x, y, z, color, type, ...) {
  radius <- abs(z * 0.45)
  for (i in seq_along(x)) {
    circle(x = x[i], y = y[i], radius = radius[i],
                     col = color[i], border = color[i])
    if (type == "ring") {
      circle(x = x[i], y = y[i], radius = 0.45,
                       col = NA, border = "black")
    }
  }
}

# Prepare ======================================================================
#' Prepare Data for Matrix Plot
#'
#' @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 diag A [`logical`] scalar indicating whether the diagonal of the
#'  matrix should be plotted. Only used if `object` is a symmetric matrix.
#' @param upper A [`logical`] scalar indicating whether the upper triangle of
#'  the matrix should be plotted. Only used if `object` is a symmetric matrix.
#' @param lower A [`logical`] scalar indicating whether the lower triangle of
#'  the matrix should be plotted. Only used if `object` is a symmetric matrix.
#' @param freq A [`logical`] scalar indicating whether conditional proportions
#'  given `margins` should be used (i.e. entries of `object`, divided by the
#'  appropriate marginal sums).
#' @param margin An [`integer`] vector giving the margins to split by:
#'  `1` indicates individuals/rows (the default), `2` indicates
#'  variables/columns. Only used if `freq` is `TRUE`.
#' @param scale A [`logical`] scalar indicating whether data should be rescaled
#'  to \eqn{[-1,1]}. Only used if `freq` if `FALSE`.
#' @param drop_zero A [`logical`] scalar: should zeros be discarded?
#' @param palette A vector of colors.
#' @param midpoint A [`numeric`] value specifying the data midpoint.
#' @param ... Currently not used.
#' @return
#'  A long [`data.frame`].
#'  \describe{
#'   \item{`row`}{}
#'   \item{`column`}{}
#'   \item{`x`,`y`}{Tile center coordinates.}
#'   \item{`z`}{Raw values.}
#'   \item{`value`}{}
#'   \item{`scaled`}{}
#'   \item{`color`}{}
#'  }
#' @keywords internal
#' @noRd
prepare <- function(object, diag = TRUE, upper = TRUE, lower = TRUE,
                    freq = FALSE, margin = 1, scale = !freq, drop_zero = FALSE,
                    palette = grDevices::hcl.colors(12, "YlOrBr", rev = TRUE),
                    midpoint = NULL, ...) {
  ## Validation
  if (!arkhe::is_symmetric(object)) {
    diag <- TRUE
    upper <- TRUE
    lower <- TRUE
  }

  ## Coerce to matrix
  object <- as.matrix(object)

  ## Relative frequencies
  val <- if (freq) prop.table(object, margin = margin) else object

  ## Rescale to [-1;1]
  sca <- if (scale) val / max(abs(val), na.rm = TRUE) else val

  val <- as.vector(val)
  sca <- as.vector(sca)
  min_val <- min(val, na.rm = TRUE)
  max_val <- max(val, na.rm = TRUE)

  ## Build a long table
  row <- row(object, as.factor = TRUE)
  col <- col(object, as.factor = TRUE)
  data <- data.frame(
    row = as.vector(row),
    column = as.vector(col),
    x = as.numeric(col),
    y = as.vector(nrow(object) - as.numeric(row) + 1), # Reverse y for plotting
    z = as.vector(object),
    value = val,
    scaled = sca
  )

  ## Map colors
  breaks <- pretty(val, n = 5)
  domain <- range(c(breaks, min_val, max_val))
  midpoint <- if (is.null(midpoint) & min_val < 0 & max_val > 0) 0 else midpoint
  pal <- khroma::palette_color_continuous(colors = palette, domain = domain, midpoint = midpoint)
  data$color <- if (length(palette) == length(val)) palette else pal(val)

  ## Clean data
  if (!upper) data <- data[!upper.tri(object), ]
  if (!lower) data <- data[!lower.tri(object), ]
  if (!diag) data <- data[data$row != data$column, ]
  if (drop_zero) data <- data[data$value != 0, ]

  ## Legend
  attr(data, "legend") <- list(
    labels = if (freq) label_percent(breaks) else breaks,
    at = breaks / max(abs(val), na.rm = TRUE),
    colors = pal(breaks)
  )
  data
}

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

#' @export
#' @rdname plot_ford
#' @aliases plot_ford,matrix-method
setMethod(
  f = "plot_ford",
  signature = signature(object = "matrix"),
  definition = function(object, weights = FALSE, EPPM = FALSE,
                        fill = "darkgrey", border = NA,
                        axes = TRUE, ...) {
    ## Prepare data
    n <- nrow(object)
    m <- ncol(object)
    seq_row <- rev(seq_len(n))
    seq_col <- seq_len(m)
    lab_row <- rownames(object) %||% seq_row
    lab_col <- colnames(object) %||% seq_col

    padding_x <- 0.05
    padding_y <- 0.5 - 0.01
    data <- prepare_ford(object, padding = padding_x)

    ## Graphical parameters
    cex.axis <- graphics::par("cex.axis")
    col.axis <- graphics::par("col.axis")
    font.axis <- graphics::par("font.axis")

    ## Save and restore
    d <- inch2line("M", cex = cex.axis)
    mfrow <- graphics::par("mfrow")
    mar <- graphics::par("mar")
    mar[1] <- 3
    mar[2] <- inch2line(lab_row, cex = cex.axis) + 0.1
    mar[3] <- inch2line(lab_col, cex = cex.axis) + 0.1
    mar[4] <- 0

    old_par <- graphics::par(mfrow = mfrow, mar = mar)
    on.exit(graphics::par(old_par))

    if (weights) {
      graphics::layout(matrix(c(1, 2), nrow = 1), widths = c(m - 1, 1))
    }

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

    ## Set plotting coordinates
    xlim <- range(data$x - data$value, data$x + data$value)
    ylim <- range(data$y) + c(-1, 1) * padding_y
    graphics::plot.window(xlim = xlim, ylim = ylim)

    ## Plot
    graphics::rect(
      xleft = data$x - data$value,
      ybottom = data$y - padding_y,
      xright = data$x + data$value,
      ytop = data$y + padding_y,
      col = fill,
      border = border
    )
    if (EPPM) {
      graphics::rect(
        xleft = data$x - data$eppm,
        ybottom = data$y - padding_y,
        xright = data$x + data$eppm,
        ytop = data$y + padding_y,
        col = "black",
        border = NA
      )
    }

    ## Construct axis
    if (axes) {
      graphics::axis(side = 2, at = seq_row, labels = lab_row, las = 2,
                     lty = 0, cex.axis = cex.axis, col.axis = col.axis,
                     font.axis = font.axis)
      graphics::axis(side = 3, at = unique(data$x), labels = lab_col, las = 2,
                     lty = 0, cex.axis = cex.axis, col.axis = col.axis,
                     font.axis = font.axis)

      x_axis <- data$x[which.max(data$value)]
      graphics::axis(side = 1, at = c(x_axis - 0.2, x_axis + 0.2), labels = FALSE)
      graphics::axis(side = 1, at = x_axis, labels = label_percent(0.2),
                     tick = FALSE)
    }

    if (weights) {
      ## Graphical parameters
      mar[2] <- 0.1
      mar[4] <- 0.1
      graphics::par(mar = mar)

      ## Open new window
      graphics::plot.new()

      ## Set plotting coordinates
      total <- rowSums(object)
      graphics::plot.window(xlim = c(0, max(total) * 1.05), ylim = ylim,
                            xaxs = "i")

      ## Plot
      y <- length(total) - seq_along(total) + 1
      graphics::rect(
        xleft = 0,
        ybottom = y - padding_y,
        xright = total,
        ytop = y + padding_y,
        col = fill,
        border = border
      )

      ## Construct axis
      if (axes) {
        graphics::segments(x0 = 0, y0 = 0, x1 = 0, y = n + 0.5,
                           col = col.axis, lwd = 1)
        graphics::axis(side = 1)
      }
    }

    invisible(object)
  }
)

#' @export
#' @rdname plot_ford
#' @aliases plot_ford,data.frame-method
setMethod(
  f = "plot_ford",
  signature = signature(object = "data.frame"),
  definition = function(object, weights = FALSE, EPPM = FALSE,
                        fill = "darkgrey", border = NA,
                        axes = TRUE, ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, weights = weights, EPPM = EPPM, fill = fill,
                         border = border, axes = axes, ...)
  }
)

#' Prepare data for Ford plot
#' @return A data.frame.
#' @keywords internal
#' @noRd
prepare_ford <- function(x, padding = 0.05) {
  ## EPPM
  EPPM <- eppm(x) / 100

  ## Relative frequencies
  freq <- x / rowSums(x)

  ## Adaptive spacing between columns
  col_max <- apply(X = freq, MARGIN = 2, FUN = max, na.rm = TRUE)
  roll_max <- roll_sum(col_max, n = 2) + padding * max(freq)
  cum_max <- c(0, cumsum(roll_max))

  ## Build a long table
  row <- row(x, as.factor = TRUE)
  col <- col(x, as.factor = TRUE)
  data <- data.frame(
    row = as.vector(row),
    column = as.vector(col),
    value = as.vector(freq),
    eppm = as.vector(EPPM)
  )

  m <- nrow(freq)
  data$x <- rep(cum_max, each = m) + 1
  data$y <- as.vector(m + 1 - as.numeric(row)) # Reverse levels order

  return(data)
}

# STATISTICS
#' @include AllClasses.R AllGenerics.R
NULL

expected <- function(x) {
  m <- nrow(x)
  p <- ncol(x)

  column_total <- matrix(colSums(x), nrow = m, ncol = p, byrow = TRUE)
  row_total <- matrix(rowSums(x), nrow = m, ncol = p, byrow = FALSE)
  grand_total <- sum(x)

  column_total * row_total / grand_total
}

#' @export
#' @rdname resample
#' @aliases resample,numeric-method
setMethod(
  f = "resample",
  signature = c(object = "numeric"),
  definition = function(object, do, n, size = sum(object), ..., f = NULL) {
    ## Validation
    arkhe::assert_count(object)

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

#' Binomial Coefficient
#'
#' Computes the number of `k`-combinations from a given set of `n` elements
#' ("`n` choose `k`").
#' @param n A length-one [`numeric`] vector.
#' @param k A length-one [`numeric`] vector.
#' @details
#'  Ramanujan approximation is used for \eqn{x!} computation if \eqn{x > 170}.
#' @return A length-one [`numeric`] vector.
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
combination <- function(n, k) {
  # Validation
  arkhe::assert_scalar(n, "numeric")
  arkhe::assert_scalar(k, "numeric")

  if (n > 170 | k > 170) {
    ## Ramanujan approximation of x!
    c <- exp(ramanujan(n) - ramanujan(k) - ramanujan(n - k))
  } else {
    c <- factorial(n) / (factorial(k) * factorial(n - k))
  }
  c
}

#' Ramanujan Factorial Approximation
#'
#' @param x A [`numeric`] vector.
#' @return A [`numeric`] vector.
#' @examples
#'  factorial(50)
#'  exp(ramanujan(50))
#' @references
#'  Ramanujan Aiyangar, S. (1988). *The lost notebook and other unpublished
#'  papers*. Berlin: Springer-Verlag.
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
ramanujan <- function(x){
  x * log(x) - x + log(x * (1 + 4 * x * (1 + 2 * x))) / 6 + log(pi) / 2
}

# BETA DIVERSITY
#' @include AllGenerics.R
NULL

# Cody =========================================================================
#' @export
#' @rdname index_cody
#' @aliases index_cody,matrix-method
setMethod(
  f = "index_cody",
  signature = signature(x = "matrix"),
  definition = function(x) {
    x <- x > 0 # presence/absence

    m <- nrow(x)
    first_transect <- x[1, ]
    last_transect <- x[m, ]

    # Number of taxa gained along the transect
    absent <- x[, !first_transect] # Keep only taxa absent in the first transect
    gained <- sum(apply(X = absent, MARGIN = 2, FUN = cumsum)[m, ] > 0)
    # Number of taxa lost
    lost <- sum(!last_transect)

    beta <- (gained + lost) / 2
    beta
  }
)

# Routledge ====================================================================
## Beta R ----------------------------------------------------------------------
#' @export
#' @rdname index_routledge
#' @aliases index_routledge1,matrix-method
setMethod(
  f = "index_routledge1",
  signature = signature(x = "matrix"),
  definition = function(x) {
    x <- x > 0 # presence/absence

    p <- ncol(x)

    # Total number of taxa recorded in the system
    S <- sum(colSums(x) > 0)
    # Works out which pairs of taxa occurs together in at least one sample
    combine <- as.matrix(occurrence(x)) > 0
    r <- sum(combine) / 2

    beta <- S^2 / (2 * r + S) - 1
    return(beta)
  }
)

## Beta I ----------------------------------------------------------------------
#' @export
#' @rdname index_routledge
#' @aliases index_routledge2,matrix-method
setMethod(
  f = "index_routledge2",
  signature = signature(x = "matrix"),
  definition = function(x) {
    x <- x > 0 # presence/absence

    # Number of samples in which each taxa is present
    e <- colSums(x)
    # Taxa richness of each sample
    alpha <- rowSums(x)
    t <- sum(e)

    beta <- log(t) - (1 / t) * sum(e * log(e)) - (1 / t) * sum(alpha * log(alpha))
    return(beta)
  }
)

## Beta E ----------------------------------------------------------------------
#' @export
#' @rdname index_routledge
#' @aliases index_routledge3,matrix-method
setMethod(
  f = "index_routledge3",
  signature = signature(x = "matrix"),
  definition = function(x) {
    I <- index_routledge2(x)
    E <- exp(I)
    return(E)
  }
)

# Whittaker ====================================================================
#' @export
#' @rdname index_whittaker
#' @aliases index_whittaker,matrix-method
setMethod(
  f = "index_whittaker",
  signature = signature(x = "matrix"),
  definition = function(x) {
    x <- x > 0 # presence/absence

    # Total number of taxa recorded in the system
    S <- sum(colSums(x) > 0)
    # Mean taxa richness
    alpha <- mean(rowSums(x))

    W <- (S / alpha) - 1
    W
  }
)

# Wilson =======================================================================
#' @export
#' @rdname index_wilson
#' @aliases index_wilson,matrix-method
setMethod(
  f = "index_wilson",
  signature = signature(x = "matrix"),
  definition = function(x) {
    x <- x > 0 # presence/absence

    m <- nrow(x)
    first_transect <- x[1, ]
    last_transect <- x[m, ]

    # Mean taxa richness
    alpha <- mean(rowSums(x))
    # Number of taxa gained along the transect
    absent <- x[, !first_transect] # Keep only taxa absent in the first transect
    gained <- sum(apply(X = absent, MARGIN = 2, FUN = cumsum)[m, ] > 0)
    # Number of taxa lost
    lost <- sum(!last_transect)

    beta <- (gained + lost) / (2 * alpha)
    beta
  }
)

# PLOT DIVERSITY
#' @include AllGenerics.R
NULL

# DiversityIndex ===============================================================
#' @export
#' @method plot DiversityIndex
plot.DiversityIndex <- function(x, log = "x",
                                col.mean = "#DDAA33", col.interval = "#004488",
                                lty.mean = "solid", lty.interval = "dashed",
                                lwd.mean = 1, lwd.interval = 1,
                                xlab = NULL, ylab = NULL,
                                main = NULL, sub = NULL,
                                ann = graphics::par("ann"),
                                axes = TRUE, frame.plot = axes,
                                panel.first = NULL, panel.last = NULL, ...) {
  ## Prepare data
  count <- as.data.frame(x)

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

  ## Set plotting coordinates
  xlim <- range(count$size)
  ylim <- range(count$index)
  if (length(x@simulation) != 0) {
    xlim <- range(x@simulation[, "size"])
    ylim <- range(x@simulation[, c("lower", "upper")])
  }
  graphics::plot.window(xlim = xlim, ylim = ylim, log = log)

  ## Evaluate pre-plot expressions
  panel.first

  ## Plot
  graphics::points(x = count$size, y = count$index, ...)

  ## Simulated assemblages
  if (length(x@simulation) != 0) {
    refined <- x@simulation
    graphics::lines(x = refined[, "size"], y = refined[, "mean"],
                    col = col.mean, lty = lty.mean, lwd = lwd.mean)
    graphics::lines(x = refined[, "size"], y = refined[, "lower"],
                    col = col.interval, lty = lty.interval, lwd = lwd.interval)
    graphics::lines(x = refined[, "size"], y = refined[, "upper"],
                    col = col.interval, lty = lty.interval, lwd = lwd.interval)
  }

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

  ## Construct axis
  if (axes) {
    graphics::axis(side = 1, las = 1)
    graphics::axis(side = 2, las = 1)
  }

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

  ## Add annotation
  if (ann) {
    y_lab <- switch (
      class(x),
      HeterogeneityIndex = tr_("Heterogeneity"),
      EvennessIndex = tr_("Evenness"),
      RichnessIndex = tr_("Richness"),
      tr_("Diversity")
    )
    xlab <- xlab %||% tr_("Sample size")
    ylab <- ylab %||% sprintf("%s (%s)", y_lab, x@method)
    graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab)
  }

  invisible(x)
}

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

# SHE ==========================================================================
#' @export
#' @rdname she
#' @aliases she,matrix-method
setMethod(
  f = "she",
  signature = c(object = "matrix"),
  definition = function(object, unbiased = FALSE,
                        xlab = NULL, ylab = NULL,
                        main = NULL, sub = NULL,
                        ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = axes,
                        panel.first = NULL, panel.last = NULL,
                        legend = list(x = "right"), ...) {
    ## Prepare data
    data <- .she(object, unbiased = unbiased)

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

    ## Set plotting coordinates
    xlim <- range(data[, 4])
    ylim <- range(data[, -4])
    graphics::plot.window(xlim = xlim, ylim = ylim)

    ## Evaluate pre-plot expressions
    panel.first

    ## Plot
    col <- c("black", "red", "blue")
    for (i in c(1, 2, 3)) {
      graphics::lines(x = data[, 4], y = data[, i], col = col[i], lty = i)
      graphics::points(x = data[, 4], y = data[, i], col = col[i], pch = 16)
    }

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

    ## Construct axis
    if (axes) {
      graphics::axis(side = 1, las = 1)
      graphics::axis(side = 2, las = 1)
    }

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

    ## Add annotation
    if (ann) {
      xlab <- xlab %||% tr_("Sample size")
      ylab <- ylab %||% tr_("Diversity")
      graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
    }

    ## Legend
    if (is.list(legend) && length(legend) > 0) {
      args <- list(legend = c("ln(S)", "H", "ln(E)"),
                   col = col, pch = 16,
                   lty = c(1, 2, 3), bty = "n")
      args <- utils::modifyList(args, legend)
      do.call(graphics::legend, args = args)
    }

    invisible(object)
  }
)

#' @export
#' @rdname she
#' @aliases she,data.frame-method
setMethod(
  f = "she",
  signature = c(object = "data.frame"),
  definition = function(object, unbiased = FALSE,
                        xlab = NULL, ylab = NULL,
                        main = NULL, sub = NULL,
                        ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = axes,
                        panel.first = NULL, panel.last = NULL,
                        legend = list(x = "right"), ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, unbiased = unbiased,
                         xlab = xlab, ylab = ylab,
                         main = main, sub = sub, ann = ann,
                         axes = axes, frame.plot = frame.plot,
                         panel.first = panel.first, panel.last = panel.last,
                         legend = legend, ...)
  }
)

#' @param object A matrix.
#' @noRd
.she <- function(object, unbiased = FALSE, ...) {

  n <- nrow(object)
  m <- ncol(object)

  SHE <- matrix(data = 0, nrow = n, ncol = 4)
  y <- numeric(m)

  for (i in seq_len(n)) {
    x <- object[i, ] + y
    n <- sum(x)

    ## Log species abundance
    S <- log(sum(x > 0)) # Observed species

    ## Shannon index
    H <- index_shannon(x, evenness = FALSE, unbiased = unbiased)

    ## log evenness
    E <- H - S

    SHE[i, ] <- c(S, H, E, n)
    y <- x
  }

  SHE
}

# Profile ======================================================================
#' @export
#' @rdname profiles
#' @aliases profiles,matrix-method
setMethod(
  f = "profiles",
  signature = c(object = "matrix"),
  definition = function(object, alpha = seq(from = 0, to = 4, by = 0.04),
                        color = NULL, symbol = FALSE,
                        xlab = NULL, ylab = NULL,
                        main = NULL, sub = NULL,
                        ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = axes,
                        panel.first = NULL, panel.last = NULL,
                        legend = list(x = "topright"), ...) {
    ## Prepare data
    alpha <- alpha[alpha != 1]
    data <- .profiles(object, alpha = alpha)
    lab <- rownames(object)
    n <- nrow(object)

    ## Graphical parameters
    lwd <- list(...)$lwd %||% graphics::par("lwd")
    if (length(lwd) == 1) lwd <- rep(lwd, length.out = n)

    lty <- list(...)$lty %||% graphics::par("lty")
    if (length(lty) == 1) lty <- rep(lty, length.out = n)
    if (!isFALSE(symbol)) lty <- khroma::palette_line(symbol)(lab)

    col <- list(...)$col %||% graphics::par("col")
    if (length(col) == 1) col <- rep(col, length.out = n)
    if (!isFALSE(color)) col <- khroma::palette_color_discrete(color)(lab)

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

    ## Set plotting coordinates
    xlim <- range(alpha)
    ylim <- range(data, finite = TRUE)
    graphics::plot.window(xlim = xlim, ylim = ylim)

    ## Evaluate pre-plot expressions
    panel.first

    ## Plot
    for (i in seq_len(n)) {
      graphics::lines(x = alpha, y = data[, i], col = col[i],
                      lty = lty[i], lwd = lwd[i])
    }

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

    ## Construct axis
    if (axes) {
      graphics::axis(side = 1, las = 1)
      graphics::axis(side = 2, las = 1)
    }

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

    ## Add annotation
    if (ann) {
      xlab <- xlab %||% "alpha"
      ylab <- ylab %||% tr_("Diversity")
      graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
    }

    ## Legend
    if (is.list(legend) && length(legend) > 0) {
      args <- list(legend = lab, col = col, lty = lty, lwd = lwd, bty = "n")
      args <- utils::modifyList(args, legend)
      do.call(graphics::legend, args = args)
    }

    invisible(object)
  }
)

#' @export
#' @rdname profiles
#' @aliases profiles,data.frame-method
setMethod(
  f = "profiles",
  signature = c(object = "data.frame"),
  definition = function(object, alpha = seq(from = 0, to = 4, by = 0.04),
                        color = NULL, symbol = FALSE,
                        xlab = NULL, ylab = NULL,
                        main = NULL, sub = NULL,
                        ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = axes,
                        panel.first = NULL, panel.last = NULL,
                        legend = list(x = "topright"), ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, alpha = seq(from = 0, to = 4, by = 0.04),
                         color = color, symbol = symbol,
                         xlab = xlab, ylab = ylab,
                         main = main, sub = sub, ann = ann,
                         axes = axes, frame.plot = frame.plot,
                         panel.first = panel.first, panel.last = panel.last,
                         legend = legend, ...)
  }
)

#' @param object A matrix.
#' @noRd
.profiles <- function(object, alpha = seq(from = 0, to = 4, by = 0.04), ...) {

  n <- nrow(object)
  m <- length(alpha)

  prof <- matrix(data = 0, nrow = m, ncol = n)

  index_renyi <- function(z, na.rm = FALSE) {
    z <- z[z > 0] # Remove unobserved species
    if (na.rm) z <- stats::na.omit(z) # Remove NAs

    function(x) {
      p <- z / sum(z)
      exp(log(sum(p^x)) / (1 - x))
    }
  }

  for (i in seq_len(n)) {
    x <- object[i, ]
    f <- index_renyi(x)
    prof[, i] <- vapply(X = alpha, FUN = f, FUN.VALUE = numeric(1), ...)
  }

  prof
}

# SERIOGRAPH
#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname seriograph
#' @aliases eppm,matrix-method
setMethod(
  f = "eppm",
  signature = signature(object = "matrix"),
  definition = function(object) {
    eppm <- (object - expected(object)) * 100 / rowSums(object)
    eppm[eppm < 0] <- 0
    dimnames(eppm) <- dimnames(object)
    eppm
  }
)

#' @export
#' @rdname seriograph
#' @aliases eppm,data.frame-method
setMethod(
  f = "eppm",
  signature = signature(object = "data.frame"),
  definition = function(object) {
    object <- data.matrix(object)
    methods::callGeneric(object)
  }
)

#' @export
#' @rdname seriograph
#' @aliases seriograph,matrix-method
setMethod(
  f = "seriograph",
  signature = signature(object = "matrix"),
  definition = function(object, weights = FALSE,
                        fill = "darkgrey", border = NA,
                        axes = TRUE, ...) {
    plot_ford(object, weights = weights, EPPM = TRUE,
              fill = fill, border = border, axes = axes)
    invisible(object)
  }
)

#' @export
#' @rdname seriograph
#' @aliases seriograph,data.frame-method
setMethod(
  f = "seriograph",
  signature = signature(object = "data.frame"),
  definition = function(object, weights = FALSE,
                        fill = "darkgrey", border = NA,
                        axes = TRUE, ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, weights = weights,
                         fill = fill, border = border, axes = axes)
  }
)

# ALPHA DIVERSITY
#' @include AllGenerics.R
NULL

# Observed =====================================================================
#' @export
#' @rdname observed
#' @aliases observed,numeric-method
setMethod(
  f = "observed",
  signature = c(x = "numeric"),
  definition = function(x, na.rm = FALSE, ...) {
    sum(x > 0, na.rm = na.rm)
  }
)

nobserved <- function(x, n, na.rm = FALSE) {
  sum(x == n, na.rm = na.rm)
}

#' @export
#' @rdname observed
#' @aliases singleton,numeric-method
setMethod(
  f = "singleton",
  signature = c(x = "numeric"),
  definition = function(x, na.rm = FALSE, ...) {
    nobserved(x, n = 1, na.rm = na.rm)
  }
)

#' @export
#' @rdname observed
#' @aliases doubleton,numeric-method
setMethod(
  f = "doubleton",
  signature = c(x = "numeric"),
  definition = function(x, na.rm = FALSE, ...) {
    nobserved(x, n = 2, na.rm = na.rm)
  }
)

# ACE ==========================================================================
#' @export
#' @rdname index_ace
#' @aliases index_ace,numeric-method
setMethod(
  f = "index_ace",
  signature = c(x = "numeric"),
  definition = function(x, k = 10, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    S_rare <- sum(x <= k) # Number of rare species
    S_abund <- sum(x > k) # Number of abundant species
    N_rare <- sum(x[x <= k]) # Number of individuals in the rare species

    F1 <- sum(x == 1) # Number of singleton species
    if (F1 == N_rare) {
      if (getOption("tabula.verbose")) {
        message(
          tr_("ACE is undefined when all rare species are singletons."), "\n",
          tr_("Consider using the bias-corrected Chao1 estimator instead.")
        )
      }
      return(NA_real_)
    }

    ## Sample coverage estimate for rare species
    ## ie. proportion of all individuals in rare species that are not singletons
    C_rare <- 1 - (F1 / N_rare)

    # Coefficient of variation
    a <- sum(vapply(
      X = seq_len(k),
      FUN = function(i, x) { i * (i - 1) * sum(x == i) },
      FUN.VALUE = double(1), x = x)
    )
    g2 <- max((S_rare / C_rare) * (a / (N_rare * (N_rare - 1))) - 1, 0)

    D <- S_abund + (S_rare / C_rare) + (F1 / C_rare) * g2
    D
  }
)

# Berger-Parker ================================================================
#' @export
#' @rdname index_berger
#' @aliases index_berger,numeric-method
setMethod(
  f = "index_berger",
  signature = c(x = "numeric"),
  definition = function(x, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    Nmax <- max(x)
    N <- sum(x)
    d <- Nmax / N
    d
  }
)

# Boone ========================================================================
#' @export
#' @rdname index_boone
#' @aliases index_boone,matrix-method
setMethod(
  f = "index_boone",
  signature = c(x = "matrix"),
  definition = function(x, j = NULL, na.rm = FALSE, ...) {
    ## Validation
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    Y <- colSums(x) # Site-wide totals of each artifact class
    if (is.null(j)) j <- which.max(Y)

    W <- Y[j] / Y # Weighting factor
    P <- W * Y / sum(W * Y)

    W <- matrix(W, nrow = nrow(x), ncol = ncol(x), byrow = TRUE)
    P <- matrix(P, nrow = nrow(x), ncol = ncol(x), byrow = TRUE)

    Wx <- W * x # Weighted counts
    px <- Wx / rowSums(Wx) # Weighted percentages

    rowSums((px - P)^2)
  }
)

# Brillouin ====================================================================
#' @export
#' @rdname index_brillouin
#' @aliases index_brillouin,numeric-method
setMethod(
  f = "index_brillouin",
  signature = c(x = "numeric"),
  definition = function(x, evenness = FALSE, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    N <- sum(x)
    bri <- (lfactorial(N) - sum(lfactorial(x))) / N

    if (evenness) {
      N <- sum(x)
      S <- length(x) # richness = number of different species
      a <- trunc(N / S)
      r <- N - S * a
      c <- (S - r) * lfactorial(a) + r * lfactorial(a + 1)
      HBmax <- (1 / N) * (lfactorial(N) - c)
      bri <- bri / HBmax
    }

    bri
  }
)

# Chao1 ========================================================================
#' @export
#' @rdname index_chao1
#' @aliases index_chao1,numeric-method
setMethod(
  f = "index_chao1",
  signature = c(x = "numeric"),
  definition = function(x, unbiased = FALSE, improved = FALSE, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    S <- length(x) # Number of observed species
    N <- sum(x) # Number of individuals
    f1 <- sum(x == 1) # Number of singleton species
    f2 <- sum(x == 2) # Number of doubleton species

    if (unbiased) {
      D <- S + (((N - 1) / N) * f1 * (f1 - 1)) / (2 * (f2 + 1))
    } else {
      if (f2 == 0) {
        D <- S + ((N - 1) / N) * f1 * ((f1 - 1) / 2)
      } else {
        D <- S + ((N - 1) / N) * (f1^2 / (2 * f2))
      }
    }
    if (improved) {
      f3 <- sum(x == 3) # Number of triple species
      f4 <- sum(x == 4) # Number of quadruple species
      if (f4 == 0) {
        if (getOption("tabula.verbose")) {
          message(
            tr_("Improved Chao1 estimator is undefined when there is no quadruple species.")
          )
        }
        return(NA_real_)
      }

      k <- f1 - ((N - 3) / (N - 1)) * ((f2 * f3) / (2 * f4))
      D <- D + ((N - 3) / N) * (f3 / (4 * f4)) * max(k, 0)
    }

    D
  }
)

# Chao2 ========================================================================
#' @export
#' @rdname index_chao2
#' @aliases index_chao2,matrix-method
setMethod(
  f = "index_chao2",
  signature = c(x = "matrix"),
  definition = function(x, unbiased = FALSE, improved = FALSE, ...) {
    x <- x > 0 # Convert to incidence

    q <- colSums(x) # Number of species in the assemblage
    q <- q[q > 0] # Remove unobserved species
    S <- length(q) # Number of observed species
    t <- nrow(x) # Total number of sampling units
    q1 <- sum(q == 1) # Number of unique species in the assemblage
    q2 <- sum(q == 2) # Number of duplicate species in the assemblage

    if (unbiased) {
      D <- S + ((t - 1) / t) * q1 * ((q1 - 1) / (2 * (q2 + 1)))
    } else {
      if (q2 == 0) {
        D <- S + ((t - 1) / t) * q1 * ((q1 - 1) / 2)
      } else {
        D <- S + ((t - 1) / t) * (q1^2 / (2 * q2))
      }
    }
    if (improved) {
      q3 <- sum(q == 3) # Number of triple species
      q4 <- sum(q == 4) # Number of quadruple species
      if (q4 == 0) {
        if (getOption("tabula.verbose")) {
          message(
            tr_("Improved Chao2 estimator is undefined when there is no quadruple species.")
          )
        }
        return(NA_real_)
      }

      k <- q1 - ((t - 3) / (t - 1)) * ((q2 * q3) / (2 * q4))
      D <- D + ((t - 3) / (4 * t)) * (q3 / q4) * max(k, 0)
    }

    D
  }
)

# ICE ==========================================================================
#' @export
#' @rdname index_ice
#' @aliases index_ice,matrix-method
setMethod(
  f = "index_ice",
  signature = c(x = "matrix"),
  definition = function(x, k = 10, ...) {
    x <- x > 0 # Convert to incidence

    q <- colSums(x) # Number of species in the assemblage
    q <- q[q > 0] # Remove unobserved species

    S_infr <- sum(q <= k) # Number of infrequent species
    S_freq <- sum(q > k) # Number of frequent species
    N_infr <- sum(q[q <= k]) # Number of incidences in the infrequent species
    # Number of sampling units that include at least one infrequent species
    t <- sum(rowSums(x[, q <= k]) != 0)

    q1 <- sum(q == 1) # Number of unique species in the assemblage
    if (q1 == N_infr) {
      if (getOption("tabula.verbose")) {
        message(
          tr_("ICE is undefined when all rare species are singletons."), "\n",
          tr_("Consider using the bias-corrected Chao2 estimator instead.")
        )
      }
      return(NA_real_)
    }

    ## Sample coverage estimate
    ## ie. proportion of all incidences of infrequent species that are not uniques
    C_infr <- 1 - (q1 / N_infr)

    ## Coefficient of variation
    a <- sum(vapply(
      X = seq_len(k),
      FUN = function(x, q) { x * (x - 1) * sum(q == x) },
      FUN.VALUE = double(1), q = q)
    )
    b <- sum(vapply(
      X = seq_len(k),
      FUN = function(x, q) { x * sum(q == x) },
      FUN.VALUE = double(1), q = q)
    )
    c <- sum(vapply(
      X = seq_len(k),
      FUN = function(x, q) { x * sum(q == x) - 1 },
      FUN.VALUE = double(1), q = q)
    )
    g2 <- max((S_infr / C_infr) * (t / (t - 1)) * (a / (b * c)) - 1, 0)

    D <- S_freq + S_infr / C_infr + q1 * g2 / C_infr
    D
  }
)

# Margalef =====================================================================
#' @export
#' @rdname index_margalef
#' @aliases index_margalef,numeric-method
setMethod(
  f = "index_margalef",
  signature = c(x = "numeric"),
  definition = function(x, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    N <- sum(x) # Number of individuals
    S <- length(x) # Number of observed species
    D <- (S - 1) / log(N, base = exp(1))
    D
  }
)

# McIntosh =====================================================================
#' @export
#' @rdname index_mcintosh
#' @aliases index_mcintosh,numeric-method
setMethod(
  f = "index_mcintosh",
  signature = c(x = "numeric"),
  definition = function(x, evenness = FALSE, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    N <- sum(x)
    S <- length(x) # richness = number of different species
    U <- sqrt(sum(x^2))

    if (evenness) {
      mac <- (N - U) / (N - (N / sqrt(S)))
    } else {
      mac <- (N - U) / (N - sqrt(N))
    }
    mac
  }
)

# Menhinick ====================================================================
#' @export
#' @rdname index_menhinick
#' @aliases index_menhinick,numeric-method
setMethod(
  f = "index_menhinick",
  signature = c(x = "numeric"),
  definition = function(x, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    N <- sum(x) # Number of individuals
    S <- length(x) # Number of observed species
    D <- S / sqrt(N)
    D
  }
)

# Shannon ======================================================================
#' @export
#' @rdname index_shannon
#' @aliases index_shannon,numeric-method
setMethod(
  f = "index_shannon",
  signature = c(x = "numeric"),
  definition = function(x, evenness = FALSE, unbiased = FALSE, ACE = FALSE,
                        base = exp(1), na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    N <- sum(x)
    S <- length(x) # richness = number of different species
    p <- x / N
    Hmax <- log(p, base)
    Hmax[is.infinite(Hmax)] <- 0

    H <- -sum(p * Hmax)
    if (unbiased) {
      if (ACE) S <- index_ace(x, ...)
      H <- H + (S - 1) / (2 * N)
    }
    if (evenness) H <- H / log(S)
    H
  }
)

# Simpson ======================================================================
#' @export
#' @rdname index_simpson
#' @aliases index_simpson,numeric-method
setMethod(
  f = "index_simpson",
  signature = c(x = "numeric"),
  definition = function(x, evenness = FALSE, unbiased = FALSE, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    N <- sum(x)
    if (unbiased) D <- sum(x * (x - 1)) / (N * (N - 1))
    else D <- sum((x / N)^2)

    if (evenness) {
      D <- 1 / D
      S <- length(x)
      D <- D / S
    }
    D
  }
)

# Squares ======================================================================
#' @export
#' @rdname index_squares
#' @aliases index_squares,numeric-method
setMethod(
  f = "index_squares",
  signature = c(x = "numeric"),
  definition = function(x, na.rm = FALSE, ...) {
    ## Validation
    x <- x[x > 0] # Remove unobserved species
    if (na.rm) x <- stats::na.omit(x) # Remove NAs
    if (anyNA(x)) return(NA)

    S <- length(x) # Number of observed species
    N <- sum(x) # Number of individuals
    f1 <- sum(x == 1) # Number of singleton species

    Ssq <- S + f1^2 * sum(x^2) / (N^2 - f1 * S)
    Ssq
  }
)

# RAREFACTION
#' @include AllGenerics.R
NULL

#' @export
#' @rdname rarefaction
#' @aliases rarefaction,matrix-method
setMethod(
  f = "rarefaction",
  signature = signature(object = "matrix"),
  definition = function(object, sample = NULL, method = c("hurlbert", "baxter"),
                        step = 1) {
    ## Validation
    method <- match.arg(method, several.ok = FALSE)

    n <- nrow(object)
    if (is.null(sample)) {
      sample <- rowSums(object)
    }
    if (length(sample) == 1) {
      sample <- rep(sample, n)
    }
    k <- seq(from = 1, to = max(sample), by = step)

    ## Matrix of results
    z <- matrix(data = NA_real_, nrow = n, ncol = length(k))
    row_names <- rownames(object) %||% paste0("S", seq_len(n)) # Fix names
    dimnames(z) <- list(row_names, k)

    for (i in seq_len(n)) {
      spl <- k[k <= sample[[i]]]
      rare <- vapply(
        X = spl,
        FUN = function(s, x, f) f(x, s),
        FUN.VALUE = numeric(1),
        x = object[i, ],
        f = get_index(method) # Select method
      )
      z[i, seq_along(rare)] <- rare
    }

    .RarefactionIndex(
      z,
      labels = row_names,
      size = as.integer(k),
      method = method
    )
  }
)

#' @export
#' @rdname rarefaction
#' @aliases rarefaction,data.frame-method
setMethod(
  f = "rarefaction",
  signature = signature(object = "data.frame"),
  definition = function(object, sample = NULL, method = c("hurlbert", "baxter"),
                        step = 1) {
    object <- data.matrix(object)
    methods::callGeneric(object, sample = sample, method = method, step = step)
  }
)

# Plot =========================================================================
#' @export
#' @method plot RarefactionIndex
plot.RarefactionIndex <- function(x, color = NULL, symbol = FALSE,
                                  xlab = NULL, ylab = NULL,
                                  main = NULL, sub = NULL,
                                  ann = graphics::par("ann"),
                                  axes = TRUE, frame.plot = axes,
                                  panel.first = NULL, panel.last = NULL,
                                  legend = list(x = "topleft"), ...) {
  ## Prepare data
  n <- nrow(x)
  lab <- labels(x)

  ## Graphical parameters
  lwd <- list(...)$lwd %||% graphics::par("lwd")
  if (length(lwd) == 1) lwd <- rep(lwd, length.out = n)

  lty <- list(...)$lty %||% graphics::par("lty")
  if (length(lty) == 1) lty <- rep(lty, length.out = n)
  if (!isFALSE(symbol)) lty <- khroma::palette_line(symbol)(lab)

  col <- list(...)$col %||% graphics::par("col")
  if (length(col) == 1) col <- rep(col, length.out = n)
  if (!isFALSE(color)) col <- khroma::palette_color_discrete(color)(lab)

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

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

  ## Evaluate pre-plot expressions
  panel.first

  ## Plot
  for (i in seq_len(n)) {
    graphics::lines(x = x@size, y = x[i, ], col = col[i],
                    lwd = lwd[i], lty = lty[i])
  }

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

  ## Construct axis
  if (axes) {
    graphics::axis(side = 1, las = 1)
    graphics::axis(side = 2, las = 1)
  }

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

  ## Add annotation
  if (ann) {
    xlab <- xlab %||% tr_("Sample size")
    ylab <- ylab %||% tr_("Expected species index")
    graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab)
  }

  ## Legend
  if (is.list(legend) && length(legend) > 0) {
    args <- list(legend = lab, col = col, lty = lty, lwd = lwd, bty = "n")
    args <- utils::modifyList(args, legend)
    do.call(graphics::legend, args = args)
  }

  invisible(x)
}

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

# Index ========================================================================
#' @export
#' @rdname index_baxter
#' @aliases index_baxter,numeric-method
setMethod(
  f = "index_baxter",
  signature = signature(x = "numeric"),
  definition = function(x, sample, ...) {
    x <- x[x > 0]
    N <- sum(x)

    E <- suppressWarnings(exp(ramanujan(N - x) + ramanujan(N - sample) -
                                ramanujan(N - x - sample) - ramanujan(N)))
    sum(1 - E, na.rm = FALSE)
  }
)

#' @export
#' @rdname index_hurlbert
#' @aliases index_hurlbert,numeric-method
setMethod(
  f = "index_hurlbert",
  signature = signature(x = "numeric"),
  definition = function(x, sample, ...) {
    x <- x[x > 0]
    N <- sum(x)

    E <- vapply(
      X = x,
      FUN = function(x, N, sample) {
        if (N - x > sample) {
          combination(N - x, sample) / combination(N, sample)
        } else {
          0
        }
      },
      FUN.VALUE = double(1),
      N, sample
    )
    sum(1 - E, na.rm = FALSE)
  }
)

# MUTATORS
#' @include AllClasses.R
NULL

# Getters ======================================================================
#' @export
#' @method labels DiversityIndex
labels.DiversityIndex <- function(object, ...) object@labels

#' @rdname mutators
#' @aliases labels,DiversityIndex-method
setMethod("labels", "DiversityIndex", labels.DiversityIndex)

#' @export
#' @method labels RarefactionIndex
labels.RarefactionIndex <- function(object, ...) object@labels

#' @rdname mutators
#' @aliases labels,RarefactionIndex-method
setMethod("labels", "RarefactionIndex", labels.RarefactionIndex)

#' @export
#' @rdname mutators
#' @aliases get_method,DiversityIndex-method
setMethod(
  f = "get_method",
  signature = "DiversityIndex",
  definition = function(x) x@method
)

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

# HELPERS

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

is_incidence <- function(x) {
  if (is.logical(x)) return(TRUE)
  x <- as.numeric(x)
  all(x == 0 | x == 1, na.rm = TRUE)
}

#' Rolling Sum
#'
#' @param x A [`numeric`] vector.
#' @param n An [`integer`] giving the rolling window size.
#' @return A [`numeric`] vector.
#' @keywords internal
#' @noRd
roll_sum <- function(x, n = 2) {
  utils::tail(cumsum(x) - cumsum(c(rep(0, n), utils::head(x, -n))), -n + 1)
}

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

#' Circle
#'
#' Draws a circle.
#' @param x,y A length-one [`numeric`] vector giving the coordinates of the
#'  center of the circle.
#' @param radius A length-one [`numeric`] vector giving the radius of the
#'  circle.
#' @param n A length-on [`integer`] vector specifying the number of vertices to
#'  draw the circle.
#' @param ... Further parameters to be passed to [graphics::polygon()].
#' @return
#'  `circle()` is called it for its side-effects: it results in a graphic
#'  being displayed.
#' @example inst/examples/ex-circle.R
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
circle <- function(x, y, radius, ..., n = 100) {
  angle.inc <- 2 * pi / n
  angles <- seq(0, 2 * pi - angle.inc, by = angle.inc)

  xv <- cos(angles) * radius + x
  yv <- sin(angles) * radius + y
  graphics::polygon(xv, yv, ...)
}

#' Label Percentages
#'
#' @param x A [`numeric`] vector.
#' @param digits An [`integer`] indicating the number of decimal places.
#'  If `NULL` (the default), breaks will have the minimum number of digits
#'  needed to show the difference between adjacent values.
#' @param trim A [`logical`] scalar. If `FALSE` (the default), values are
#'  right-justified to a common width (see [base::format()]).
#' @return A [`character`] vector.
#' @keywords internal
#' @noRd
label_percent <- function(x, digits = NULL, trim = FALSE) {
  i <- !is.na(x)
  y <- x[i]
  y <- abs(y) * 100
  y <- format(y, trim = trim, digits = digits)
  y <- paste0(y, "%")
  x[i] <- y
  x
}

# GENERIC METHODS
#' @include AllClasses.R
NULL

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

# Extract ======================================================================
## Mutators --------------------------------------------------------------------
#' Get or Set Parts of an Object
#'
#' Getters and setters to extract or replace parts of an object.
#' @param object,x An \R object from which to get or set element(s).
#' @param ... Currently not used.
# @param value A possible value for the element(s) of `object` (see below).
#' @return
#'  * `labels()` returns a suitable set of labels from an object for use in
#'    printing or plotting.
# @example inst/examples/ex-mutator.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name mutators
#' @rdname mutators
#' @aliases get set
NULL

#' @rdname mutators
#' @aliases get_method-method
setGeneric(
  name = "get_method",
  def = function(x) standardGeneric("get_method")
)

## Coerce ----------------------------------------------------------------------
#' Coerce to a Data Frame
#'
#' @param x An object.
#' @param row.names,optional Currently not used.
#' @param ... Currently not used.
#' @return
#'  A [`data.frame`].
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name data.frame
#' @rdname data.frame
NULL

# Statistic ====================================================================
#' Bootstrap Estimation
#'
#' Samples randomly from the elements of `object` with replacement.
#' @param object An \R object (typically a [DiversityIndex-class] object).
#' @param n A non-negative [`integer`] giving the number of bootstrap
#'  replications.
#' @param f A [`function`] that takes a single numeric vector (the result of
#'  `do`) as argument.
#' @return
#'  If `f` is `NULL` (the default), `bootstrap()` returns a named `numeric`
#'  vector with the following elements:
#'  \describe{
#'   \item{`original`}{The observed value of `do` applied to `object`.}
#'   \item{`mean`}{The bootstrap estimate of mean of `do`.}
#'   \item{`bias`}{The bootstrap estimate of bias of `do`.}
#'   \item{`error`}{he bootstrap estimate of standard error of `do`.}
#'  }
#'
#'  If `f` is a `function`, `bootstrap()` returns the result of `f` applied to
#'  the `n` values of `do`.
#' @example inst/examples/ex-bootstrap.R
#' @author N. Frerebeau
#' @docType methods
#' @family resampling methods
#' @name bootstrap
#' @rdname bootstrap
NULL

#' Jackknife Estimation
#'
#' @param object An \R object (typically a [DiversityIndex-class] object).
#' @param f A [`function`] that takes a single numeric vector (the leave-one-out
#'  values of `do`) as argument.
#' @return
#'  If `f` is `NULL` (the default), `jackknife()` returns a named `numeric`
#'  vector with the following elements:
#'  \describe{
#'   \item{`original`}{The observed value of `do` applied to `object`.}
#'   \item{`mean`}{The jackknife estimate of mean of `do`.}
#'   \item{`bias`}{The jackknife estimate of bias of `do`.}
#'   \item{`error`}{he jackknife estimate of standard error of `do`.}
#'  }
#'
#'  If `f` is a `function`, `jackknife()` returns the result of `f` applied to
#'  the leave-one-out values of `do`.
#' @example inst/examples/ex-jackknife.R
#' @author N. Frerebeau
#' @docType methods
#' @family resampling methods
#' @name jackknife
#' @rdname jackknife
NULL

#' Resample
#'
#' Simulates observations from a multinomial distribution.
#' @param object A [`numeric`] vector of count data (absolute frequencies).
#' @param do A [`function`] that takes `object` 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()]
#' @example inst/examples/ex-resample.R
#' @author N. Frerebeau
#' @docType methods
#' @family resampling methods
#' @aliases resample-method
setGeneric(
  name = "resample",
  def = function(object, ...) standardGeneric("resample")
)

# Diversity ====================================================================
## Heterogeneity ---------------------------------------------------------------
#' Heterogeneity and Evenness
#'
#' @description
#'  * `heterogeneity()` computes an heterogeneity or dominance index.
#'  * `evenness()` computes an evenness measure.
#' @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 method A [`character`] string specifying the index to be computed
#'  (see details). Any unambiguous substring can be given.
#' @param evenness A [`logical`] scalar: should an evenness measure be computed
#'  instead of an heterogeneity/dominance index?
#' @param ... Further arguments to be passed to internal methods (see below).
#' @details
#'  *Diversity* measurement assumes that all individuals in a specific
#'  taxa are equivalent and that all types are equally different from each
#'  other (Peet 1974). A measure of diversity can be achieved by using indices
#'  built on the relative abundance of taxa. These indices (sometimes referred
#'  to as non-parametric indices) benefit from not making assumptions about the
#'  underlying distribution of taxa abundance: they only take relative
#'  abundances of the species that are present and species richness into
#'  account. Peet (1974) refers to them as indices of *heterogeneity*.
#'
#'  Diversity indices focus on one aspect of the taxa abundance and emphasize
#'  either *[richness][richness()]* (weighting towards uncommon taxa) or
#'  *dominance* (weighting towards abundant taxa; Magurran 1988).
#'
#'  *Evenness* is a measure of how evenly individuals are distributed across the
#'  sample.
#' @section Heterogeneity and Evenness Measures:
#'  The following heterogeneity index and corresponding evenness measures
#'  are available (see Magurran 1988 for details):
#'  \describe{
#'   \item{`berger`}{[Berger-Parker dominance index][index_berger()].}
#'   \item{`boone`}{[Boone heterogeneity measure][index_boone()].}
#'   \item{`brillouin`}{[Brillouin diversity index][index_brillouin()].}
#'   \item{`mcintosh`}{[McIntosh dominance index][index_mcintosh()].}
#'   \item{`shannon`}{[Shannon-Wiener diversity index][index_shannon()].}
#'   \item{`simpson`}{[Simpson dominance index][index_simpson()].}
#'  }
#'
#'  The `berger`, `mcintosh` and `simpson` methods return a *dominance* index,
#'  not the reciprocal or inverse form usually adopted, so that an increase in
#'  the value of the index accompanies a decrease in diversity.
#' @return
#'  * `heterogeneity()` returns an [HeterogeneityIndex-class] object.
#'  * `evenness()` returns an [EvennessIndex-class] object.
#' @seealso [index_berger()], [index_boone()], [index_brillouin()],
#'  [index_mcintosh()], [index_shannon()], [index_simpson()]
#' @references
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press.
#'  \doi{10.1007/978-94-015-7358-0}.
#'
#'  Peet, R. K. (1974). The Measurement of Species Diversity. *Annual Review of
#'  Ecology and Systematics*, 5(1), 285-307.
#'  \doi{10.1146/annurev.es.05.110174.001441}.
#' @example inst/examples/ex-diversity.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @aliases heterogeneity-method
setGeneric(
  name = "heterogeneity",
  def = function(object, ...) standardGeneric("heterogeneity"),
  valueClass = "HeterogeneityIndex"
)

#' @rdname heterogeneity
#' @aliases evenness-method
setGeneric(
  name = "evenness",
  def = function(object, ...) standardGeneric("evenness"),
  valueClass = "EvennessIndex"
)

#' Berger-Parker Dominance Index
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @details
#'  The Berger-Parker index expresses the proportional importance of the most
#'  abundant type. This metric is highly biased by sample size and richness,
#'  moreover it does not make use of all the information available from sample.
#'
#'  This is a *dominance* index, so that an increase in the value of the index
#'  accompanies a decrease in diversity.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Berger, W. H. & Parker, F. L. (1970). Diversity of Planktonic Foraminifera
#'  in Deep-Sea Sediments. *Science*, 168(3937), 1345-1347.
#'  \doi{10.1126/science.168.3937.1345}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_berger-method
setGeneric(
  name = "index_berger",
  def = function(x, ...) standardGeneric("index_berger")
)

#' Boone Heterogeneity Measure
#'
#' @param x A \eqn{m \times p}{m x p} `numeric` [`matrix`] of count data
#'  (absolute frequencies, i.e. a contingency table).
#' @param j An [`integer`] giving the index of the reference type/taxa.
#'  If `NULL` (the default), the most frequent type/taxa in any assemblage will
#'  be used.
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Boone, J. L. (1987). Defining and Measuring Midden Catchment. *American
#'  Antiquity*, 52(2), 336-45. \doi{10.2307/281785}.
#'
#'  Kintigh, K. W. (1989). Sample Size, Significance, and Measures of
#'  Diversity. In Leonard, R. D. and Jones, G. T., *Quantifying Diversity
#'  in Archaeology*. New Directions in Archaeology. Cambridge:
#'  Cambridge University Press, p. 25-36.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_boone-method
setGeneric(
  name = "index_boone",
  def = function(x, ...) standardGeneric("index_boone")
)

#' Brillouin Diversity Index.
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param evenness A [`numeric`] scalar: should evenness be computed?
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @details
#'  The Brillouin index describes a known collection: it does not assume random
#'  sampling in an infinite population. Pielou (1975) and Laxton (1978) argues
#'  for the use of the Brillouin index in all circumstances, especially in
#'  preference to the Shannon index.
#' @note
#'  Ramanujan approximation is used for \eqn{x!} computation if \eqn{x > 170}.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Brillouin, L. (1956). *Science and information theory*. New York:
#'  Academic Press.
#'
#'  Laxton, R. R. (1978). The measure of diversity. *Journal of Theoretical
#'  Biology*, 70(1), 51-67.
#'  \doi{10.1016/0022-5193(78)90302-8}.
#'
#'  Pielou, E. C. (1975). *Ecological Diversity*. New York: Wiley.
#'  \doi{10.4319/lo.1977.22.1.0174b}
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_brillouin-method
setGeneric(
  name = "index_brillouin",
  def = function(x, ...) standardGeneric("index_brillouin")
)

#' McIntosh Dominance Index.
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param evenness A [`numeric`] scalar: should evenness be computed?
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @details
#'  The McIntosh index expresses the heterogeneity of a sample in geometric
#'  terms. It describes the sample as a point of a \eqn{S}-dimensional
#'  hypervolume and uses the Euclidean distance of this point from the origin.
#'
#'  This is a *dominance* index, so that an increase in the value of the index
#'  accompanies a decrease in diversity.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  McIntosh, R. P. (1967). An Index of Diversity and the Relation of Certain
#'  Concepts to Diversity. *Ecology*, 48(3), 392-404.
#'  \doi{10.2307/1932674}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_mcintosh-method
setGeneric(
  name = "index_mcintosh",
  def = function(x, ...) standardGeneric("index_mcintosh")
)

#' Shannon-Wiener Diversity Index
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param evenness A [`numeric`] scalar: should evenness be computed?
#' @param unbiased A [`logical`] scalar: should the bias-corrected estimator be
#'  used?
#' @param ACE A [`logical`] scalar: should the ACE species richness estimator be
#'  used in the bias correction?
#' @param base A positive [`numeric`] value specifying the base with respect to
#'  which logarithms are computed.
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @details
#'  The Shannon index assumes that individuals are randomly sampled from an
#'  infinite population and that all taxa are represented in the sample (it
#'  does not reflect the sample size). The main source of error arises from the
#'  failure to include all taxa in the sample: this error increases as the
#'  proportion of species discovered in the sample declines (Peet 1974,
#'  Magurran 1988). The maximum likelihood estimator (MLE) is used for the
#'  relative abundance, this is known to be negatively biased by sample size.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Peet, R. K. (1974). The Measurement of Species Diversity. *Annual Review of
#'  Ecology and Systematics*, 5(1), 285-307.
#'  \doi{10.1146/annurev.es.05.110174.001441}.
#'
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press.
#'  \doi{10.1007/978-94-015-7358-0}.
#'
#'  Shannon, C. E. (1948). A Mathematical Theory of Communication. *The
#'  Bell System Technical Journal*, 27, 379-423.
#'  \doi{10.1002/j.1538-7305.1948.tb01338.x}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_shannon-method
setGeneric(
  name = "index_shannon",
  def = function(x, ...) standardGeneric("index_shannon")
)

#' Simpson Dominance Index
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param evenness A [`numeric`] scalar: should evenness be computed?
#' @param unbiased A [`logical`] scalar: should the bias-corrected estimator be
#'  used?
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @details
#'  The Simpson index expresses the probability that two individuals randomly
#'  picked from a finite sample belong to two different types. It can be
#'  interpreted as the weighted mean of the proportional abundances. This
#'  metric is a true probability value, it ranges from \eqn{0} (all taxa are
#'  equally present) to \eqn{1} (one taxon dominates the community completely).
#'
#'  This is a *dominance* index, so that an increase in the value of the index
#'  accompanies a decrease in diversity.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Simpson, E. H. (1949). Measurement of Diversity. *Nature*, 163(4148),
#'  688-688. \doi{10.1038/163688a0}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_simpson-method
setGeneric(
  name = "index_simpson",
  def = function(x, ...) standardGeneric("index_simpson")
)

## Richness --------------------------------------------------------------------
#' Richness
#'
#' @description
#'  * `richness()` computes sample richness.
#'  * `composition()` computes asymptotic species richness.
#' @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 method A [`character`] string or vector of strings specifying the
#' index to be computed (see details). Any unambiguous substring can be given.
#' @param ... Further arguments to be passed to internal methods (see below).
#' @section Details:
#'  The number of observed taxa, provides an instantly comprehensible
#'  expression of diversity. While the number of taxa within a sample
#'  is easy to ascertain, as a term, it makes little sense: some taxa
#'  may not have been seen, or there may not be a fixed number of taxa
#'  (e.g. in an open system; Peet 1974). As an alternative, *richness*
#'  (\eqn{S}) can be used for the concept of taxa number (McIntosh 1967).
#'
#'  It is not always possible to ensure that all sample sizes are equal
#'  and the number of different taxa increases with sample size and
#'  sampling effort (Magurran 1988). Then, *[rarefaction][rarefaction()]*
#'  (\eqn{E(S)}) is the number of taxa expected if all samples were of a
#'  standard size (i.e. taxa per fixed number of individuals).
#'  Rarefaction assumes that imbalances between taxa are due to sampling and
#'  not to differences in actual abundances.
#' @section Richness Measures:
#'  The following richness measures are available for count data:
#'  \describe{
#'   \item{`observed`}{Number of observed taxa/types.}
#'   \item{`margalef`}{[Margalef richness index][index_margalef()].}
#'   \item{`menhinick`}{[Menhinick richness index][index_menhinick()].}
#'  }
#'
#' @section Asymptotic Species Richness:
#'  The following measures are available for count data:
#'  \describe{
#'   \item{`ace`}{[Abundance-based Coverage Estimator][index_ace()].}
#'   \item{`chao1`}{(improved/unbiased) [Chao1 estimator][index_chao1()].}
#'   \item{`squares`}{[Squares estimator][index_squares()].}
#'  }
#'
#'  The following measures are available for replicated incidence data:
#'  \describe{
#'   \item{`ice`}{[Incidence-based Coverage Estimator][index_ice()].}
#'   \item{`chao2`}{(improved/unbiased) [Chao2 estimator][index_chao2()].}
#'  }
#' @return
#'  * `richness()` returns a [RichnessIndex-class] object.
#'  * `composition()` returns a [CompositionIndex-class] object.
#' @seealso [index_margalef()], [index_menhinick()], [index_ace()],
#'  [index_chao1()], [index_squares()], [index_ice()], [index_chao2()]
#' @references
#'  Kintigh, K. W. (1989). Sample Size, Significance, and Measures of
#'  Diversity. In Leonard, R. D. and Jones, G. T., *Quantifying Diversity
#'  in Archaeology*. New Directions in Archaeology. Cambridge:
#'  Cambridge University Press, p. 25-36.
#'
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press. \doi{10.1007/978-94-015-7358-0}.
#'
#'  Magurran, A E. & Brian J. McGill (2011). *Biological Diversity:
#'  Frontiers in Measurement and Assessment*. Oxford: Oxford University Press.
#'
#'  McIntosh, R. P. (1967). An Index of Diversity and the Relation of Certain
#'  Concepts to Diversity. *Ecology*, 48(3), 392-404. \doi{10.2307/1932674}.
#'
#'  Peet, R. K. (1974). The Measurement of Species Diversity. *Annual Review of
#'  Ecology and Systematics*, 5(1), 285-307.
#'  \doi{10.1146/annurev.es.05.110174.001441}.
#' @seealso [`plot()`][plot_diversity]
#' @example inst/examples/ex-richness.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @aliases richness-method
setGeneric(
  name = "richness",
  def = function(object, ...) standardGeneric("richness"),
  valueClass = "RichnessIndex"
)

#' @rdname richness
#' @aliases composition-method
setGeneric(
  name = "composition",
  def = function(object, ...) standardGeneric("composition"),
  valueClass = "CompositionIndex"
)

#' Number of Observed Species
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @family alpha diversity measures
#' @docType methods
#' @aliases observed-method
setGeneric(
  name = "observed",
  def = function(x, ...) standardGeneric("observed")
)

#' @rdname observed
#' @aliases singleton-method
setGeneric(
  name = "singleton",
  def = function(x, ...) standardGeneric("singleton")
)

#' @rdname observed
#' @aliases doubleton-method
setGeneric(
  name = "doubleton",
  def = function(x, ...) standardGeneric("doubleton")
)

#' Abundance-based Coverage Estimator
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param k A length-one [`numeric`] vector giving the threshold between
#'  rare/infrequent and abundant/frequent species.
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Chao, A. & Lee, S.-M. (1992). Estimating the Number of Classes via Sample
#'  Coverage. *Journal of the American Statistical Association*, 87(417),
#'  210-217. \doi{10.1080/01621459.1992.10475194}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_ace-method
setGeneric(
  name = "index_ace",
  def = function(x, ...) standardGeneric("index_ace")
)

#' Chao1 Estimator
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param unbiased A [`logical`] scalar: should the bias-corrected estimator be
#'  used?
#' @param improved A [`logical`] scalar: should the improved estimator be used?
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Chao, A. (1984). Nonparametric Estimation of the Number of Classes in a
#'  Population. *Scandinavian Journal of Statistics*, 11(4), 265-270.
#'
#'  Chiu, C.-H., Wang, Y.-T., Walther, B. A. & Chao, A. (2014). An improved
#'  nonparametric lower bound of species richness via a modified good-turing
#'  frequency formula. *Biometrics*, 70(3), 671-682.
#'  \doi{10.1111/biom.12200}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_chao1-method
setGeneric(
  name = "index_chao1",
  def = function(x, ...) standardGeneric("index_chao1")
)

#' Chao2 Estimator
#'
#' @param x A \eqn{m \times p}{m x p} [`matrix`] of presence/absence data
#'  (incidence).
#' @param unbiased A [`logical`] scalar: should the bias-corrected estimator be
#'  used?
#' @param improved A [`logical`] scalar: should the improved estimator be used?
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Chao, A. (1987). Estimating the Population Size for Capture-Recapture Data
#'  with Unequal Catchability. *Biometrics* 43(4), 783-791.
#'
#'  Chiu, C.-H., Wang, Y.-T., Walther, B. A. & Chao, A. (2014). An improved
#'  nonparametric lower bound of species richness via a modified good-turing
#'  frequency formula. *Biometrics*, 70(3), 671-682.
#'  \doi{10.2307/2531532}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_chao2-method
setGeneric(
  name = "index_chao2",
  def = function(x, ...) standardGeneric("index_chao2")
)

#' Incidence-based Coverage Estimator
#'
#' @param x A \eqn{m \times p}{m x p} [`matrix`] of presence/absence data
#'  (incidence).
#' @param k A length-one [`numeric`] vector giving the threshold between
#'  rare/infrequent and abundant/frequent species.
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Chao, A. & Chiu, C.-H. (2016). Species Richness: Estimation and Comparison.
#'  *In* Balakrishnan, N., Colton, T., Everitt, B., Piegorsch, B., Ruggeri,
#'  F. & Teugels, J. L. (Eds.), *Wiley StatsRef: Statistics Reference Online*.
#'  Chichester, UK: John Wiley & Sons, Ltd., 1-26.
#'  \doi{10.1002/9781118445112.stat03432.pub2}
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_ice-method
setGeneric(
  name = "index_ice",
  def = function(x, ...) standardGeneric("index_ice")
)

#' Margalef Richness Index
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Margalef, R. (1958). Information Theory in Ecology. *General Systems*,
#'  3, 36-71.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_margalef-method
setGeneric(
  name = "index_margalef",
  def = function(x, ...) standardGeneric("index_margalef")
)

#' Menhinick Richness Index
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Menhinick, E. F. (1964). A Comparison of Some Species-Individuals Diversity
#'  Indices Applied to Samples of Field Insects. *Ecology*, 45(4), 859-861.
#'  \doi{10.2307/1934933}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_menhinick-method
setGeneric(
  name = "index_menhinick",
  def = function(x, ...) standardGeneric("index_menhinick")
)

#' Squares Estimator
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param na.rm A [`numeric`] scalar: should missing values (including `NaN`) be
#'  removed?
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Alroy, J. (2018). Limits to Species Richness in Terrestrial Communities.
#'  *Ecology Letters*, 21(12), 1781-1789. \doi{10.1111/ele.13152}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_squares-method
setGeneric(
  name = "index_squares",
  def = function(x, ...) standardGeneric("index_squares")
)

## Rarefaction -----------------------------------------------------------------
#' Rarefaction
#'
#' @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 sample A length-one [`numeric`] vector giving the sub-sample size.
#'  The size of sample should be smaller than total community size.
#' @param method A [`character`] string or vector of strings specifying the
#' index to be computed (see details). Any unambiguous substring can be given.
#' @param step An [`integer`] giving the increment of the sample size.
#' @param ... Currently not used.
#' @inheritSection richness Details
#' @section Rarefaction Measures:
#'  The following rarefaction measures are available for count data:
#'  \describe{
#'   \item{`baxter`}{[Baxter's rarefaction][index_baxter()].}
#'   \item{`hurlbert`}{[Hurlbert's unbiased estimate][index_hurlbert()] of
#'   Sander's rarefaction.}
#'  }
#' @return
#'  A [RarefactionIndex-class] object.
#' @example inst/examples/ex-rarefaction.R
#' @seealso [index_baxter()], [index_hurlbert()], [`plot()`][plot_rarefaction]
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @aliases rarefaction-method
setGeneric(
  name = "rarefaction",
  def = function(object, ...) standardGeneric("rarefaction")
)

#' Baxter's Rarefaction
#'
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param sample A length-one [`numeric`] vector giving the sub-sample size.
#'  The size of sample should be smaller than total community size.
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Baxter, M. J. (2001). Methodological Issues in the Study of Assemblage
#'  Diversity. *American Antiquity*, 66(4), 715-725. \doi{10.2307/2694184}.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_baxter-method
setGeneric(
  name = "index_baxter",
  def = function(x, ...) standardGeneric("index_baxter")
)

#' Hurlbert's Rarefaction
#'
#' Hurlbert's unbiased estimate of Sander's rarefaction.
#' @param x A [`numeric`] vector of count data (absolute frequencies).
#' @param sample A length-one [`numeric`] vector giving the sub-sample size.
#'  The size of sample should be smaller than total community size.
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Hurlbert, S. H. (1971). The Nonconcept of Species Diversity: A Critique and
#'  Alternative Parameters. *Ecology*, 52(4), 577-586.
#'  \doi{10.2307/1934145}.
#'
#'  Sander, H. L. (1968). Marine Benthic Diversity: A Comparative Study.
#'  *The American Naturalist*, 102(925), 243-282.
#' @author N. Frerebeau
#' @family alpha diversity measures
#' @docType methods
#' @aliases index_hurlbert-method
setGeneric(
  name = "index_hurlbert",
  def = function(x, ...) standardGeneric("index_hurlbert")
)

## Similarity ------------------------------------------------------------------
#' Similarity
#'
#' @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 method A [`character`] string specifying the method to be
#'  used (see details). Any unambiguous substring can be given.
#' @param ... Currently not used.
#' @details
#'  \eqn{\beta}-diversity can be measured by addressing *similarity*
#'  between pairs of samples/cases.
#'
#'  `bray`, `jaccard`, `morisita` and `sorenson` indices provide a scale of
#'  similarity from \eqn{0}-\eqn{1} where \eqn{1} is perfect similarity and
#'  \eqn{0} is no similarity.
#'  `brainerd` is scaled between \eqn{0} and \eqn{200}.
#'
#'  \describe{
#'   \item{`brainerd`}{[Brainerd-Robinson quantitative index][index_brainerd()].}
#'   \item{`bray`}{[Bray-Curtis similarity (Sorenson quantitative index)][index_bray()].}
#'   \item{`jaccard`}{[Jaccard qualitative index][index_jaccard()].}
#'   \item{`morisita`}{[Morisita-Horn quantitative index][index_morisita()].}
#'   \item{`sorenson`}{[Dice-Sorenson index (Sorenson qualitative index)][index_sorenson()].}
#'  }
#'
#'  For `jaccard` and `sorenson`, data are standardized on a presence/absence
#'  scale (\eqn{0}/\eqn{1}) beforehand.
#' @return
#'  A [stats::dist] object.
#' @seealso [index_binomial()], [index_brainerd()], [index_bray()],
#'  [index_jaccard()], [index_morisita()], [index_sorenson()]
#' @references
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press. \doi{10.1007/978-94-015-7358-0}.
#' @example inst/examples/ex-similarity.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @aliases similarity-method
setGeneric(
  name = "similarity",
  def = function(object, ...) standardGeneric("similarity")
)

#' Jaccard Index
#'
#' @param x,y A [`numeric`] vector.
#' @param ... Currently not used.
#' @details
#'  Data are standardized on a presence/absence scale (\eqn{0}/\eqn{1})
#'  beforehand.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press. \doi{10.1007/978-94-015-7358-0}.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_jaccard-method
setGeneric(
  name = "index_jaccard",
  def = function(x, y, ...) standardGeneric("index_jaccard")
)

#' Dice-Sorenson Index
#'
#' @param x,y A [`numeric`] vector.
#' @param ... Currently not used.
#' @details
#'  Data are standardized on a presence/absence scale (\eqn{0}/\eqn{1})
#'  beforehand.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press. \doi{10.1007/978-94-015-7358-0}.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_sorenson-method
setGeneric(
  name = "index_sorenson",
  def = function(x, y, ...) standardGeneric("index_sorenson")
)

#' Bray-Curtis Similarity
#'
#' Bray and Curtis modified version of the Sorenson index.
#' @param x,y A [`numeric`] vector.
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Bray, J. R. & Curtis, J. T. (1957). An Ordination of the Upland Forest
#'  Communities of Southern Wisconsin. *Ecological Monographs*, 27(4),
#'  325-349. \doi{10.2307/1942268}.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_bray-method
setGeneric(
  name = "index_bray",
  def = function(x, y, ...) standardGeneric("index_bray")
)

#' Morisita-Horn Quantitative Index
#'
#' @param x,y A [`numeric`] vector.
#' @param ... Currently not used.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press. \doi{10.1007/978-94-015-7358-0}.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_morisita-method
setGeneric(
  name = "index_morisita",
  def = function(x, y, ...) standardGeneric("index_morisita")
)

#' Brainerd-Robinson Quantitative Index
#'
#' @param x,y A [`numeric`] vector.
#' @param ... Currently not used.
#' @details
#'  A city-block metric of similarity between pairs of samples/cases.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Brainerd, G. W. (1951). The Place of Chronological Ordering in
#'  Archaeological Analysis. *American Antiquity*, 16(04), 301-313.
#'  \doi{10.2307/276979}.
#'
#'  Robinson, W. S. (1951). A Method for Chronologically Ordering Archaeological
#'  Deposits. *American Antiquity*, 16(04), 293-301. \doi{10.2307/276978}.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_brainerd-method
setGeneric(
  name = "index_brainerd",
  def = function(x, y, ...) standardGeneric("index_brainerd")
)

## Turnover --------------------------------------------------------------------
#' Turnover
#'
#' Returns the degree of turnover in taxa composition along a gradient or
#' transect.
#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
#'  [`data.frame`] of count data or incidence data. A [`data.frame`]
#'  will be coerced to a `numeric` `matrix` via [data.matrix()].
#' @param method A [`character`] string specifying the method to be
#'  used (see details). Any unambiguous substring can be given.
#' @param ... Further arguments to be passed to internal methods.
#' @details
#'  The following methods can be used to ascertain the degree of *turnover*
#'  in taxa composition along a gradient (\eqn{\beta}-diversity) on qualitative
#'  (presence/absence) data:
#'
#'  \describe{
#'   \item{`cody`}{[Cody measure][index_cody()].}
#'   \item{`routledge1`}{[Routledge first measure][index_routledge].}
#'   \item{`routledge2`}{[Routledge second measure][index_routledge].}
#'   \item{`routledge3`}{[Routledge third measure][index_routledge] (exponential
#'   form of the second measure).}
#'   \item{`whittaker`}{[Whittaker measure][index_whittaker()].}
#'   \item{`wilson`}{[Wilson measure][index_wilson()].}
#'  }
#'
#'  This assumes that the order of the matrix rows (from \eqn{1} to \eqn{n})
#'  follows the progression along the gradient/transect.
#' @return
#'  A [`numeric`] vector.
#' @seealso [index_cody()], [index_routledge1()], [index_routledge2()],
#'  [index_routledge3()], [index_whittaker()], [index_wilson()]
#' @example inst/examples/ex-turnover.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @aliases turnover-method
setGeneric(
  name = "turnover",
  def = function(object, ...) standardGeneric("turnover")
)

#' Cody Measure
#'
#' @param x A \eqn{m \times p}{m x p} `numeric` [`matrix`] of count data
#'  (absolute frequencies, i.e. a contingency table).
#' @param ... Currently not used.
#' @details
#'  This assumes that the order of the matrix rows (from \eqn{1} to \eqn{n})
#'  follows the progression along the gradient/transect.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Cody, M. L. (1975). Towards a theory of continental species diversity: Bird
#'  distributions over Mediterranean habitat gradients. *In* M. L. Cody &
#'  J. M. Diamond (Eds.), *Ecology and Evolution of Communities*.
#'  Cambridge, MA: Harvard University Press, p. 214-257.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_cody-method
setGeneric(
  name = "index_cody",
  def = function(x, ...) standardGeneric("index_cody")
)

#' Routledge Measures
#'
#' @param x A \eqn{m \times p}{m x p} `numeric` [`matrix`] of count data
#'  (absolute frequencies, i.e. a contingency table).
#' @param ... Currently not used.
#' @details
#'  This assumes that the order of the matrix rows (from \eqn{1} to \eqn{n})
#'  follows the progression along the gradient/transect.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Routledge, R. D. (1977). On Whittaker's Components of Diversity.
#'  *Ecology*, 58(5), 1120-1127. \doi{10.2307/1936932}.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @name index_routledge
#' @rdname index_routledge
NULL

#' @rdname index_routledge
#' @aliases index_routledge1-method
setGeneric(
  name = "index_routledge1",
  def = function(x, ...) standardGeneric("index_routledge1")
)

#' @rdname index_routledge
#' @aliases index_routledge2-method
setGeneric(
  name = "index_routledge2",
  def = function(x, ...) standardGeneric("index_routledge2")
)

#' @rdname index_routledge
#' @aliases index_routledge3-method
setGeneric(
  name = "index_routledge3",
  def = function(x, ...) standardGeneric("index_routledge3")
)

#' Whittaker Measure
#'
#' @param x A \eqn{m \times p}{m x p} `numeric` [`matrix`] of count data
#'  (absolute frequencies, i.e. a contingency table).
#' @param ... Currently not used.
#' @details
#'  This assumes that the order of the matrix rows (from \eqn{1} to \eqn{n})
#'  follows the progression along the gradient/transect.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Whittaker, R. H. (1960). Vegetation of the Siskiyou Mountains, Oregon and
#'  California. *Ecological Monographs*, 30(3), 279-338.
#'  \doi{10.2307/1943563}.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_whittaker-method
setGeneric(
  name = "index_whittaker",
  def = function(x, ...) standardGeneric("index_whittaker")
)

#' Wilson Measure
#'
#' @param x A \eqn{m \times p}{m x p} `numeric` [`matrix`] of count data
#'  (absolute frequencies, i.e. a contingency table).
#' @param ... Currently not used.
#' @details
#'  This assumes that the order of the matrix rows (from \eqn{1} to \eqn{n})
#'  follows the progression along the gradient/transect.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Wilson, M. V., & Shmida, A. (1984). Measuring Beta Diversity with
#'  Presence-Absence Data. *The Journal of Ecology*, 72(3), 1055-1064.
#'  \doi{10.2307/2259551}.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_wilson-method
setGeneric(
  name = "index_wilson",
  def = function(x, ...) standardGeneric("index_wilson")
)

## Co-Occurrence ---------------------------------------------------------------
#' Co-Occurrence
#'
#' @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 method A [`character`] string specifying the method to be
#'  used. It must be one of "`absolute`", "`relative`" or "`binomial`"
#'  (see details). Any unambiguous substring can be given.
#' @param ... Currently not used.
#' @details
#'  \describe{
#'   \item{`absolute`}{Count how many times each pairs of types occur together
#'   in at least one sample (absolute frequencies).}
#'   \item{`relative`}{Count how many times each pairs of types occur together
#'   in at least one sample (relative frequencies).}
#'   \item{`binomial`}{[Binomial co-occurrence assessment][index_binomial()].}
#'  }
#' @return
#'  A [stats::dist] object.
#' @example inst/examples/ex-occurrence.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @aliases occurrence-method
setGeneric(
  name = "occurrence",
  def = function(object, ...) standardGeneric("occurrence")
)

#' Binomial Co-Occurrence Assessment
#'
#' @param x,y A [`numeric`] vector.
#' @param ... Currently not used.
#' @details
#'  This assesses the degree of co-occurrence between taxa/types within a
#'  dataset. The strongest associations are shown by large positive numbers,
#'  the strongest segregations by large negative numbers.
#'
#'  The Binomial co-occurrence assessment approximates a Z-score.
#' @return
#'  A [`numeric`] vector.
#' @references
#'  Kintigh, K. (2006). Ceramic Dating and Type Associations. In J. Hantman and
#'  R. Most (eds.), *Managing Archaeological Data: Essays in Honor of
#'  Sylvia W. Gaines*. Anthropological Research Paper, 57. Tempe, AZ: Arizona
#'  State University, p. 17-26.
#' @author N. Frerebeau
#' @family beta diversity measures
#' @docType methods
#' @aliases index_binomial-method
setGeneric(
  name = "index_binomial",
  def = function(x, y, ...) standardGeneric("index_binomial")
)

## Simulate --------------------------------------------------------------------
#' Measure Diversity by Comparing to Simulated Assemblages
#'
#' @param object A [DiversityIndex-class] object.
#' @param interval A [`character`] string giving the type of confidence
#'  interval to be returned. It must be one "`percentiles`" (sample quantiles,
#'  as described in Kintigh 1984; the default), "`student`" or "`normal`".
#'  Any unambiguous substring can be given.
#' @param level A length-one [`numeric`] vector giving the confidence level.
#' @param step An [`integer`] giving the increment of the sample size.
#' @param n A non-negative [`integer`] giving the number of bootstrap
#'  replications.
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @return
#'  Returns a [DiversityIndex-class] object.
#' @references
#'  Baxter, M. J. (2001). Methodological Issues in the Study of Assemblage
#'  Diversity. *American Antiquity*, 66(4), 715-725. \doi{10.2307/2694184}.
#'
#'  Kintigh, K. W. (1984). Measuring Archaeological Diversity by Comparison
#'  with Simulated Assemblages. *American Antiquity*, 49(1), 44-54.
#'  \doi{10.2307/280511}.
#' @seealso [`plot()`][plot_diversity], [resample()]
#' @example inst/examples/ex-plot_diversity.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @name simulate
#' @rdname simulate
NULL

## Diversity Test --------------------------------------------------------------
#' Diversity Test
#'
#' Compares Shannon/Simpson diversity between samples.
#' @param x,y A [`numeric`] vector, a \eqn{m \times p}{m x p} [`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 adjust A [`character`] string specifying the method for
#'  adjusting \eqn{p} values (see [stats::p.adjust()]).
#' @param ... Further arguments to be passed to internal methods.
#' @return
#'  If `x` and `y` are `numeric` vectors, returns a [`list`] containing the
#'  following components:
#'  \describe{
#'   \item{`statistic`}{The value of the t-statistic.}
#'   \item{`parameter`}{The degrees of freedom for the t-statistic.}
#'   \item{`p.value`}{The p-value for the test.}
#'  }
#'
#'  If `x` is a `matrix` or a `data.frame`, returns a table of adjusted p-values
#'  in lower triangular form.
#' @example inst/examples/ex-test.R
#' @author N. Frerebeau
#' @references
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press. \doi{10.1007/978-94-015-7358-0}.
#' @family statistics
#' @docType methods
#' @name test
#' @rdname test
NULL

#' @rdname test
#' @aliases test_shannon-method
setGeneric(
  name = "test_shannon",
  def = function(x, y, ...) standardGeneric("test_shannon")
)

#' @rdname test
#' @aliases test_simpson-method
setGeneric(
  name = "test_simpson",
  def = function(x, y, ...) standardGeneric("test_simpson")
)

# Plot =========================================================================
## Diversity -------------------------------------------------------------------
#' Diversity Plot
#'
#' @param x A [DiversityIndex-class] object to be plotted.
#' @param log A [`character`] string indicating which axes should be in log
#'  scale. Defaults to `x`.
#' @param col.mean,col.interval A [`character`] string specifying the
#'  color of the lines.
#' @param lty.mean,lty.interval A [`character`] string or [`numeric`]
#'  value specifying the line types.
#' @param lwd.mean,lwd.interval A non-negative [`numeric`] value specifying
#'  the line widths.
#' @param xlab,ylab A [`character`] vector giving the x and y axis labels.
#' @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] to be passed to
#'  [graphics::points()], particularly, `cex`, `col` and `pch`.
#' @return
#'  `plot()` is called for its side-effects: it results in a graphic being
#'  displayed (invisibly returns `x`).
#' @example inst/examples/ex-plot_diversity.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @name plot_diversity
#' @rdname plot_diversity
NULL

#' Rarefaction Plot
#'
#' @param x A [RarefactionIndex-class] object to be plotted.
#' @param color A vector of colors (will be mapped to the rownames of `object`).
#'  If `color` is a named a named vector, then the colors will be associated
#'  with the rownames of `object`. Ignored if set to `FALSE`.
#' @param symbol A specification for the line type (will be mapped to
#'  the names of `x`). If `symbol` is a named a named vector, then the
#'  line types will be associated with the names of `x`.
#'  Ignored if set to `FALSE`.
#' @param xlab,ylab A [`character`] vector giving the x and y axis labels.
#' @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 legend A [`list`] of additional arguments to be passed to
#'  [graphics::legend()]; names of the list are used as argument names.
#'  If `NULL`, no legend is displayed.
#' @param ... Further [graphical parameters][graphics::par] to be passed to
#'  [graphics::lines()].
#' @return
#'  `plot()` is called for its side-effects: it results in a graphic being
#'  displayed (invisibly returns `x`).
#' @example inst/examples/ex-rarefaction.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @name plot_rarefaction
#' @rdname plot_rarefaction
NULL

#' SHE Analysis
#'
#' @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 unbiased A [`logical`] scalar: should the bias-corrected estimator be
#'  used (see [index_shannon()])?
#' @param xlab,ylab A [`character`] vector giving the x and y axis labels.
#' @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 legend A [`list`] of additional arguments to be passed to
#'  [graphics::legend()]; names of the list are used as argument names.
#'  If `NULL`, no legend is displayed.
#' @param ... Further [graphical parameters][graphics::par] to be passed to
#'  [graphics::lines()] and [graphics::points()].
#' @return
#'  `she()` is called for its side-effects: it results in a graphic being
#'  displayed (invisibly returns `object`).
#' @details
#'  If samples are taken along a gradient or stratigraphic section, breaks in
#'  the curve may be used to infer discontinuities.
#'
#'  This assumes that the order of the matrix rows (from \eqn{1} to \eqn{n})
#'  follows the progression along the gradient/transect.
#' @references
#'  Buzas, M. A. & Hayek, L.-A. C. (1998). SHE analysis for biofacies
#'  identification. *Journal of Foraminiferal Research*, 1998, 28(3), 233-239.
#'
#'  Hayek, L.-A. C. & Buzas, M. A. (2010). *Surveying Natural Populations:
#'  Quantitative Tools for Assessing Biodiversity*. Second edition.
#'  New York: Columbia University Press.
#' @example inst/examples/ex-she.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @aliases she-method
setGeneric(
  name = "she",
  def = function(object, ...) standardGeneric("she")
)

#' Diversity Profiles
#'
#' @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 alpha A [`numeric`] vector giving the values of the alpha parameter.
#' @param color A vector of colors (will be mapped to the rownames of `object`).
#'  If `color` is a named a named vector, then the colors will be associated
#'  with the rownames of `object`. Ignored if set to `FALSE`.
#' @param symbol A specification for the line type (will be mapped to
#'  the rownames of `object`). If `symbol` is a named a named vector, then the
#'  line types will be associated with the rownames of `object`.
#'  Ignored if set to `FALSE`.
#' @param xlab,ylab A [`character`] vector giving the x and y axis labels.
#' @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 legend A [`list`] of additional arguments to be passed to
#'  [graphics::legend()]; names of the list are used as argument names.
#'  If `NULL`, no legend is displayed.
#' @param ... Further [graphical parameters][graphics::par] to be passed to
#'  [graphics::lines()]
#' @details
#'  If the profiles cross, the diversities are non-comparable across samples.
#' @return
#'  `profiles()` is called for its side-effects: it results in a graphic being
#'  displayed (invisibly returns `object`).
#' @references
#'  Tóthmérész, B. (1995). Comparison of Different Methods for Diversity
#'  Ordering. *Journal of Vegetation Science*, 6(2), 283-290.
#'  \doi{10.2307/3236223}.
#' @example inst/examples/ex-profiles.R
#' @author N. Frerebeau
#' @family diversity measures
#' @docType methods
#' @aliases profiles-method
setGeneric(
  name = "profiles",
  def = function(object, ...) standardGeneric("profiles")
)

## Matrix plot -----------------------------------------------------------------
### Spot Plot ------------------------------------------------------------------
#' Spot Plot
#'
#' Plots a spot matrix.
#' @inheritParams plot_matrix
#' @param type A [`character`] string specifying the graph to be plotted.
#'  It must be one of "`ring`" (the default) or "`plain`". Any unambiguous
#'  substring can be given.
#' @param ... Currently not used.
#' @details
#'  The spot matrix can be considered as a variant of the
#'  [Bertin diagram][plot_bertin()] where the data are first transformed to
#'  relative frequencies.
#' @return
#'  `plot_spot()` is called for its side-effects: it results in a graphic
#'  being displayed (invisibly returns `object`).
#' @note
#'  Adapted from Dan Gopstein's original
#'  [idea](https://dan.gop/articles/spot-matrix/).
#' @example inst/examples/ex-plot_spot.R
#' @author N. Frerebeau
#' @family plot methods
#' @docType methods
#' @aliases plot_spot-method
setGeneric(
  name = "plot_spot",
  def = function(object, ...) standardGeneric("plot_spot")
)

### Heatmap --------------------------------------------------------------------
#' Heatmap
#'
#' Plots a heatmap.
#' @inheritParams plot_matrix
#' @param fixed_ratio A [`logical`] scalar: should a fixed aspect ratio (1) be
#' used?
#' @param ... Currently not used.
#' @return
#'  `plot_heatmap()` is called for its side-effects: it results in a graphic
#'  being displayed (invisibly returns `object`).
#' @example inst/examples/ex-plot_heatmap.R
#' @author N. Frerebeau
#' @family plot methods
#' @docType methods
#' @aliases plot_heatmap-method
setGeneric(
  name = "plot_heatmap",
  def = function(object, ...) standardGeneric("plot_heatmap")
)

### Matrigraph -----------------------------------------------------------------
#' Matrigraph
#'
#' @description
#'  * `matrigraph()` produces a heatmap highlighting the deviations from
#'    independence.
#'  * `pvi()` computes for each cell of a numeric matrix the percentage to the
#'    column theoretical independence value.
#' @inheritParams plot_matrix
#' @param reverse A [`logical`] scalar: should negative deviations be centered
#'  (see details)?
#' @param ... Currently not used.
#' @details
#'  PVI (in french "pourcentages de valeur d'indépendance") is calculated for
#'  each cell as the percentage to the column theoretical independence value:
#'  PVI greater than \eqn{1} represent positive deviations from the
#'  independence, whereas PVI smaller than \eqn{1} represent negative
#'  deviations (Desachy 2004).
#'
#'  The PVI matrix allows to explore deviations from independence (an
#'  intuitive approach to \eqn{\chi^2}{Chi-squared}), in such a way that a
#'  high-contrast matrix has quite significant deviations,
#'  with a low risk of being due to randomness (Desachy 2004).
#'
#'  `matrigraph()` displays the deviations from independence:
#'
#'  * If the PVI is equal to \eqn{1} (statistical independence), the cell of the
#'    matrix is filled in grey.
#'  * If the PVI is less than \eqn{1} (negative deviation from independence),
#'    the size of the grey square is proportional to the PVI (the white margin
#'    thus represents the fraction of negative deviation).
#'  * If the PVI is greater than \eqn{1} (positive deviation), a black
#'    square representing the fraction of positive deviations is
#'    superimposed. For large positive deviations (PVI greater than \eqn{2}),
#'    the cell in filled in black.
#'
#'  If `reverse` is `TRUE`, the fraction of negative deviations is displayed
#'  as a white square.
#' @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}.
#' @return
#'  * `matrigraph()` is called for its side-effects: it results in a graphic
#'    being displayed (invisibly returns `object`).
#'  * `pvi()` returns a [`numeric`] [`matrix`].
#' @example inst/examples/ex-matrigraph.R
#' @author N. Frerebeau
#' @seealso [plot_heatmap()]
#' @family plot methods
#' @docType methods
#' @aliases matrigraph-method
setGeneric(
  name = "matrigraph",
  def = function(object, ...) standardGeneric("matrigraph")
)

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

## Bar Plot --------------------------------------------------------------------
### Bertin ---------------------------------------------------------------------
#' Bertin Diagram
#'
#' Plots a Bertin diagram.
#' @inheritParams plot_matrix
#' @param threshold A [`function`] that takes a numeric vector as argument and
#'  returns a numeric threshold value (see below). If `NULL` (the default), no
#'  threshold is computed. Only used if `freq` is `FALSE`.
#' @param flip A [`logical`] scalar: should `x` and `y` axis be flipped?
#'  Defaults to `TRUE`.
#' @param ... Currently not used.
#' @section Bertin Matrix:
#'  As de Falguerolles *et al.* (1997) points out:
#'  "In abstract terms, a Bertin matrix is a matrix of  displays. [...] To fix
#'  ideas, think of a data matrix, variable by case, with real valued variables.
#'  For each variable, draw a bar chart of variable value by case. High-light
#'  all bars representing a value above some sample threshold for that
#'  variable."
#' @return
#'  `plot_bertin()` is called for its side-effects: it results in a graphic
#'  being displayed (invisibly returns `object`).
#' @references
#'  Bertin, J. (1977). *La graphique et le traitement graphique de
#'  l'information*. Paris: Flammarion. Nouvelle Bibliothèque Scientifique.
#'
#'  de Falguerolles, A., Friedrich, F. & Sawitzki, G. (1997). A Tribute to J.
#'  Bertin's Graphical Data Analysis. In W. Badilla & F. Faulbaum (eds.),
#'  *SoftStat '97: Advances in Statistical Software 6*. Stuttgart: Lucius
#'  & Lucius, p. 11-20.
#' @example inst/examples/ex-plot_bertin.R
#' @author N. Frerebeau
#' @family plot methods
#' @docType methods
#' @aliases plot_bertin-method Bertin
setGeneric(
  name = "plot_bertin",
  def = function(object, ...) standardGeneric("plot_bertin")
)

### Ford -----------------------------------------------------------------------
#' Ford Diagram
#'
#' Plots a Ford (battleship curve) 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).
#' @param weights A [`logical`] scalar: should the row sums be displayed?
#' @param EPPM A [`logical`] scalar: should the EPPM be drawn?
#'  See `seriograph()`.
#' @param fill The color for filling the bars.
#' @param border The color to draw the borders.
#' @param axes A [`logical`] scalar: should axes be drawn on the plot? It will
#'  omit labels where they would abut or overlap previously drawn labels.
#' @param ... Currently not used.
#' @return
#'  `plot_ford()` is called for its side-effects: it results in a graphic
#'  being displayed (invisibly returns `object`).
#' @references
#'  Ford, J. A. (1962). *A quantitative method for deriving cultural
#'  chronology*. Washington, DC: Pan American Union. Technical manual 1.
#' @example inst/examples/ex-plot_ford.R
#' @author N. Frerebeau
#' @family plot methods
#' @docType methods
#' @aliases plot_ford-method Ford
setGeneric(
  name = "plot_ford",
  def = function(object, ...) standardGeneric("plot_ford")
)

### Seriograph -----------------------------------------------------------------
#' Seriograph
#'
#' @description
#'  * `seriograph()` produces a Ford diagram highlighting the relationships
#'    between rows and columns.
#'  * `eppm()` computes for each cell of a numeric matrix the positive
#'    difference from the column mean percentage.
#' @inheritParams plot_ford
#' @details
#'  The positive difference from the column mean percentage (in french "écart
#'  positif au pourcentage moyen", EPPM) represents a deviation from the
#'  situation of statistical independence. As independence can be interpreted as
#'  the absence of relationships between types and the chronological order of
#'  the assemblages, EPPM is a useful tool to explore significance
#'  of relationship between rows and columns related to seriation (Desachy
#'  2004).
#'
#'  `seriograph()` superimposes the frequencies (grey) and EPPM values (black)
#'  for each row-column pair in a Ford diagram.
#' @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}.
#' @return
#'  * `seriograph()` is called for its side-effects: it results in a graphic
#'    being displayed (invisibly returns `object`).
#'  * `eppm()` returns a [`numeric`] [`matrix`].
#' @example inst/examples/ex-seriograph.R
#' @author N. Frerebeau
#' @seealso [plot_ford()]
#' @family plot methods
#' @docType methods
#' @aliases seriograph-method
setGeneric(
  name = "seriograph",
  def = function(object, ...) standardGeneric("seriograph")
)

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

### Dice-Leraas ----------------------------------------------------------------
#' Dice-Leraas Diagram
#'
#' Plots a Dice-Leraas 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 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].
#' @details
#'  In a Dice-Leraas diagram, the horizontal line represents the range of data
#'  (min-max) and the small vertical line indicates the mean. The black
#'  rectangle is twice the standard error on the mean, while the white rectangle
#'  is one standard deviation on either side of the mean.
#' @references
#'  Dice, L. R., & Leraas, H. J. (1936). A Graphic Method for Comparing Several
#'  Sets of Measurements. *Contributions from the Laboratory of Vertebrate
#'  Genetics*, 3: 1-3.
#'
#'  Hubbs, C. L., & C. Hubbs (1953). An Improved Graphical Analysis and
#'  Comparison of Series of Samples. *Systematic Biology*, 2(2): 49-56.
#'  \doi{10.2307/sysbio/2.2.49}.
#'
#'  Simpson, G. G., Roe, A., & Lewontin, R. C. *Quantitative Zoology*.
#'  New York: Harcourt, Brace and Company, 1960.
#' @return
#'  `plot_diceleraas()` is called for its side-effects: it results in a
#'  graphic being displayed (invisibly returns `object`).
#' @example inst/examples/ex-diceleraas.R
#' @author N. Frerebeau
#' @family plot methods
#' @docType methods
#' @aliases plot_diceleraas-method
setGeneric(
  name = "plot_diceleraas",
  def = function(object, ...) standardGeneric("plot_diceleraas")
)

## Line Plot -------------------------------------------------------------------
#' Rank Plot
#'
#' Plots a rank *vs* relative abundance 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 color A vector of colors (will be mapped to the rownames of `object`).
#'  If `color` is a named a named vector, then the colors will be associated
#'  with the rownames of `object`. Ignored if set to `FALSE`.
#' @param symbol A specification for the line type (will be mapped to
#'  the rownames of `object`). If `symbol` is a named a named vector, then the
#'  line types will be associated with the rownames of `object`.
#'  Ignored if set to `FALSE`.
#' @param log A [`character`] string which contains "`x`" if the x axis is to be
#'  logarithmic, "`y`" if the y axis is to be logarithmic and "`xy`" or "`yx`"
#'  if both axes are to be logarithmic (base 10).
# @param facet A [`logical`] scalar: should a matrix of panels defined by
#  case/sample be drawn?
#' @param xlab,ylab A [`character`] vector giving the x and y axis labels.
#' @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 legend A [`list`] of additional arguments to be passed to
#'  [graphics::legend()]; names of the list are used as argument names.
#'  If `NULL`, no legend is displayed.
#' @param ... Further [graphical parameters][graphics::par].
#' @return
#'  `plot_rank()` is called for its side-effects: it results in a graphic
#'  being displayed (invisibly returns `object`).
#' @references
#'  Magurran, A. E. (1988). *Ecological Diversity and its Measurement*.
#'  Princeton, NJ: Princeton University Press. \doi{10.1007/978-94-015-7358-0}.
#' @example inst/examples/ex-plot_rank.R
#' @author N. Frerebeau
#' @family plot methods
#' @docType methods
#' @aliases plot_rank-method
setGeneric(
  name = "plot_rank",
  def = function(object, ...) standardGeneric("plot_rank")
)

# SIMILARITY
#' @include AllGenerics.R
NULL

# Similarity ===================================================================
#' @export
#' @rdname similarity
#' @aliases similarity,matrix-method
setMethod(
  f = "similarity",
  signature = c(object = "matrix"),
  definition = function(object, method = c("brainerd", "bray", "jaccard",
                                           "morisita", "sorenson")) {
    ## Validation
    if (any(method == "binomial")) {
      stop("Use occurrence(x, method = \"bionomial\") instead.", call. = FALSE)
    }
    method <- match.arg(method, several.ok = FALSE)
    ## Select method
    fun <- get_index(method)

    ## Pairwise comparison
    ## Sample/case comparisons
    m <- nrow(object)
    labels <- rownames(object)
    beta <- function(x, f) f(object[x[1], ], object[x[2], ])

    cbn <- utils::combn(seq_len(m), m = 2, simplify = TRUE)
    index <- apply(X = cbn, MARGIN = 2, FUN = beta, f = fun)

    ## Matrix of results
    diag_value <- ifelse(method == "brainerd", 200, 1)
    mtx <- matrix(data = diag_value, nrow = m, ncol = m,
                  dimnames = list(labels, labels))
    mtx[lower.tri(mtx, diag = FALSE)] <- index
    mtx <- t(mtx)
    mtx[lower.tri(mtx, diag = FALSE)] <- index

    sim <- stats::as.dist(mtx)
    attr(sim, "method") <- method
    sim
  }
)

#' @export
#' @rdname similarity
#' @aliases similarity,data.frame-method
setMethod(
  f = "similarity",
  signature = c(object = "data.frame"),
  definition = function(object, method = c("brainerd", "bray", "jaccard",
                                           "morisita", "sorenson")) {
    object <- data.matrix(object)
    methods::callGeneric(object, method = method)
  }
)

# Index ========================================================================
## Qualitative index -----------------------------------------------------------
#' @export
#' @rdname index_jaccard
#' @aliases index_jaccard,character,character-method
setMethod(
  f = "index_jaccard",
  signature = c(x = "character", y = "character"),
  definition = function(x, y) {
    inter <- length(intersect(x, y))
    union <- length(x) + length(y) - inter
    inter / union
  }
)

#' @export
#' @rdname index_jaccard
#' @aliases index_jaccard,logical,logical-method
setMethod(
  f = "index_jaccard",
  signature = c(x = "logical", y = "logical"),
  definition = function(x, y) {
    ## Validation
    arkhe::assert_length(y, length(x))

    a <- sum(x)
    b <- sum(y)
    j <- sum((x + y) == 2)
    Cj <- j / (a + b - j)
    Cj
  }
)

#' @export
#' @rdname index_jaccard
#' @aliases index_jaccard,numeric,numeric-method
setMethod(
  f = "index_jaccard",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    ## presence/absence
    x <- x > 0
    y <- y > 0
    methods::callGeneric(x, y)
  }
)

#' @export
#' @rdname index_sorenson
#' @aliases index_sorenson,logical,logical-method
setMethod(
  f = "index_sorenson",
  signature = c(x = "logical", y = "logical"),
  definition = function(x, y) {
    ## Validation
    arkhe::assert_length(y, length(x))

    a <- sum(x)
    b <- sum(y)
    j <- sum((x + y) == 2)
    Cs <- 2 * j / (a + b)
    Cs
  }
)

#' @export
#' @rdname index_sorenson
#' @aliases index_sorenson,numeric,numeric-method
setMethod(
  f = "index_sorenson",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    ## presence/absence
    x <- x > 0
    y <- y > 0
    methods::callGeneric(x, y)
  }
)

## Quantitative index ----------------------------------------------------------
#' @export
#' @rdname index_brainerd
#' @aliases index_brainerd,numeric,numeric-method
setMethod(
  f = "index_brainerd",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    ## Validation
    arkhe::assert_length(y, length(x))

    a <- x / sum(x)
    b <- y / sum(y)
    Cb <- 2 - sum(abs(a - b))
    Cb * 100
  }
)

#' @export
#' @rdname index_bray
#' @aliases index_bray,numeric,numeric-method
setMethod(
  f = "index_bray",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    ## Validation
    arkhe::assert_length(y, length(x))

    a <- sum(x) # Number of individuals in site A
    b <- sum(y) # Number of individuals in site B

    Cs <- 1 - sum(abs(y - x)) / sum(x + y)
    Cs
  }
)

#' @export
#' @rdname index_morisita
#' @aliases index_morisita,numeric,numeric-method
setMethod(
  f = "index_morisita",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    ## Validation
    arkhe::assert_length(y, length(x))

    a <- sum(x) # Number of individuals in site A
    b <- sum(y) # Number of individuals in site B
    da <- sum(x^2) / a^2
    db <- sum(y^2) / b^2

    Cm <- (2 * sum(x * y)) / ((da + db) * a * b)
    Cm
  }
)

# DIVERSITY
#' @include AllGenerics.R
NULL

# Index ========================================================================
index_observed <- function(x, ...) {
  sum(x > 0, ...)  # Number of observed species
}
get_index <- function(x) {
  match.fun(sprintf("index_%s", x))
}
do_index <- function(x, method, ...) {
  f <- get_index(method)
  f(x, ...)
}

#' Compute a Diversity Index
#'
#' @param x A [`numeric`] [`matrix`].
#' @param method A [`character`] string specifying the measure to be computed.
#' @param by_row A [`logical`] scalar: should `method` be computed for each row?
#' @param ... Further parameters to be passed to `method`.
#' @return A [DiversityIndex-class] object.
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
index_diversity <- function(x, method, ..., by_row = TRUE) {
  fun <- get_index(method)
  if (by_row) {
    idx <- apply(X = x, MARGIN = 1, FUN = fun, ...)
  } else {
    idx <- fun(x, ...)
  }

  ## Fix names
  row_names <- rownames(x) %||% paste0("S", seq_along(idx))

  .DiversityIndex(
    idx,
    labels = row_names,
    size = as.integer(rowSums(x)),
    data = x,
    method = method
  )
}

# Heterogeneity ================================================================
#' @export
#' @rdname heterogeneity
#' @aliases heterogeneity,matrix-method
setMethod(
  f = "heterogeneity",
  signature = c(object = "matrix"),
  definition = function(object, ...,
                        method = c("berger", "boone", "brillouin",
                                   "mcintosh", "shannon", "simpson")) {
    method <- match.arg(method, several.ok = FALSE)
    by_row <- method != "boone"
    index <- index_diversity(object, method, ..., evenness = FALSE,
                             by_row = by_row)
    .HeterogeneityIndex(index)
  }
)

#' @export
#' @rdname heterogeneity
#' @aliases heterogeneity,data.frame-method
setMethod(
  f = "heterogeneity",
  signature = c(object = "data.frame"),
  definition = function(object, ...,
                        method = c("berger", "boone", "brillouin",
                                   "mcintosh", "shannon", "simpson")) {
    object <- data.matrix(object)
    methods::callGeneric(object, ..., method = method)
  }
)

# Evenness =====================================================================
#' @export
#' @rdname heterogeneity
#' @aliases evenness,matrix-method
setMethod(
  f = "evenness",
  signature = c(object = "matrix"),
  definition = function(object, ...,
                        method = c("shannon", "brillouin",
                                   "mcintosh", "simpson")) {
    method <- match.arg(method, several.ok = FALSE)
    index <- index_diversity(object, method, ..., evenness = TRUE)
    .EvennessIndex(index)
  }
)

#' @export
#' @rdname heterogeneity
#' @aliases evenness,data.frame-method
setMethod(
  f = "evenness",
  signature = c(object = "data.frame"),
  definition = function(object, ...,
                        method = c("shannon", "brillouin",
                                   "mcintosh", "simpson")) {
    object <- data.matrix(object)
    methods::callGeneric(object, ..., method = method)
  }
)

# Richness =====================================================================
#' @export
#' @rdname richness
#' @aliases richness,matrix-method
setMethod(
  f = "richness",
  signature = c(object = "matrix"),
  definition = function(object, ..., method = c("observed", "margalef", "menhinick")) {
    ## Backward compatibility
    if (method == "count") method <- "observed"

    method <- match.arg(method, several.ok = FALSE)
    index <- index_diversity(object, method, ...)
    .RichnessIndex(index)
  }
)

#' @export
#' @rdname richness
#' @aliases richness,data.frame-method
setMethod(
  f = "richness",
  signature = c(object = "data.frame"),
  definition = function(object, ..., method = c("observed", "margalef", "menhinick")) {
    object <- data.matrix(object)
    methods::callGeneric(object, ..., method = method)
  }
)

# Composition ==================================================================
#' @export
#' @rdname richness
#' @aliases composition,matrix-method
setMethod(
  f = "composition",
  signature = c(object = "matrix"),
  definition = function(object, ...,
                        method = c("chao1", "ace", "squares", "chao2", "ice")) {
    method <- match.arg(method, several.ok = FALSE)
    by_row <- any(method == c("chao1", "ace", "squares"))
    index <- index_diversity(object, method, ..., by_row = by_row)
    .CompositionIndex(index)
  }
)

#' @export
#' @rdname richness
#' @aliases composition,data.frame-method
setMethod(
  f = "composition",
  signature = c(object = "data.frame"),
  definition = function(object, ...,
                        method = c("chao1", "ace", "squares", "chao2", "ice")) {
    object <- data.matrix(object)
    methods::callGeneric(object, ..., method = method)
  }
)

# Turnover =====================================================================
#' @export
#' @rdname turnover
#' @aliases turnover,matrix-method
setMethod(
  f = "turnover",
  signature = c(object = "matrix"),
  definition = function(object, ...,
                        method = c("whittaker", "cody", "routledge1",
                                   "routledge2", "routledge3", "wilson")) {
    method <- match.arg(method, several.ok = FALSE)
    fun <- get_index(method)
    fun(object)
  }
)

#' @export
#' @rdname turnover
#' @aliases turnover,data.frame-method
setMethod(
  f = "turnover",
  signature = c(object = "data.frame"),
  definition = function(object, ...,
                        method = c("whittaker", "cody", "routledge1",
                                   "routledge2", "routledge3", "wilson")) {
    object <- data.matrix(object)
    methods::callGeneric(object, ..., method = method)
  }
)

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

    w <- object@data
    m <- nrow(w)
    method <- object@method

    results <- vector(mode = "list", length = m)
    for (i in seq_len(m)) {
      results[[i]] <- arkhe::bootstrap(
        object = w[i, ],
        do = do_index,
        n = n,
        method = method,
        evenness = methods::is(object, "EvennessIndex"),
        f = f
      )
    }
    results <- do.call(rbind, results)
    rownames(results) <- rownames(w)
    as.data.frame(results)
  }
)

## Jackknife -------------------------------------------------------------------
#' @export
#' @rdname jackknife
#' @aliases jackknife,DiversityIndex-method
setMethod(
  f = "jackknife",
  signature = c(object = "DiversityIndex"),
  definition = function(object, f = NULL) {

    w <- object@data
    m <- nrow(w)
    method <- object@method

    results <- vector(mode = "list", length = m)
    for (i in seq_len(m)) {
      results[[i]] <- arkhe::jackknife(
        object = w[i, ],
        do = do_index,
        method = method,
        evenness = methods::is(object, "EvennessIndex"),
        f = f
      )
    }
    results <- do.call(rbind, results)
    rownames(results) <- rownames(w)
    as.data.frame(results)
  }
)

## Simulate --------------------------------------------------------------------
#' @export
#' @rdname simulate
#' @aliases simulate,DiversityIndex-method
setMethod(
  f = "simulate",
  signature = c(object = "DiversityIndex"),
  definition = function(object, n = 1000, step = 1,
                        interval = c("percentiles", "student", "normal"),
                        level = 0.80, progress = getOption("tabula.progress")) {
    ## Simulate
    ## Specify the probability for the classes
    data <- object@data
    method <- object@method # Select method

    ## Sample size
    size <- max(rowSums(data))
    sample_sizes <- seq(from = 1, to = size * 1.05, by = step)

    m <- length(sample_sizes)
    k <- seq_len(m)

    simulated <- vector(mode = "list", length = m)
    fun <- function(x) conf(x, level = level, type = interval)

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

    for (i in k) {
      simulated[[i]] <- resample(
        object = colSums(data),
        do = do_index,
        method = method,
        evenness = methods::is(object, "EvennessIndex"),
        n = n,
        size = sample_sizes[[i]],
        f = fun
      )
      if (progress_bar) utils::setTxtProgressBar(pbar, i)
    }

    if (progress_bar) close(pbar)

    simulated <- do.call(rbind, simulated)
    simulated <- cbind(size = sample_sizes, simulated)

    methods::initialize(object, simulation = simulated)
  }
)

conf <- function(x, type = c("percentiles", "student", "normal"),
                 level = 0.80) {

  type <- match.arg(type, several.ok = FALSE)

  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
}

# PLOT BERTIN
#' @include AllGenerics.R
NULL

#' @export
#' @rdname plot_bertin
#' @aliases plot_bertin,matrix-method
setMethod(
  f = "plot_bertin",
  signature = signature(object = "matrix"),
  definition = function(object, threshold = NULL, freq = FALSE, margin = 1,
                        color = c("white", "black"), flip = TRUE,
                        axes = TRUE, ...) {
    ## Conditional proportions
    if (freq) object <- prop.table(object, margin = margin)

    ## Compute threshold for each variable
    if (!freq && is.function(threshold)) {
      thr <- apply(X = object, MARGIN = 2, FUN = threshold)
      thr <- matrix(thr, nrow = nrow(object), ncol = ncol(object), byrow = TRUE)
      thr <- ifelse(object > thr, color[length(color)], color[1L])
      color <- as.vector(t(thr))
    }

    ## /!\ Bertin plot flips x and y axis /!\
    data <- if (flip) t(object) else object

    plot_matrix(data, panel = panel_bertin, color = color,
                axes = axes, legend = FALSE, asp = NA)

    ## Legend
    # TODO

    invisible(object)
  }
)

#' @export
#' @rdname plot_bertin
#' @aliases plot_bertin,data.frame-method
setMethod(
  f = "plot_bertin",
  signature = signature(object = "data.frame"),
  definition = function(object, threshold = NULL, freq = FALSE, margin = 1,
                        color = c("white", "black"), flip = TRUE,
                        axes = TRUE, ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, threshold = threshold,
                         freq = freq, margin = margin,
                         color = color, flip = flip, axes = axes)
  }
)

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

#' @export
#' @rdname plot_rank
#' @aliases plot_rank,matrix-method
setMethod(
  f = "plot_rank",
  signature = c(object = "matrix"),
  definition = function(object, log = NULL,
                        color = NULL, symbol = FALSE,
                        xlab = NULL, ylab = NULL,
                        main = NULL, sub = NULL,
                        ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = axes,
                        panel.first = NULL, panel.last = NULL,
                        legend = list(x = "topright"), ...) {
    ## Prepare data
    n <- nrow(object)
    lab <- rownames(object)
    object <- object / rowSums(object)
    rk <- t(apply(X = object, MARGIN = 1, FUN = function(x) rank(-x)))

    ## Graphical parameters
    lwd <- list(...)$lwd %||% graphics::par("lwd")
    lty <- list(...)$lty %||% graphics::par("lty")
    cex <- list(...)$cex %||% graphics::par("cex")
    if (length(lwd) == 1) lwd <- rep(lwd, length.out = n)
    if (length(lty) == 1) lty <- rep(lty, length.out = n)
    if (length(cex) == 1) cex <- rep(cex, length.out = n)
    if (is.null(log)) log <- ""

    pch <- list(...)$pch %||% 16
    if (length(pch) == 1) pch <- rep(pch, length.out = n)
    if (!isFALSE(symbol)) pch <- khroma::palette_shape(symbol)(lab)

    col <- list(...)$col %||% graphics::par("col")
    if (length(col) == 1) col <- rep(col, length.out = n)
    if (!isFALSE(color)) col <- khroma::palette_color_discrete(color)(lab)

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

    ## Set plotting coordinates
    xlim <- range(rk)
    ylim <- range(object)
    graphics::plot.window(xlim = xlim, ylim = ylim, log = log)

    ## Evaluate pre-plot expressions
    panel.first

    ## Plot
    for (i in seq_len(n)) {
      z <- order(rk[i, ])
      graphics::lines(x = rk[i, z], y = object[i, z],
                      col = col[i], lty = lty[i], lwd = lwd[i])
      graphics::points(x = rk[i, z], y = object[i, z],
                       col = col[i], pch = pch[i], cex = cex[i])
    }

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

    ## Construct axis
    if (axes) {
      graphics::axis(side = 1, las = 1)
      graphics::axis(side = 2, las = 1)
    }

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

    ## Add annotation
    if (ann) {
      xlab <- xlab %||% tr_("Rank")
      ylab <- ylab %||% tr_("Frequency")
      graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
    }

    ## Legend
    if (is.list(legend) && length(legend) > 0) {
      args <- list(legend = lab, col = col, pch = pch,
                   lty = lty, lwd = lwd, bty = "n")
      args <- utils::modifyList(args, legend)
      do.call(graphics::legend, args = args)
    }

    invisible(object)
  }
)

#' @export
#' @rdname plot_rank
#' @aliases plot_rank,data.frame-method
setMethod(
  f = "plot_rank",
  signature = signature(object = "data.frame"),
  definition = function(object, log = NULL,
                        xlab = NULL, ylab = NULL,
                        main = NULL, sub = NULL,
                        ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = axes,
                        panel.first = NULL, panel.last = NULL,
                        legend = list(x = "topright"), ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, log = log,
                         xlab = xlab, ylab = ylab,
                         main = main, sub = sub, ann = ann, axes = axes,
                         frame.plot = frame.plot, panel.first = panel.first,
                         panel.last = panel.last, legend = legend, ...)
  }
)

# PLOT DICE-LERASS
#' @include AllClasses.R AllGenerics.R
NULL

#' @export
#' @rdname plot_diceleraas
#' @aliases plot_diceleraas,matrix-method
setMethod(
  f = "plot_diceleraas",
  signature = signature(object = "matrix"),
  definition = function(object,
                        main = NULL, sub = NULL,
                        ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = FALSE,
                        panel.first = NULL, panel.last = NULL, ...) {
    ## Prepare data
    object <- t(object)
    object[object == 0] <- NA
    n <- nrow(object)
    id <- seq_len(n)
    lab <- rownames(object) %||% id

    ## Mean
    moy <- rowMeans(object, na.rm = TRUE)
    ## Standard deviation
    ec <- apply(X = object, MARGIN = 1, FUN = stats::sd, na.rm = TRUE)
    ## Standard error
    se <- ec / sqrt(rowSums(!is.na(object)))
    ## Range
    ran <- t(apply(X = object, MARGIN = 1, FUN = range, na.rm = TRUE))
    colnames(ran) <- c("min", "max")

    data <- data.frame(y = rev(id), mean = moy, sd = ec, se = se * 2, ran)

    ## Graphical parameters
    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")
    lwd <- list(...)$lwd %||% graphics::par("lwd")
    lty <- list(...)$lty %||% graphics::par("lty")
    col <- list(...)$col %||% c("black")
    if (length(lwd) < n) lwd <- rep(lwd, length.out = n)
    if (length(lty) < n) lty <- rep(lty, length.out = n)
    if (length(col) < n) col <- rep(col, length.out = n)

    ## Save and restore
    mar <- graphics::par("mar")
    mar[2] <- inch2line(lab, 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
    xlim <- range(data$min, data$max, data$mean + data$sd,
                  data$mean - data$sd, na.rm = TRUE)
    ylim <- c(1, n + 1)
    graphics::plot.window(xlim = xlim, ylim = ylim)

    ## Evaluate pre-plot expressions
    panel.first

    ## Plot
    for (i in id) {
      tmp <- data[i, ]
      graphics::polygon(
        x = tmp$sd * c(-1, -1, 1, 1) + tmp$mean,
        y = tmp$y + (1 / 3) * c(0, 1, 1, 0),
        col = "white", border = NA
      )
      graphics::lines(
        x = tmp$sd * c(-1, -1, 1, 1) + tmp$mean,
        y = tmp$y + (1 / 3) * c(0, 1, 1, 0),
        col = col[i], lty = lty[i]
      )
      graphics::polygon(
        x = tmp$se * c(-1, -1, 1, 1) + tmp$mean,
        y = tmp$y + (1 / 3) * c(0, 1, 1, 0),
        col = col[i], border = NA
      )
      graphics::segments(
        x0 = tmp$mean, y0 = tmp$y,
        x1 = tmp$mean, y1 = tmp$y + 0.5,
        col = col[i], lty = lty[i]
      )
      graphics::segments(
        x0 = tmp$min, y0 = tmp$y,
        x1 = tmp$max, y1 = tmp$y,
        col = col[i], lty = lty[i]
      )
    }

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

    ## Construct axis
    if (axes) {
      graphics::axis(side = 1, las = 1)
      graphics::mtext(rev(lab), side = 2, at = id, las = 1, padj = 0, line = 0,
                      cex = cex.axis, col.axis = col.axis, font = font.axis)
    }

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

    ## Add annotation
    if (ann) {
      graphics::title(main = main, sub = sub,
                      xlab = tr_("Absolute frequency"), ylab = NULL, ...)
    }

    invisible(object)
  }
)

#' @export
#' @rdname plot_diceleraas
#' @aliases plot_diceleraas,data.frame-method
setMethod(
  f = "plot_diceleraas",
  signature = signature(object = "data.frame"),
  definition = function(object,
                        main = NULL, sub = NULL,
                        ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = FALSE,
                        panel.first = NULL, panel.last = NULL, ...) {
    object <- data.matrix(object)
    methods::callGeneric(object,
                         main = main, sub = sub, ann = ann,
                         axes = axes, frame.plot = frame.plot,
                         panel.first = panel.first, panel.last = panel.last, ...)
  }
)

# MATRIGRAPH
#' @include AllGenerics.R
NULL

#' @export
#' @rdname matrigraph
#' @aliases pvi,matrix-method
setMethod(
  f = "pvi",
  signature = signature(object = "matrix"),
  definition = function(object) {
    pvi <- object * 100 / expected(object)
    dimnames(pvi) <- dimnames(object)
    pvi
  }
)

#' @export
#' @rdname matrigraph
#' @aliases pvi,data.frame-method
setMethod(
  f = "pvi",
  signature = signature(object = "data.frame"),
  definition = function(object) {
    object <- data.matrix(object)
    methods::callGeneric(object)
  }
)

#' @export
#' @rdname matrigraph
#' @aliases matrigraph,matrix-method
setMethod(
  f = "matrigraph",
  signature = signature(object = "matrix"),
  definition = function(object, reverse = FALSE, axes = TRUE, ...) {
    plot_matrix(pvi(object), panel = panel_matrigraph, scale = FALSE,
                axes = axes, legend = FALSE, asp = 1, reverse = reverse)
    invisible(object)
  }
)

#' @export
#' @rdname matrigraph
#' @aliases matrigraph,data.frame-method
setMethod(
  f = "matrigraph",
  signature = signature(object = "data.frame"),
  definition = function(object, reverse = FALSE, ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, reverse = reverse, ...)
  }
)

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

#' @export
#' @rdname plot_spot
#' @aliases plot_spot,matrix-method
setMethod(
  f = "plot_spot",
  signature = signature(object = "matrix"),
  definition = function(object, type = c("ring", "plain"),
                        color = NULL,
                        diag = TRUE, upper = TRUE, lower = TRUE,
                        freq = FALSE, margin = 1,
                        axes = TRUE, legend = TRUE, ...) {
    ## Validation
    type <- match.arg(type, several.ok = FALSE)

    ## Backward compatibility
    col <- list(...)$col
    if (!is.null(col) && is.null(color)) color <- col

    plot_matrix(object, panel = panel_spot, color = color,
                 diag = diag, upper = upper, lower = lower,
                 freq = freq, margin = margin, drop_zero = TRUE,
                 axes = axes, legend = legend, type = type)

    invisible(object)
  }
)

#' @export
#' @rdname plot_spot
#' @aliases plot_spot,data.frame-method
setMethod(
  f = "plot_spot",
  signature = signature(object = "data.frame"),
  definition = function(object, type = c("ring", "plain"),
                        color = NULL,
                        diag = TRUE, upper = TRUE, lower = TRUE,
                        freq = FALSE, margin = 1,
                        axes = TRUE, legend = TRUE, ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, type = type, color = color,
                         diag = diag, upper = upper, lower = lower,
                         freq = freq, margin = margin,
                         axes = axes, legend = legend, ...)
  }
)

#' @export
#' @rdname plot_spot
#' @aliases plot_spot,dist-method
setMethod(
  f = "plot_spot",
  signature = signature(object = "dist"),
  definition = function(object, type = c("ring", "plain"),
                        color = NULL,
                        diag = FALSE, upper = FALSE, lower = !upper,
                        axes = TRUE, legend = TRUE, ...) {
    # index_name <- attr(object, "method")
    object <- as.matrix(object)
    methods::callGeneric(object, type = type, color = color,
                         diag = diag, upper = upper, lower = lower,
                         axes = axes, legend = legend)
  }
)

# PLOT HEATMAP
#' @include AllGenerics.R
NULL

#' @export
#' @rdname plot_heatmap
#' @aliases plot_heatmap,matrix-method
setMethod(
  f = "plot_heatmap",
  signature = signature(object = "matrix"),
  definition = function(object, color = NULL,
                        diag = TRUE, upper = TRUE, lower = TRUE,
                        freq = FALSE, margin = 1, fixed_ratio = TRUE,
                        axes = TRUE, legend = TRUE, ...) {
    ## Backward compatibility
    col <- list(...)$col
    if (!is.null(col) && is.null(color)) color <- col

    plot_matrix(object, panel = panel_tiles, color = color,
                 diag = diag, upper = upper, lower = lower,
                 freq = freq, margin = margin, drop_zero = FALSE,
                 axes = axes, legend = legend, asp = fixed_ratio)

    invisible(object)
  }
)

#' @export
#' @rdname plot_heatmap
#' @aliases plot_heatmap,data.frame-method
setMethod(
  f = "plot_heatmap",
  signature = signature(object = "data.frame"),
  definition = function(object, color = NULL,
                        diag = TRUE, upper = TRUE, lower = TRUE,
                        freq = FALSE, margin = 1, fixed_ratio = TRUE,
                        axes = TRUE, legend = TRUE, ...) {
    object <- data.matrix(object)
    methods::callGeneric(object, color = color,
                         diag = diag, upper = upper, lower = lower,
                         freq = freq, margin = margin,
                         fixed_ratio = fixed_ratio,
                         axes = axes, legend = legend)
  }
)

#' @export
#' @rdname plot_heatmap
#' @aliases plot_heatmap,dist-method
setMethod(
  f = "plot_heatmap",
  signature = signature(object = "dist"),
  definition = function(object, color = NULL,
                        diag = FALSE, upper = FALSE, lower = !upper,
                        axes = TRUE, legend = TRUE, ...) {
    object <- as.matrix(object)
    methods::callGeneric(object, color = color,
                         diag = diag, upper = upper, lower = lower,
                         axes = axes, legend = legend)
  }
)

# CO-OCCURRENCE
#' @include AllGenerics.R
NULL

#' @export
#' @rdname occurrence
#' @aliases occurrence,matrix-method
setMethod(
  f = "occurrence",
  signature = c(object = "matrix"),
  definition = function(object, method = c("absolute", "relative", "binomial")) {
    ## Validation
    method <- match.arg(method, several.ok = FALSE)

    ## Pairwise comparison
    p <- ncol(object)
    labels <- colnames(object)

    if (method == "absolute" || method == "relative") {
      incid <- object > 0
      fun <- function(x) sum(incid[, x[1]] + incid[, x[2]] == 2)
    }
    if (method == "binomial") {
      fun <- function(x) index_binomial(object[, x[1]], object[, x[2]])
    }

    cbn <- utils::combn(p, m = 2, simplify = TRUE)
    index <- apply(X = cbn, MARGIN = 2, FUN = fun)

    ## Matrix of results
    mtx <- matrix(data = 0, nrow = p, ncol = p,
                  dimnames = list(labels, labels))
    mtx[lower.tri(mtx, diag = FALSE)] <- index
    mtx <- t(mtx)
    mtx[lower.tri(mtx, diag = FALSE)] <- index

    if (method == "relative") {
      mtx <- mtx / nrow(object)
    }

    occ <- stats::as.dist(mtx)
    attr(occ, "total") <- nrow(object)
    occ
  }
)

#' @export
#' @rdname occurrence
#' @aliases occurrence,data.frame-method
setMethod(
  f = "occurrence",
  signature = c(object = "data.frame"),
  definition = function(object, method = c("absolute", "relative", "binomial")) {
    object <- data.matrix(object)
    methods::callGeneric(object, method = method)
  }
)

## Binomial co-occurrence ------------------------------------------------------
#' @export
#' @rdname index_binomial
#' @aliases index_binomial,numeric,numeric-method
setMethod(
  f = "index_binomial",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    ## Validation
    arkhe::assert_length(y, length(x))

    ## Total number of assemblages
    N <- length(x)
    ## Expected proportion of co-occurrences for artifact classes
    p <- sum(x > 0) * sum(y > 0) / N^2
    ## Number of observed co-occurence for artifact classes
    o <- sum((x > 0) + (y > 0) == 2)
    if (p == 1) {
      ## Avoid NaN generation
      ## TODO: warning ?
      Cbi <- 0
    } else {
      Cbi <- (o - N * p) / sqrt(N * p * (1 - p))
    }
    Cbi
  }
)

# COERCION
#' @include AllGenerics.R AllClasses.R
NULL

# To data.frame ================================================================
#' @method as.data.frame DiversityIndex
#' @export
as.data.frame.DiversityIndex <- function(x, ...) {
  data.frame(
    observed = apply(X = x@data, MARGIN = 1, FUN = observed),
    singleton = apply(X = x@data, MARGIN = 1, FUN = singleton),
    doubleton = apply(X = x@data, MARGIN = 1, FUN = doubleton),
    size = x@size,
    index = x@.Data,
    row.names = labels(x),
    stringsAsFactors = FALSE
  )
}

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

# SHOW METHODS
#' @include AllClasses.R
NULL

setMethod(
  f = "show",
  signature = "DiversityIndex",
  definition = function(object) {
    methods::callGeneric(object@.Data)
  }
)
setMethod(
  f = "show",
  signature = "RarefactionIndex",
  definition = function(object) {
    methods::callGeneric(object@.Data)
  }
)

1		#' @include AllGenerics.R AllClasses.R
2		NULL
3
4		# Shannon ======================================================================
5		variance_shannon <- function(x, base = exp(1), na.rm = FALSE, ...) {
6		## Validation
7	5x	x <- x[x > 0] # Remove unobserved species
8	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
9	!	if (anyNA(x)) return(NA)
10
11	5x	N <- sum(x)
12	5x	S <- length(x)
13	5x	p <- x / N
14
15	5x	a <- sum(p * (log(p, base = base))^2)
16	5x	b <- sum(p * log(p, base = base))^2
17
18	5x	var <- ((a - b) / N) + ((S - 1) / (2 * N^2))
19	5x	var
20		}
21
22		#' @export
23		#' @rdname test
24		#' @aliases test_shannon,numeric,numeric-method
25		setMethod(
26		f = "test_shannon",
27		signature = c(x = "numeric", y = "numeric"),
28		definition = function(x, y, ...) {
29		## Validation
30	2x	arkhe::assert_length(y, length(x))
31
32		## Calculate the number of individuals
33	2x	Nx <- sum(x, na.rm = TRUE)
34	2x	Ny <- sum(y, na.rm = TRUE)
35
36		## Calculate Shannon diversity
37		## See PAST documentation, p. 208
38	2x	d <- (sum(x > 0) - 1) / (2 * sum(x))
39	2x	Hx <- index_shannon(x, ...)
40	2x	Hy <- index_shannon(y, ...)
41
42		## Calculate Shannon variance
43	2x	Vx <- variance_shannon(x, ...)
44	2x	Vy <- variance_shannon(y, ...)
45
46		## t test statistic
47	2x	tt <- (Hx - Hy) / sqrt(Vx + Vy)
48
49		## Degrees of freedom
50	2x	df <- (Vx + Vy)^2 / sum(c(Vx, Vy)^2 / c(Nx, Ny))
51
52		## p value
53	2x	p <- 2 * (1 - stats::pt(q = abs(tt), df = df))
54
55	2x	list(
56	2x	statistic = tt,
57	2x	parameter = df,
58	2x	p.value = p
59		)
60		}
61		)
62
63		#' @export
64		#' @describeIn test Produces two sided pairwise comparisons.
65		#' @aliases test_shannon,matrix,missing-method
66		setMethod(
67		f = "test_shannon",
68		signature = c(x = "matrix", y = "missing"),
69		definition = function(x, adjust = "holm", ...) {
70		## Get the names of the assemblages
71	1x	row_names <- rownames(x)
72	1x	if (length(row_names) != 0) {
73	1x	row_names <- factor(row_names, levels = unique(row_names))
74		} else {
75	!	row_names <- factor(seq_len(nrow(x)))
76		}
77
78		## Compute t test
79	1x	compare <- function(i, j) {
80	1x	test_shannon(x[i, ], x[j, ], ...)$p.value
81		}
82
83	1x	result <- stats::pairwise.table(
84	1x	compare.levels = compare,
85	1x	level.names = row_names,
86	1x	p.adjust.method = adjust
87		)
88	1x	result
89		}
90		)
91
92		#' @export
93		#' @describeIn test Produces two sided pairwise comparisons.
94		#' @aliases test_shannon,data.frame,missing-method
95		setMethod(
96		f = "test_shannon",
97		signature = c(x = "data.frame", y = "missing"),
98		definition = function(x, adjust = "holm", ...) {
99	!	x <- data.matrix(x)
100	!	methods::callGeneric(x, adjust = adjust, ...)
101		}
102		)
103
104		# Simpson ======================================================================
105		variance_simpson <- function(x, na.rm = FALSE, ...) {
106		## Validation
107	4x	x <- x[x > 0] # Remove unobserved species
108	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
109	!	if (anyNA(x)) return(NA)
110
111	4x	N <- sum(x)
112	4x	S <- length(x)
113	4x	p <- x / N
114
115	4x	a <- 4 * N * (N - 1) * (N - 2) * sum(p^3)
116	4x	b <- 2 * N * (N - 1) * sum(p^2)
117	4x	c <- 2 * N * (N - 1) * (2 * N - 3) * sum(p^2)^2
118
119	4x	var <- (a + b - c) / (N^2 * (N - 1)^2)
120	4x	var
121		}
122
123		#' @export
124		#' @rdname test
125		#' @aliases test_simpson,numeric,numeric-method
126		setMethod(
127		f = "test_simpson",
128		signature = c(x = "numeric", y = "numeric"),
129		definition = function(x, y, adjust = "holm", ...) {
130		## Validation
131	2x	arkhe::assert_length(y, length(x))
132
133		## Calculate the number of individuals
134	2x	Nx <- sum(x, na.rm = TRUE)
135	2x	Ny <- sum(y, na.rm = TRUE)
136
137		## Calculate Shannon diversity
138	2x	Hx <- index_simpson(x, ...)
139	2x	Hy <- index_simpson(y, ...)
140
141		## Calculate Shannon variance
142	2x	Vx <- variance_simpson(x, ...)
143	2x	Vy <- variance_simpson(y, ...)
144
145		## t test statistic
146	2x	tt <- (Hx - Hy) / sqrt(Vx + Vy)
147
148		## Degrees of freedom
149	2x	df <- (Vx + Vy)^2 / sum(c(Vx, Vy)^2 / c(Nx, Ny))
150
151		## p value
152	2x	p <- 2 * (1 - stats::pt(q = abs(tt), df = df))
153
154	2x	list(
155	2x	statistic = tt,
156	2x	parameter = df,
157	2x	p.value = p
158		)
159		}
160		)
161
162		#' @export
163		#' @describeIn test Produces two sided pairwise comparisons.
164		#' @aliases test_simpson,matrix,missing-method
165		setMethod(
166		f = "test_simpson",
167		signature = c(x = "matrix", y = "missing"),
168		definition = function(x, adjust = "holm", ...) {
169		## Get the names of the assemblages
170	1x	row_names <- rownames(x)
171	1x	if (length(row_names) != 0) {
172	1x	row_names <- factor(row_names, levels = unique(row_names))
173		} else {
174	!	row_names <- factor(seq_len(nrow(x)))
175		}
176
177		## Compute t test
178	1x	compare <- function(i, j) {
179	1x	test_simpson(x[i, ], x[j, ])$p.value
180		}
181
182	1x	result <- stats::pairwise.table(
183	1x	compare.levels = compare,
184	1x	level.names = row_names,
185	1x	p.adjust.method = adjust
186		)
187	1x	result
188		}
189		)
190
191		#' @export
192		#' @describeIn test Produces two sided pairwise comparisons.
193		#' @aliases test_simpson,data.frame,missing-method
194		setMethod(
195		f = "test_simpson",
196		signature = c(x = "data.frame", y = "missing"),
197		definition = function(x, adjust = "holm", ...) {
198	!	x <- data.matrix(x)
199	!	methods::callGeneric(x, adjust = adjust)
200		}
201		)

1		# MATRIX PLOT
2
3		# Plot =========================================================================
4		#' Matrix Plot
5		#'
6		#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
7		#' [`data.frame`] of count data (absolute frequencies giving the number of
8		#' individuals for each category, i.e. a contingency table).
9		#' @param panel A [`function`] in the form `function(x, y, z, color, ...)`
10		#' which gives the action to be carried out in each panel of the display.
11		#' @param diag A [`logical`] scalar indicating whether the diagonal of the
12		#' matrix should be plotted. Only used if `object` is a symmetric matrix.
13		#' @param upper A [`logical`] scalar indicating whether the upper triangle of
14		#' the matrix should be plotted. Only used if `object` is a symmetric matrix.
15		#' @param lower A [`logical`] scalar indicating whether the lower triangle of
16		#' the matrix should be plotted. Only used if `object` is a symmetric matrix.
17		#' @param freq A [`logical`] scalar indicating whether conditional proportions
18		#' given `margins` should be used (i.e. entries of `object`, divided by the
19		#' appropriate marginal sums).
20		#' @param margin An [`integer`] vector giving the margins to split by:
21		#' `1` indicates individuals/rows (the default), `2` indicates
22		#' variables/columns. Only used if `freq` is `TRUE`.
23		#' @param scale A [`logical`] scalar indicating whether data should be rescaled
24		#' to \eqn{[-1,1]}. Only used if `freq` if `FALSE`.
25		#' @param drop_zero A [`logical`] scalar: should zeros be discarded?
26		#' @param color A vector of colors or a [`function`] that when called with a
27		#' single argument (an integer specifying the number of colors) returns a
28		#' vector of colors.
29		#' @param midpoint A [`numeric`] value specifying the data midpoint.
30		#' @param axes A [`logical`] scalar: should axes be drawn on the plot? It will
31		#' omit labels where they would abut or overlap previously drawn labels.
32		#' @param legend A [`logical`] scalar: should a legend be displayed?
33		#' @param asp A length-one [`numeric`] vector, giving the aspect ratio
34		#' \eqn{y/x}.
35		#' @param ... Further arguments to be passed to `panel`.
36		#' @keywords internal
37		plot_matrix <- function(object, panel, diag = TRUE, upper = TRUE, lower = TRUE,
38		freq = FALSE, margin = 1, scale = TRUE, drop_zero = TRUE,
39		color = graphics::par("fg"), midpoint = NULL,
40		axes = TRUE, legend = TRUE, asp = 1, ...) {
41		## Validation
42	21x	if (is_incidence(object)) {
43	1x	object[] <- as.numeric(object)
44	1x	legend <- FALSE
45		}
46
47		## Prepare data
48	21x	n <- nrow(object)
49	21x	m <- ncol(object)
50	21x	seq_row <- rev(seq_len(n))
51	21x	seq_col <- seq_len(m)
52	21x	lab_row <- rownames(object) %\|\|% seq_row
53	21x	lab_col <- colnames(object) %\|\|% seq_col
54
55	21x	data <- prepare(object, diag = diag, upper = upper, lower = lower,
56	21x	freq = freq, margin = margin, scale = scale,
57	21x	drop_zero = drop_zero,
58	21x	palette = color, midpoint = midpoint)
59
60		## Graphical parameters
61	21x	cex.axis <- graphics::par("cex.axis")
62	21x	col.axis <- graphics::par("col.axis")
63	21x	font.axis <- graphics::par("font.axis")
64
65		## Save and restore
66	21x	d <- inch2line("M", cex = cex.axis)
67	21x	old_par <- graphics::par("mar", "plt")
68	21x	on.exit(graphics::par(old_par))
69
70	21x	mar_left <- inch2line(lab_row, cex = cex.axis)
71	21x	mar_top <- inch2line(lab_col, cex = cex.axis)
72	21x	mar_right <- if (legend) inch2line("999%", cex = cex.axis) else d
73	21x	graphics::par(mar = c(d, mar_left, mar_top, mar_right))
74
75		## Open new window
76	21x	grDevices::dev.hold()
77	21x	on.exit(grDevices::dev.flush(), add = TRUE)
78	21x	graphics::plot.new()
79
80		## Squish plotting area
81	21x	pin <- graphics::par("pin")
82	21x	plt <- graphics::par("plt")
83
84		## Add horizontal space for the legend (5%)
85	21x	legend_width <- if (legend) max(1, m / 20) else 0
86
87	6x	if (isTRUE(asp)) asp <- 1
88	21x	if (!isFALSE(asp) && !is.na(asp)) {
89	15x	aspect_ratio <- n / (m + legend_width)
90	15x	pin_y <- pin[1] * aspect_ratio * asp
91
92	15x	if (pin_y < pin[2]) {
93		## Squish vertically
94	11x	graphics::par(pin = c(pin[1], pin_y))
95	11x	graphics::par(plt = c(plt[1:2], graphics::par('plt')[3:4]))
96		} else {
97		## Squish horizontally
98	4x	pin_x <- pin[2] / aspect_ratio / asp
99	4x	graphics::par(pin = c(pin_x, pin[2]))
100	4x	graphics::par(plt = c(graphics::par('plt')[1:2], plt[3:4]))
101		}
102		}
103
104		## Set plotting coordinates
105	21x	xlim <- c(0, m + legend_width) + 0.5
106	21x	ylim <- c(0, n) + 0.5
107	21x	graphics::plot.window(xlim = xlim, ylim = ylim, xaxs = "i", yaxs = "i", asp = asp)
108
109		## Plot
110	21x	panel(x = data$x, y = data$y, z = data$scaled, color = data$color, ...)
111
112		## Construct axis
113	21x	if (axes) {
114	21x	graphics::axis(side = 2, at = seq_row, labels = lab_row, las = 2,
115	21x	lty = 0, cex.axis = cex.axis, col.axis = col.axis,
116	21x	font.axis = font.axis)
117	21x	graphics::axis(side = 3, at = seq_col, labels = lab_col, las = 2,
118	21x	lty = 0, cex.axis = cex.axis, col.axis = col.axis,
119	21x	font.axis = font.axis)
120		}
121
122		## Legend
123	21x	if (legend) {
124	13x	lgd <- attr(data, "legend")
125	13x	legend_gradient(x = m, y = n, labels = lgd$labels,
126	13x	at = lgd$at, width = legend_width, col = lgd$colors)
127		}
128		}
129
130		legend_gradient <- function(x, y, labels, at, width, col) {
131	13x	legend_image <- grDevices::as.raster(col)
132	13x	legend_y <- (at - min(at)) * y / diff(range(at)) + 0.5
133
134	13x	graphics::rasterImage(legend_image, xleft = x + 1, ybottom = max(legend_y),
135	13x	xright = x + 1 + width, ytop = min(legend_y))
136	13x	graphics::segments(x0 = x + 1, y0 = legend_y,
137	13x	x1 = x + 1 + width, y1 = legend_y,
138	13x	col = "white")
139	13x	graphics::polygon(x = c(x, x + width, x + width, x) + 1,
140	13x	y = c(0.5, 0.5, max(legend_y), max(legend_y)),
141	13x	col = NA, border = "black")
142	13x	graphics::axis(side = 4, at = legend_y, labels = labels, las = 2)
143		}
144
145		# Panels =======================================================================
146		panel_bertin <- function(x, y, z, color, ..., space = 0.05) {
147	5x	y_bottom <- y - 0.5
148	5x	y_top <- y - 0.5 + z - space * 2
149	5x	y_top[y_top < y_bottom] <- y_bottom[y_top < y_bottom]
150	5x	graphics::rect(
151	5x	xleft = x - 0.5 + space,
152	5x	ybottom = y_bottom,
153	5x	xright = x + 0.5 - space,
154	5x	ytop = y_top,
155	5x	col = color,
156	5x	border = "black"
157		)
158		}
159		panel_matrigraph <- function(x, y, z, ..., reverse = FALSE) {
160	2x	pvi <- data.frame(x = x, y = y, z = z)
161	2x	pvi$z <- pvi$z / 100
162	2x	pvi_plus <- pvi[pvi$z > 1, ]
163	2x	pvi_plus$z <- pvi_plus$z - 1
164	2x	pvi_plus$z[pvi_plus$z > 1] <- 1
165
166	2x	if (reverse) {
167	1x	pvi_minus <- pvi[pvi$z < 1, ]
168	1x	pvi_minus$z <- 1 - pvi_minus$z
169		} else {
170	1x	pvi_minus <- pvi
171	1x	pvi_minus$z[pvi_minus$z > 1] <- 1
172		}
173
174	2x	pvi_minus$z <- pvi_minus$z * 0.5
175	2x	pvi_plus$z <- pvi_plus$z * 0.5
176
177	2x	col_bkg <- if (reverse) "darkgrey" else "white"
178	2x	col_minus <- if (reverse) "white" else "darkgrey"
179
180	2x	graphics::rect(xleft = pvi$x - 0.5, ybottom = pvi$y - 0.5,
181	2x	xright = pvi$x + 0.5, ytop = pvi$y + 0.5,
182	2x	col = col_bkg, border = NA)
183	2x	graphics::rect(xleft = pvi_minus$x - pvi_minus$z,
184	2x	ybottom = pvi_minus$y - pvi_minus$z,
185	2x	xright = pvi_minus$x + pvi_minus$z,
186	2x	ytop = pvi_minus$y + pvi_minus$z,
187	2x	col = col_minus, border = NA)
188	2x	graphics::rect(xleft = pvi_plus$x - pvi_plus$z,
189	2x	ybottom = pvi_plus$y - pvi_plus$z,
190	2x	xright = pvi_plus$x + pvi_plus$z,
191	2x	ytop = pvi_plus$y + pvi_plus$z,
192	2x	col = "black", border = NA)
193		}
194		panel_tiles <- function(x, y, color, ...) {
195	7x	width <- 0.5
196	7x	graphics::rect(
197	7x	xleft = x - width,
198	7x	ybottom = y - width,
199	7x	xright = x + width,
200	7x	ytop = y + width,
201	7x	col = color,
202	7x	border = NA
203		)
204		}
205		panel_spot <- function(x, y, z, color, type, ...) {
206	7x	radius <- abs(z * 0.45)
207	7x	for (i in seq_along(x)) {
208	185x	circle(x = x[i], y = y[i], radius = radius[i],
209	185x	col = color[i], border = color[i])
210	185x	if (type == "ring") {
211	123x	circle(x = x[i], y = y[i], radius = 0.45,
212	123x	col = NA, border = "black")
213		}
214		}
215		}
216
217		# Prepare ======================================================================
218		#' Prepare Data for Matrix Plot
219		#'
220		#' @param object A \eqn{m \times p}{m x p} `numeric` [`matrix`] or
221		#' [`data.frame`] of count data (absolute frequencies giving the number of
222		#' individuals for each category, i.e. a contingency table).
223		#' @param diag A [`logical`] scalar indicating whether the diagonal of the
224		#' matrix should be plotted. Only used if `object` is a symmetric matrix.
225		#' @param upper A [`logical`] scalar indicating whether the upper triangle of
226		#' the matrix should be plotted. Only used if `object` is a symmetric matrix.
227		#' @param lower A [`logical`] scalar indicating whether the lower triangle of
228		#' the matrix should be plotted. Only used if `object` is a symmetric matrix.
229		#' @param freq A [`logical`] scalar indicating whether conditional proportions
230		#' given `margins` should be used (i.e. entries of `object`, divided by the
231		#' appropriate marginal sums).
232		#' @param margin An [`integer`] vector giving the margins to split by:
233		#' `1` indicates individuals/rows (the default), `2` indicates
234		#' variables/columns. Only used if `freq` is `TRUE`.
235		#' @param scale A [`logical`] scalar indicating whether data should be rescaled
236		#' to \eqn{[-1,1]}. Only used if `freq` if `FALSE`.
237		#' @param drop_zero A [`logical`] scalar: should zeros be discarded?
238		#' @param palette A vector of colors.
239		#' @param midpoint A [`numeric`] value specifying the data midpoint.
240		#' @param ... Currently not used.
241		#' @return
242		#' A long [`data.frame`].
243		#' \describe{
244		#' \item{`row`}{}
245		#' \item{`column`}{}
246		#' \item{`x`,`y`}{Tile center coordinates.}
247		#' \item{`z`}{Raw values.}
248		#' \item{`value`}{}
249		#' \item{`scaled`}{}
250		#' \item{`color`}{}
251		#' }
252		#' @keywords internal
253		#' @noRd
254		prepare <- function(object, diag = TRUE, upper = TRUE, lower = TRUE,
255		freq = FALSE, margin = 1, scale = !freq, drop_zero = FALSE,
256		palette = grDevices::hcl.colors(12, "YlOrBr", rev = TRUE),
257		midpoint = NULL, ...) {
258		## Validation
259	21x	if (!arkhe::is_symmetric(object)) {
260	16x	diag <- TRUE
261	16x	upper <- TRUE
262	16x	lower <- TRUE
263		}
264
265		## Coerce to matrix
266	21x	object <- as.matrix(object)
267
268		## Relative frequencies
269	21x	val <- if (freq) prop.table(object, margin = margin) else object
270
271		## Rescale to [-1;1]
272	21x	sca <- if (scale) val / max(abs(val), na.rm = TRUE) else val
273
274	21x	val <- as.vector(val)
275	21x	sca <- as.vector(sca)
276	21x	min_val <- min(val, na.rm = TRUE)
277	21x	max_val <- max(val, na.rm = TRUE)
278
279		## Build a long table
280	21x	row <- row(object, as.factor = TRUE)
281	21x	col <- col(object, as.factor = TRUE)
282	21x	data <- data.frame(
283	21x	row = as.vector(row),
284	21x	column = as.vector(col),
285	21x	x = as.numeric(col),
286	21x	y = as.vector(nrow(object) - as.numeric(row) + 1), # Reverse y for plotting
287	21x	z = as.vector(object),
288	21x	value = val,
289	21x	scaled = sca
290		)
291
292		## Map colors
293	21x	breaks <- pretty(val, n = 5)
294	21x	domain <- range(c(breaks, min_val, max_val))
295	21x	midpoint <- if (is.null(midpoint) & min_val < 0 & max_val > 0) 0 else midpoint
296	21x	pal <- khroma::palette_color_continuous(colors = palette, domain = domain, midpoint = midpoint)
297	21x	data$color <- if (length(palette) == length(val)) palette else pal(val)
298
299		## Clean data
300	3x	if (!upper) data <- data[!upper.tri(object), ]
301	2x	if (!lower) data <- data[!lower.tri(object), ]
302	5x	if (!diag) data <- data[data$row != data$column, ]
303	14x	if (drop_zero) data <- data[data$value != 0, ]
304
305		## Legend
306	21x	attr(data, "legend") <- list(
307	21x	labels = if (freq) label_percent(breaks) else breaks,
308	21x	at = breaks / max(abs(val), na.rm = TRUE),
309	21x	colors = pal(breaks)
310		)
311	21x	data
312		}

1		# PLOT FORD
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname plot_ford
7		#' @aliases plot_ford,matrix-method
8		setMethod(
9		f = "plot_ford",
10		signature = signature(object = "matrix"),
11		definition = function(object, weights = FALSE, EPPM = FALSE,
12		fill = "darkgrey", border = NA,
13		axes = TRUE, ...) {
14		## Prepare data
15	3x	n <- nrow(object)
16	3x	m <- ncol(object)
17	3x	seq_row <- rev(seq_len(n))
18	3x	seq_col <- seq_len(m)
19	3x	lab_row <- rownames(object) %\|\|% seq_row
20	3x	lab_col <- colnames(object) %\|\|% seq_col
21
22	3x	padding_x <- 0.05
23	3x	padding_y <- 0.5 - 0.01
24	3x	data <- prepare_ford(object, padding = padding_x)
25
26		## Graphical parameters
27	3x	cex.axis <- graphics::par("cex.axis")
28	3x	col.axis <- graphics::par("col.axis")
29	3x	font.axis <- graphics::par("font.axis")
30
31		## Save and restore
32	3x	d <- inch2line("M", cex = cex.axis)
33	3x	mfrow <- graphics::par("mfrow")
34	3x	mar <- graphics::par("mar")
35	3x	mar[1] <- 3
36	3x	mar[2] <- inch2line(lab_row, cex = cex.axis) + 0.1
37	3x	mar[3] <- inch2line(lab_col, cex = cex.axis) + 0.1
38	3x	mar[4] <- 0
39
40	3x	old_par <- graphics::par(mfrow = mfrow, mar = mar)
41	3x	on.exit(graphics::par(old_par))
42
43	3x	if (weights) {
44	1x	graphics::layout(matrix(c(1, 2), nrow = 1), widths = c(m - 1, 1))
45		}
46
47		## Open new window
48	3x	grDevices::dev.hold()
49	3x	on.exit(grDevices::dev.flush(), add = TRUE)
50	3x	graphics::plot.new()
51
52		## Set plotting coordinates
53	3x	xlim <- range(data$x - data$value, data$x + data$value)
54	3x	ylim <- range(data$y) + c(-1, 1) * padding_y
55	3x	graphics::plot.window(xlim = xlim, ylim = ylim)
56
57		## Plot
58	3x	graphics::rect(
59	3x	xleft = data$x - data$value,
60	3x	ybottom = data$y - padding_y,
61	3x	xright = data$x + data$value,
62	3x	ytop = data$y + padding_y,
63	3x	col = fill,
64	3x	border = border
65		)
66	3x	if (EPPM) {
67	1x	graphics::rect(
68	1x	xleft = data$x - data$eppm,
69	1x	ybottom = data$y - padding_y,
70	1x	xright = data$x + data$eppm,
71	1x	ytop = data$y + padding_y,
72	1x	col = "black",
73	1x	border = NA
74		)
75		}
76
77		## Construct axis
78	3x	if (axes) {
79	3x	graphics::axis(side = 2, at = seq_row, labels = lab_row, las = 2,
80	3x	lty = 0, cex.axis = cex.axis, col.axis = col.axis,
81	3x	font.axis = font.axis)
82	3x	graphics::axis(side = 3, at = unique(data$x), labels = lab_col, las = 2,
83	3x	lty = 0, cex.axis = cex.axis, col.axis = col.axis,
84	3x	font.axis = font.axis)
85
86	3x	x_axis <- data$x[which.max(data$value)]
87	3x	graphics::axis(side = 1, at = c(x_axis - 0.2, x_axis + 0.2), labels = FALSE)
88	3x	graphics::axis(side = 1, at = x_axis, labels = label_percent(0.2),
89	3x	tick = FALSE)
90		}
91
92	3x	if (weights) {
93		## Graphical parameters
94	1x	mar[2] <- 0.1
95	1x	mar[4] <- 0.1
96	1x	graphics::par(mar = mar)
97
98		## Open new window
99	1x	graphics::plot.new()
100
101		## Set plotting coordinates
102	1x	total <- rowSums(object)
103	1x	graphics::plot.window(xlim = c(0, max(total) * 1.05), ylim = ylim,
104	1x	xaxs = "i")
105
106		## Plot
107	1x	y <- length(total) - seq_along(total) + 1
108	1x	graphics::rect(
109	1x	xleft = 0,
110	1x	ybottom = y - padding_y,
111	1x	xright = total,
112	1x	ytop = y + padding_y,
113	1x	col = fill,
114	1x	border = border
115		)
116
117		## Construct axis
118	1x	if (axes) {
119	1x	graphics::segments(x0 = 0, y0 = 0, x1 = 0, y = n + 0.5,
120	1x	col = col.axis, lwd = 1)
121	1x	graphics::axis(side = 1)
122		}
123		}
124
125	3x	invisible(object)
126		}
127		)
128
129		#' @export
130		#' @rdname plot_ford
131		#' @aliases plot_ford,data.frame-method
132		setMethod(
133		f = "plot_ford",
134		signature = signature(object = "data.frame"),
135		definition = function(object, weights = FALSE, EPPM = FALSE,
136		fill = "darkgrey", border = NA,
137		axes = TRUE, ...) {
138	2x	object <- data.matrix(object)
139	2x	methods::callGeneric(object, weights = weights, EPPM = EPPM, fill = fill,
140	2x	border = border, axes = axes, ...)
141		}
142		)
143
144		#' Prepare data for Ford plot
145		#' @return A data.frame.
146		#' @keywords internal
147		#' @noRd
148		prepare_ford <- function(x, padding = 0.05) {
149		## EPPM
150	3x	EPPM <- eppm(x) / 100
151
152		## Relative frequencies
153	3x	freq <- x / rowSums(x)
154
155		## Adaptive spacing between columns
156	3x	col_max <- apply(X = freq, MARGIN = 2, FUN = max, na.rm = TRUE)
157	3x	roll_max <- roll_sum(col_max, n = 2) + padding * max(freq)
158	3x	cum_max <- c(0, cumsum(roll_max))
159
160		## Build a long table
161	3x	row <- row(x, as.factor = TRUE)
162	3x	col <- col(x, as.factor = TRUE)
163	3x	data <- data.frame(
164	3x	row = as.vector(row),
165	3x	column = as.vector(col),
166	3x	value = as.vector(freq),
167	3x	eppm = as.vector(EPPM)
168		)
169
170	3x	m <- nrow(freq)
171	3x	data$x <- rep(cum_max, each = m) + 1
172	3x	data$y <- as.vector(m + 1 - as.numeric(row)) # Reverse levels order
173
174	3x	return(data)
175		}

1		# STATISTICS
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		expected <- function(x) {
6	7x	m <- nrow(x)
7	7x	p <- ncol(x)
8
9	7x	column_total <- matrix(colSums(x), nrow = m, ncol = p, byrow = TRUE)
10	7x	row_total <- matrix(rowSums(x), nrow = m, ncol = p, byrow = FALSE)
11	7x	grand_total <- sum(x)
12
13	7x	column_total * row_total / grand_total
14		}
15
16		#' @export
17		#' @rdname resample
18		#' @aliases resample,numeric-method
19		setMethod(
20		f = "resample",
21		signature = c(object = "numeric"),
22		definition = function(object, do, n, size = sum(object), ..., f = NULL) {
23		## Validation
24	318x	arkhe::assert_count(object)
25
26	318x	prob <- object / sum(object)
27	318x	replicates <- stats::rmultinom(n, size = size, prob = prob)
28	318x	values <- apply(X = replicates, MARGIN = 2, FUN = do, ...)
29	318x	if (is.function(f)) values <- f(values)
30	318x	values
31		}
32		)
33
34		#' Binomial Coefficient
35		#'
36		#' Computes the number of `k`-combinations from a given set of `n` elements
37		#' ("`n` choose `k`").
38		#' @param n A length-one [`numeric`] vector.
39		#' @param k A length-one [`numeric`] vector.
40		#' @details
41		#' Ramanujan approximation is used for \eqn{x!} computation if \eqn{x > 170}.
42		#' @return A length-one [`numeric`] vector.
43		#' @author N. Frerebeau
44		#' @keywords internal
45		#' @noRd
46		combination <- function(n, k) {
47		# Validation
48	373x	arkhe::assert_scalar(n, "numeric")
49	373x	arkhe::assert_scalar(k, "numeric")
50
51	372x	if (n > 170 \| k > 170) {
52		## Ramanujan approximation of x!
53	1x	c <- exp(ramanujan(n) - ramanujan(k) - ramanujan(n - k))
54		} else {
55	371x	c <- factorial(n) / (factorial(k) * factorial(n - k))
56		}
57	372x	c
58		}
59
60		#' Ramanujan Factorial Approximation
61		#'
62		#' @param x A [`numeric`] vector.
63		#' @return A [`numeric`] vector.
64		#' @examples
65		#' factorial(50)
66		#' exp(ramanujan(50))
67		#' @references
68		#' Ramanujan Aiyangar, S. (1988). *The lost notebook and other unpublished
69		#' papers*. Berlin: Springer-Verlag.
70		#' @author N. Frerebeau
71		#' @keywords internal
72		#' @noRd
73		ramanujan <- function(x){
74	467x	x * log(x) - x + log(x * (1 + 4 * x * (1 + 2 * x))) / 6 + log(pi) / 2
75		}

1		# BETA DIVERSITY
2		#' @include AllGenerics.R
3		NULL
4
5		# Cody =========================================================================
6		#' @export
7		#' @rdname index_cody
8		#' @aliases index_cody,matrix-method
9		setMethod(
10		f = "index_cody",
11		signature = signature(x = "matrix"),
12		definition = function(x) {
13	2x	x <- x > 0 # presence/absence
14
15	2x	m <- nrow(x)
16	2x	first_transect <- x[1, ]
17	2x	last_transect <- x[m, ]
18
19		# Number of taxa gained along the transect
20	2x	absent <- x[, !first_transect] # Keep only taxa absent in the first transect
21	2x	gained <- sum(apply(X = absent, MARGIN = 2, FUN = cumsum)[m, ] > 0)
22		# Number of taxa lost
23	2x	lost <- sum(!last_transect)
24
25	2x	beta <- (gained + lost) / 2
26	2x	beta
27		}
28		)
29
30		# Routledge ====================================================================
31		## Beta R ----------------------------------------------------------------------
32		#' @export
33		#' @rdname index_routledge
34		#' @aliases index_routledge1,matrix-method
35		setMethod(
36		f = "index_routledge1",
37		signature = signature(x = "matrix"),
38		definition = function(x) {
39	2x	x <- x > 0 # presence/absence
40
41	2x	p <- ncol(x)
42
43		# Total number of taxa recorded in the system
44	2x	S <- sum(colSums(x) > 0)
45		# Works out which pairs of taxa occurs together in at least one sample
46	2x	combine <- as.matrix(occurrence(x)) > 0
47	2x	r <- sum(combine) / 2
48
49	2x	beta <- S^2 / (2 * r + S) - 1
50	2x	return(beta)
51		}
52		)
53
54		## Beta I ----------------------------------------------------------------------
55		#' @export
56		#' @rdname index_routledge
57		#' @aliases index_routledge2,matrix-method
58		setMethod(
59		f = "index_routledge2",
60		signature = signature(x = "matrix"),
61		definition = function(x) {
62	4x	x <- x > 0 # presence/absence
63
64		# Number of samples in which each taxa is present
65	4x	e <- colSums(x)
66		# Taxa richness of each sample
67	4x	alpha <- rowSums(x)
68	4x	t <- sum(e)
69
70	4x	beta <- log(t) - (1 / t) * sum(e * log(e)) - (1 / t) * sum(alpha * log(alpha))
71	4x	return(beta)
72		}
73		)
74
75		## Beta E ----------------------------------------------------------------------
76		#' @export
77		#' @rdname index_routledge
78		#' @aliases index_routledge3,matrix-method
79		setMethod(
80		f = "index_routledge3",
81		signature = signature(x = "matrix"),
82		definition = function(x) {
83	2x	I <- index_routledge2(x)
84	2x	E <- exp(I)
85	2x	return(E)
86		}
87		)
88
89		# Whittaker ====================================================================
90		#' @export
91		#' @rdname index_whittaker
92		#' @aliases index_whittaker,matrix-method
93		setMethod(
94		f = "index_whittaker",
95		signature = signature(x = "matrix"),
96		definition = function(x) {
97	2x	x <- x > 0 # presence/absence
98
99		# Total number of taxa recorded in the system
100	2x	S <- sum(colSums(x) > 0)
101		# Mean taxa richness
102	2x	alpha <- mean(rowSums(x))
103
104	2x	W <- (S / alpha) - 1
105	2x	W
106		}
107		)
108
109		# Wilson =======================================================================
110		#' @export
111		#' @rdname index_wilson
112		#' @aliases index_wilson,matrix-method
113		setMethod(
114		f = "index_wilson",
115		signature = signature(x = "matrix"),
116		definition = function(x) {
117	2x	x <- x > 0 # presence/absence
118
119	2x	m <- nrow(x)
120	2x	first_transect <- x[1, ]
121	2x	last_transect <- x[m, ]
122
123		# Mean taxa richness
124	2x	alpha <- mean(rowSums(x))
125		# Number of taxa gained along the transect
126	2x	absent <- x[, !first_transect] # Keep only taxa absent in the first transect
127	2x	gained <- sum(apply(X = absent, MARGIN = 2, FUN = cumsum)[m, ] > 0)
128		# Number of taxa lost
129	2x	lost <- sum(!last_transect)
130
131	2x	beta <- (gained + lost) / (2 * alpha)
132	2x	beta
133		}
134		)

1		# PLOT DIVERSITY
2		#' @include AllGenerics.R
3		NULL
4
5		# DiversityIndex ===============================================================
6		#' @export
7		#' @method plot DiversityIndex
8		plot.DiversityIndex <- function(x, log = "x",
9		col.mean = "#DDAA33", col.interval = "#004488",
10		lty.mean = "solid", lty.interval = "dashed",
11		lwd.mean = 1, lwd.interval = 1,
12		xlab = NULL, ylab = NULL,
13		main = NULL, sub = NULL,
14		ann = graphics::par("ann"),
15		axes = TRUE, frame.plot = axes,
16		panel.first = NULL, panel.last = NULL, ...) {
17		## Prepare data
18	2x	count <- as.data.frame(x)
19
20		## Open new window
21	2x	grDevices::dev.hold()
22	2x	on.exit(grDevices::dev.flush(), add = TRUE)
23	2x	graphics::plot.new()
24
25		## Set plotting coordinates
26	2x	xlim <- range(count$size)
27	2x	ylim <- range(count$index)
28	2x	if (length(x@simulation) != 0) {
29	2x	xlim <- range(x@simulation[, "size"])
30	2x	ylim <- range(x@simulation[, c("lower", "upper")])
31		}
32	2x	graphics::plot.window(xlim = xlim, ylim = ylim, log = log)
33
34		## Evaluate pre-plot expressions
35	2x	panel.first
36
37		## Plot
38	2x	graphics::points(x = count$size, y = count$index, ...)
39
40		## Simulated assemblages
41	2x	if (length(x@simulation) != 0) {
42	2x	refined <- x@simulation
43	2x	graphics::lines(x = refined[, "size"], y = refined[, "mean"],
44	2x	col = col.mean, lty = lty.mean, lwd = lwd.mean)
45	2x	graphics::lines(x = refined[, "size"], y = refined[, "lower"],
46	2x	col = col.interval, lty = lty.interval, lwd = lwd.interval)
47	2x	graphics::lines(x = refined[, "size"], y = refined[, "upper"],
48	2x	col = col.interval, lty = lty.interval, lwd = lwd.interval)
49		}
50
51		## Evaluate post-plot and pre-axis expressions
52	2x	panel.last
53
54		## Construct axis
55	2x	if (axes) {
56	2x	graphics::axis(side = 1, las = 1)
57	2x	graphics::axis(side = 2, las = 1)
58		}
59
60		## Plot frame
61	2x	if (frame.plot) {
62	2x	graphics::box()
63		}
64
65		## Add annotation
66	2x	if (ann) {
67	2x	y_lab <- switch (
68	2x	class(x),
69	2x	HeterogeneityIndex = tr_("Heterogeneity"),
70	2x	EvennessIndex = tr_("Evenness"),
71	2x	RichnessIndex = tr_("Richness"),
72	2x	tr_("Diversity")
73		)
74	2x	xlab <- xlab %\|\|% tr_("Sample size")
75	2x	ylab <- ylab %\|\|% sprintf("%s (%s)", y_lab, x@method)
76	2x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab)
77		}
78
79	2x	invisible(x)
80		}
81
82		#' @export
83		#' @rdname plot_diversity
84		#' @aliases plot,DiversityIndex,missing-method
85		setMethod("plot", c(x = "DiversityIndex", y = "missing"), plot.DiversityIndex)
86
87		# SHE ==========================================================================
88		#' @export
89		#' @rdname she
90		#' @aliases she,matrix-method
91		setMethod(
92		f = "she",
93		signature = c(object = "matrix"),
94		definition = function(object, unbiased = FALSE,
95		xlab = NULL, ylab = NULL,
96		main = NULL, sub = NULL,
97		ann = graphics::par("ann"),
98		axes = TRUE, frame.plot = axes,
99		panel.first = NULL, panel.last = NULL,
100		legend = list(x = "right"), ...) {
101		## Prepare data
102	1x	data <- .she(object, unbiased = unbiased)
103
104		## Open new window
105	1x	grDevices::dev.hold()
106	1x	on.exit(grDevices::dev.flush(), add = TRUE)
107	1x	graphics::plot.new()
108
109		## Set plotting coordinates
110	1x	xlim <- range(data[, 4])
111	1x	ylim <- range(data[, -4])
112	1x	graphics::plot.window(xlim = xlim, ylim = ylim)
113
114		## Evaluate pre-plot expressions
115	1x	panel.first
116
117		## Plot
118	1x	col <- c("black", "red", "blue")
119	1x	for (i in c(1, 2, 3)) {
120	3x	graphics::lines(x = data[, 4], y = data[, i], col = col[i], lty = i)
121	3x	graphics::points(x = data[, 4], y = data[, i], col = col[i], pch = 16)
122		}
123
124		## Evaluate post-plot and pre-axis expressions
125	1x	panel.last
126
127		## Construct axis
128	1x	if (axes) {
129	1x	graphics::axis(side = 1, las = 1)
130	1x	graphics::axis(side = 2, las = 1)
131		}
132
133		## Plot frame
134	1x	if (frame.plot) {
135	1x	graphics::box()
136		}
137
138		## Add annotation
139	1x	if (ann) {
140	1x	xlab <- xlab %\|\|% tr_("Sample size")
141	1x	ylab <- ylab %\|\|% tr_("Diversity")
142	1x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
143		}
144
145		## Legend
146	1x	if (is.list(legend) && length(legend) > 0) {
147	1x	args <- list(legend = c("ln(S)", "H", "ln(E)"),
148	1x	col = col, pch = 16,
149	1x	lty = c(1, 2, 3), bty = "n")
150	1x	args <- utils::modifyList(args, legend)
151	1x	do.call(graphics::legend, args = args)
152		}
153
154	1x	invisible(object)
155		}
156		)
157
158		#' @export
159		#' @rdname she
160		#' @aliases she,data.frame-method
161		setMethod(
162		f = "she",
163		signature = c(object = "data.frame"),
164		definition = function(object, unbiased = FALSE,
165		xlab = NULL, ylab = NULL,
166		main = NULL, sub = NULL,
167		ann = graphics::par("ann"),
168		axes = TRUE, frame.plot = axes,
169		panel.first = NULL, panel.last = NULL,
170		legend = list(x = "right"), ...) {
171	1x	object <- data.matrix(object)
172	1x	methods::callGeneric(object, unbiased = unbiased,
173	1x	xlab = xlab, ylab = ylab,
174	1x	main = main, sub = sub, ann = ann,
175	1x	axes = axes, frame.plot = frame.plot,
176	1x	panel.first = panel.first, panel.last = panel.last,
177	1x	legend = legend, ...)
178		}
179		)
180
181		#' @param object A matrix.
182		#' @noRd
183		.she <- function(object, unbiased = FALSE, ...) {
184
185	1x	n <- nrow(object)
186	1x	m <- ncol(object)
187
188	1x	SHE <- matrix(data = 0, nrow = n, ncol = 4)
189	1x	y <- numeric(m)
190
191	1x	for (i in seq_len(n)) {
192	5x	x <- object[i, ] + y
193	5x	n <- sum(x)
194
195		## Log species abundance
196	5x	S <- log(sum(x > 0)) # Observed species
197
198		## Shannon index
199	5x	H <- index_shannon(x, evenness = FALSE, unbiased = unbiased)
200
201		## log evenness
202	5x	E <- H - S
203
204	5x	SHE[i, ] <- c(S, H, E, n)
205	5x	y <- x
206		}
207
208	1x	SHE
209		}
210
211		# Profile ======================================================================
212		#' @export
213		#' @rdname profiles
214		#' @aliases profiles,matrix-method
215		setMethod(
216		f = "profiles",
217		signature = c(object = "matrix"),
218		definition = function(object, alpha = seq(from = 0, to = 4, by = 0.04),
219		color = NULL, symbol = FALSE,
220		xlab = NULL, ylab = NULL,
221		main = NULL, sub = NULL,
222		ann = graphics::par("ann"),
223		axes = TRUE, frame.plot = axes,
224		panel.first = NULL, panel.last = NULL,
225		legend = list(x = "topright"), ...) {
226		## Prepare data
227	1x	alpha <- alpha[alpha != 1]
228	1x	data <- .profiles(object, alpha = alpha)
229	1x	lab <- rownames(object)
230	1x	n <- nrow(object)
231
232		## Graphical parameters
233	1x	lwd <- list(...)$lwd %\|\|% graphics::par("lwd")
234	1x	if (length(lwd) == 1) lwd <- rep(lwd, length.out = n)
235
236	1x	lty <- list(...)$lty %\|\|% graphics::par("lty")
237	1x	if (length(lty) == 1) lty <- rep(lty, length.out = n)
238	!	if (!isFALSE(symbol)) lty <- khroma::palette_line(symbol)(lab)
239
240	1x	col <- list(...)$col %\|\|% graphics::par("col")
241	1x	if (length(col) == 1) col <- rep(col, length.out = n)
242	1x	if (!isFALSE(color)) col <- khroma::palette_color_discrete(color)(lab)
243
244		## Open new window
245	1x	grDevices::dev.hold()
246	1x	on.exit(grDevices::dev.flush(), add = TRUE)
247	1x	graphics::plot.new()
248
249		## Set plotting coordinates
250	1x	xlim <- range(alpha)
251	1x	ylim <- range(data, finite = TRUE)
252	1x	graphics::plot.window(xlim = xlim, ylim = ylim)
253
254		## Evaluate pre-plot expressions
255	1x	panel.first
256
257		## Plot
258	1x	for (i in seq_len(n)) {
259	5x	graphics::lines(x = alpha, y = data[, i], col = col[i],
260	5x	lty = lty[i], lwd = lwd[i])
261		}
262
263		## Evaluate post-plot and pre-axis expressions
264	1x	panel.last
265
266		## Construct axis
267	1x	if (axes) {
268	1x	graphics::axis(side = 1, las = 1)
269	1x	graphics::axis(side = 2, las = 1)
270		}
271
272		## Plot frame
273	1x	if (frame.plot) {
274	1x	graphics::box()
275		}
276
277		## Add annotation
278	1x	if (ann) {
279	1x	xlab <- xlab %\|\|% "alpha"
280	1x	ylab <- ylab %\|\|% tr_("Diversity")
281	1x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
282		}
283
284		## Legend
285	1x	if (is.list(legend) && length(legend) > 0) {
286	1x	args <- list(legend = lab, col = col, lty = lty, lwd = lwd, bty = "n")
287	1x	args <- utils::modifyList(args, legend)
288	1x	do.call(graphics::legend, args = args)
289		}
290
291	1x	invisible(object)
292		}
293		)
294
295		#' @export
296		#' @rdname profiles
297		#' @aliases profiles,data.frame-method
298		setMethod(
299		f = "profiles",
300		signature = c(object = "data.frame"),
301		definition = function(object, alpha = seq(from = 0, to = 4, by = 0.04),
302		color = NULL, symbol = FALSE,
303		xlab = NULL, ylab = NULL,
304		main = NULL, sub = NULL,
305		ann = graphics::par("ann"),
306		axes = TRUE, frame.plot = axes,
307		panel.first = NULL, panel.last = NULL,
308		legend = list(x = "topright"), ...) {
309	1x	object <- data.matrix(object)
310	1x	methods::callGeneric(object, alpha = seq(from = 0, to = 4, by = 0.04),
311	1x	color = color, symbol = symbol,
312	1x	xlab = xlab, ylab = ylab,
313	1x	main = main, sub = sub, ann = ann,
314	1x	axes = axes, frame.plot = frame.plot,
315	1x	panel.first = panel.first, panel.last = panel.last,
316	1x	legend = legend, ...)
317		}
318		)
319
320		#' @param object A matrix.
321		#' @noRd
322		.profiles <- function(object, alpha = seq(from = 0, to = 4, by = 0.04), ...) {
323
324	1x	n <- nrow(object)
325	1x	m <- length(alpha)
326
327	1x	prof <- matrix(data = 0, nrow = m, ncol = n)
328
329	1x	index_renyi <- function(z, na.rm = FALSE) {
330	5x	z <- z[z > 0] # Remove unobserved species
331	!	if (na.rm) z <- stats::na.omit(z) # Remove NAs
332
333	5x	function(x) {
334	500x	p <- z / sum(z)
335	500x	exp(log(sum(p^x)) / (1 - x))
336		}
337		}
338
339	1x	for (i in seq_len(n)) {
340	5x	x <- object[i, ]
341	5x	f <- index_renyi(x)
342	5x	prof[, i] <- vapply(X = alpha, FUN = f, FUN.VALUE = numeric(1), ...)
343		}
344
345	1x	prof
346		}

1		# SERIOGRAPH
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname seriograph
7		#' @aliases eppm,matrix-method
8		setMethod(
9		f = "eppm",
10		signature = signature(object = "matrix"),
11		definition = function(object) {
12	4x	eppm <- (object - expected(object)) * 100 / rowSums(object)
13	4x	eppm[eppm < 0] <- 0
14	4x	dimnames(eppm) <- dimnames(object)
15	4x	eppm
16		}
17		)
18
19		#' @export
20		#' @rdname seriograph
21		#' @aliases eppm,data.frame-method
22		setMethod(
23		f = "eppm",
24		signature = signature(object = "data.frame"),
25		definition = function(object) {
26	1x	object <- data.matrix(object)
27	1x	methods::callGeneric(object)
28		}
29		)
30
31		#' @export
32		#' @rdname seriograph
33		#' @aliases seriograph,matrix-method
34		setMethod(
35		f = "seriograph",
36		signature = signature(object = "matrix"),
37		definition = function(object, weights = FALSE,
38		fill = "darkgrey", border = NA,
39		axes = TRUE, ...) {
40	1x	plot_ford(object, weights = weights, EPPM = TRUE,
41	1x	fill = fill, border = border, axes = axes)
42	1x	invisible(object)
43		}
44		)
45
46		#' @export
47		#' @rdname seriograph
48		#' @aliases seriograph,data.frame-method
49		setMethod(
50		f = "seriograph",
51		signature = signature(object = "data.frame"),
52		definition = function(object, weights = FALSE,
53		fill = "darkgrey", border = NA,
54		axes = TRUE, ...) {
55	1x	object <- data.matrix(object)
56	1x	methods::callGeneric(object, weights = weights,
57	1x	fill = fill, border = border, axes = axes)
58		}
59		)

1		# ALPHA DIVERSITY
2		#' @include AllGenerics.R
3		NULL
4
5		# Observed =====================================================================
6		#' @export
7		#' @rdname observed
8		#' @aliases observed,numeric-method
9		setMethod(
10		f = "observed",
11		signature = c(x = "numeric"),
12		definition = function(x, na.rm = FALSE, ...) {
13	10x	sum(x > 0, na.rm = na.rm)
14		}
15		)
16
17		nobserved <- function(x, n, na.rm = FALSE) {
18	20x	sum(x == n, na.rm = na.rm)
19		}
20
21		#' @export
22		#' @rdname observed
23		#' @aliases singleton,numeric-method
24		setMethod(
25		f = "singleton",
26		signature = c(x = "numeric"),
27		definition = function(x, na.rm = FALSE, ...) {
28	10x	nobserved(x, n = 1, na.rm = na.rm)
29		}
30		)
31
32		#' @export
33		#' @rdname observed
34		#' @aliases doubleton,numeric-method
35		setMethod(
36		f = "doubleton",
37		signature = c(x = "numeric"),
38		definition = function(x, na.rm = FALSE, ...) {
39	10x	nobserved(x, n = 2, na.rm = na.rm)
40		}
41		)
42
43		# ACE ==========================================================================
44		#' @export
45		#' @rdname index_ace
46		#' @aliases index_ace,numeric-method
47		setMethod(
48		f = "index_ace",
49		signature = c(x = "numeric"),
50		definition = function(x, k = 10, na.rm = FALSE, ...) {
51		## Validation
52	3x	x <- x[x > 0] # Remove unobserved species
53	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
54	!	if (anyNA(x)) return(NA)
55
56	3x	S_rare <- sum(x <= k) # Number of rare species
57	3x	S_abund <- sum(x > k) # Number of abundant species
58	3x	N_rare <- sum(x[x <= k]) # Number of individuals in the rare species
59
60	3x	F1 <- sum(x == 1) # Number of singleton species
61	3x	if (F1 == N_rare) {
62	!	if (getOption("tabula.verbose")) {
63	!	message(
64	!	tr_("ACE is undefined when all rare species are singletons."), "\n",
65	!	tr_("Consider using the bias-corrected Chao1 estimator instead.")
66		)
67		}
68	!	return(NA_real_)
69		}
70
71		## Sample coverage estimate for rare species
72		## ie. proportion of all individuals in rare species that are not singletons
73	3x	C_rare <- 1 - (F1 / N_rare)
74
75		# Coefficient of variation
76	3x	a <- sum(vapply(
77	3x	X = seq_len(k),
78	30x	FUN = function(i, x) { i * (i - 1) * sum(x == i) },
79	3x	FUN.VALUE = double(1), x = x)
80		)
81	3x	g2 <- max((S_rare / C_rare) * (a / (N_rare * (N_rare - 1))) - 1, 0)
82
83	3x	D <- S_abund + (S_rare / C_rare) + (F1 / C_rare) * g2
84	3x	D
85		}
86		)
87
88		# Berger-Parker ================================================================
89		#' @export
90		#' @rdname index_berger
91		#' @aliases index_berger,numeric-method
92		setMethod(
93		f = "index_berger",
94		signature = c(x = "numeric"),
95		definition = function(x, na.rm = FALSE, ...) {
96		## Validation
97	7x	x <- x[x > 0] # Remove unobserved species
98	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
99	!	if (anyNA(x)) return(NA)
100
101	7x	Nmax <- max(x)
102	7x	N <- sum(x)
103	7x	d <- Nmax / N
104	7x	d
105		}
106		)
107
108		# Boone ========================================================================
109		#' @export
110		#' @rdname index_boone
111		#' @aliases index_boone,matrix-method
112		setMethod(
113		f = "index_boone",
114		signature = c(x = "matrix"),
115		definition = function(x, j = NULL, na.rm = FALSE, ...) {
116		## Validation
117	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
118	!	if (anyNA(x)) return(NA)
119
120	!	Y <- colSums(x) # Site-wide totals of each artifact class
121	!	if (is.null(j)) j <- which.max(Y)
122
123	!	W <- Y[j] / Y # Weighting factor
124	!	P <- W * Y / sum(W * Y)
125
126	!	W <- matrix(W, nrow = nrow(x), ncol = ncol(x), byrow = TRUE)
127	!	P <- matrix(P, nrow = nrow(x), ncol = ncol(x), byrow = TRUE)
128
129	!	Wx <- W * x # Weighted counts
130	!	px <- Wx / rowSums(Wx) # Weighted percentages
131
132	!	rowSums((px - P)^2)
133		}
134		)
135
136		# Brillouin ====================================================================
137		#' @export
138		#' @rdname index_brillouin
139		#' @aliases index_brillouin,numeric-method
140		setMethod(
141		f = "index_brillouin",
142		signature = c(x = "numeric"),
143		definition = function(x, evenness = FALSE, na.rm = FALSE, ...) {
144		## Validation
145	13x	x <- x[x > 0] # Remove unobserved species
146	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
147	!	if (anyNA(x)) return(NA)
148
149	13x	N <- sum(x)
150	13x	bri <- (lfactorial(N) - sum(lfactorial(x))) / N
151
152	13x	if (evenness) {
153	6x	N <- sum(x)
154	6x	S <- length(x) # richness = number of different species
155	6x	a <- trunc(N / S)
156	6x	r <- N - S * a
157	6x	c <- (S - r) * lfactorial(a) + r * lfactorial(a + 1)
158	6x	HBmax <- (1 / N) * (lfactorial(N) - c)
159	6x	bri <- bri / HBmax
160		}
161
162	13x	bri
163		}
164		)
165
166		# Chao1 ========================================================================
167		#' @export
168		#' @rdname index_chao1
169		#' @aliases index_chao1,numeric-method
170		setMethod(
171		f = "index_chao1",
172		signature = c(x = "numeric"),
173		definition = function(x, unbiased = FALSE, improved = FALSE, na.rm = FALSE, ...) {
174		## Validation
175	8x	x <- x[x > 0] # Remove unobserved species
176	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
177	!	if (anyNA(x)) return(NA)
178
179	8x	S <- length(x) # Number of observed species
180	8x	N <- sum(x) # Number of individuals
181	8x	f1 <- sum(x == 1) # Number of singleton species
182	8x	f2 <- sum(x == 2) # Number of doubleton species
183
184	8x	if (unbiased) {
185	4x	D <- S + (((N - 1) / N) * f1 * (f1 - 1)) / (2 * (f2 + 1))
186		} else {
187	4x	if (f2 == 0) {
188	!	D <- S + ((N - 1) / N) * f1 * ((f1 - 1) / 2)
189		} else {
190	4x	D <- S + ((N - 1) / N) * (f1^2 / (2 * f2))
191		}
192		}
193	8x	if (improved) {
194	3x	f3 <- sum(x == 3) # Number of triple species
195	3x	f4 <- sum(x == 4) # Number of quadruple species
196	3x	if (f4 == 0) {
197	1x	if (getOption("tabula.verbose")) {
198	!	message(
199	!	tr_("Improved Chao1 estimator is undefined when there is no quadruple species.")
200		)
201		}
202	1x	return(NA_real_)
203		}
204
205	2x	k <- f1 - ((N - 3) / (N - 1)) * ((f2 * f3) / (2 * f4))
206	2x	D <- D + ((N - 3) / N) * (f3 / (4 * f4)) * max(k, 0)
207		}
208
209	7x	D
210		}
211		)
212
213		# Chao2 ========================================================================
214		#' @export
215		#' @rdname index_chao2
216		#' @aliases index_chao2,matrix-method
217		setMethod(
218		f = "index_chao2",
219		signature = c(x = "matrix"),
220		definition = function(x, unbiased = FALSE, improved = FALSE, ...) {
221	6x	x <- x > 0 # Convert to incidence
222
223	6x	q <- colSums(x) # Number of species in the assemblage
224	6x	q <- q[q > 0] # Remove unobserved species
225	6x	S <- length(q) # Number of observed species
226	6x	t <- nrow(x) # Total number of sampling units
227	6x	q1 <- sum(q == 1) # Number of unique species in the assemblage
228	6x	q2 <- sum(q == 2) # Number of duplicate species in the assemblage
229
230	6x	if (unbiased) {
231	2x	D <- S + ((t - 1) / t) * q1 * ((q1 - 1) / (2 * (q2 + 1)))
232		} else {
233	4x	if (q2 == 0) {
234	1x	D <- S + ((t - 1) / t) * q1 * ((q1 - 1) / 2)
235		} else {
236	3x	D <- S + ((t - 1) / t) * (q1^2 / (2 * q2))
237		}
238		}
239	6x	if (improved) {
240	2x	q3 <- sum(q == 3) # Number of triple species
241	2x	q4 <- sum(q == 4) # Number of quadruple species
242	2x	if (q4 == 0) {
243	!	if (getOption("tabula.verbose")) {
244	!	message(
245	!	tr_("Improved Chao2 estimator is undefined when there is no quadruple species.")
246		)
247		}
248	!	return(NA_real_)
249		}
250
251	2x	k <- q1 - ((t - 3) / (t - 1)) * ((q2 * q3) / (2 * q4))
252	2x	D <- D + ((t - 3) / (4 * t)) * (q3 / q4) * max(k, 0)
253		}
254
255	6x	D
256		}
257		)
258
259		# ICE ==========================================================================
260		#' @export
261		#' @rdname index_ice
262		#' @aliases index_ice,matrix-method
263		setMethod(
264		f = "index_ice",
265		signature = c(x = "matrix"),
266		definition = function(x, k = 10, ...) {
267	2x	x <- x > 0 # Convert to incidence
268
269	2x	q <- colSums(x) # Number of species in the assemblage
270	2x	q <- q[q > 0] # Remove unobserved species
271
272	2x	S_infr <- sum(q <= k) # Number of infrequent species
273	2x	S_freq <- sum(q > k) # Number of frequent species
274	2x	N_infr <- sum(q[q <= k]) # Number of incidences in the infrequent species
275		# Number of sampling units that include at least one infrequent species
276	2x	t <- sum(rowSums(x[, q <= k]) != 0)
277
278	2x	q1 <- sum(q == 1) # Number of unique species in the assemblage
279	2x	if (q1 == N_infr) {
280	!	if (getOption("tabula.verbose")) {
281	!	message(
282	!	tr_("ICE is undefined when all rare species are singletons."), "\n",
283	!	tr_("Consider using the bias-corrected Chao2 estimator instead.")
284		)
285		}
286	!	return(NA_real_)
287		}
288
289		## Sample coverage estimate
290		## ie. proportion of all incidences of infrequent species that are not uniques
291	2x	C_infr <- 1 - (q1 / N_infr)
292
293		## Coefficient of variation
294	2x	a <- sum(vapply(
295	2x	X = seq_len(k),
296	20x	FUN = function(x, q) { x * (x - 1) * sum(q == x) },
297	2x	FUN.VALUE = double(1), q = q)
298		)
299	2x	b <- sum(vapply(
300	2x	X = seq_len(k),
301	20x	FUN = function(x, q) { x * sum(q == x) },
302	2x	FUN.VALUE = double(1), q = q)
303		)
304	2x	c <- sum(vapply(
305	2x	X = seq_len(k),
306	20x	FUN = function(x, q) { x * sum(q == x) - 1 },
307	2x	FUN.VALUE = double(1), q = q)
308		)
309	2x	g2 <- max((S_infr / C_infr) * (t / (t - 1)) * (a / (b * c)) - 1, 0)
310
311	2x	D <- S_freq + S_infr / C_infr + q1 * g2 / C_infr
312	2x	D
313		}
314		)
315
316		# Margalef =====================================================================
317		#' @export
318		#' @rdname index_margalef
319		#' @aliases index_margalef,numeric-method
320		setMethod(
321		f = "index_margalef",
322		signature = c(x = "numeric"),
323		definition = function(x, na.rm = FALSE, ...) {
324		## Validation
325	7x	x <- x[x > 0] # Remove unobserved species
326	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
327	!	if (anyNA(x)) return(NA)
328
329	7x	N <- sum(x) # Number of individuals
330	7x	S <- length(x) # Number of observed species
331	7x	D <- (S - 1) / log(N, base = exp(1))
332	7x	D
333		}
334		)
335
336		# McIntosh =====================================================================
337		#' @export
338		#' @rdname index_mcintosh
339		#' @aliases index_mcintosh,numeric-method
340		setMethod(
341		f = "index_mcintosh",
342		signature = c(x = "numeric"),
343		definition = function(x, evenness = FALSE, na.rm = FALSE, ...) {
344		## Validation
345	12x	x <- x[x > 0] # Remove unobserved species
346	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
347	!	if (anyNA(x)) return(NA)
348
349	12x	N <- sum(x)
350	12x	S <- length(x) # richness = number of different species
351	12x	U <- sqrt(sum(x^2))
352
353	12x	if (evenness) {
354	6x	mac <- (N - U) / (N - (N / sqrt(S)))
355		} else {
356	6x	mac <- (N - U) / (N - sqrt(N))
357		}
358	12x	mac
359		}
360		)
361
362		# Menhinick ====================================================================
363		#' @export
364		#' @rdname index_menhinick
365		#' @aliases index_menhinick,numeric-method
366		setMethod(
367		f = "index_menhinick",
368		signature = c(x = "numeric"),
369		definition = function(x, na.rm = FALSE, ...) {
370		## Validation
371	7x	x <- x[x > 0] # Remove unobserved species
372	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
373	!	if (anyNA(x)) return(NA)
374
375	7x	N <- sum(x) # Number of individuals
376	7x	S <- length(x) # Number of observed species
377	7x	D <- S / sqrt(N)
378	7x	D
379		}
380		)
381
382		# Shannon ======================================================================
383		#' @export
384		#' @rdname index_shannon
385		#' @aliases index_shannon,numeric-method
386		setMethod(
387		f = "index_shannon",
388		signature = c(x = "numeric"),
389		definition = function(x, evenness = FALSE, unbiased = FALSE, ACE = FALSE,
390		base = exp(1), na.rm = FALSE, ...) {
391		## Validation
392	159787x	x <- x[x > 0] # Remove unobserved species
393	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
394	!	if (anyNA(x)) return(NA)
395
396	159787x	N <- sum(x)
397	159787x	S <- length(x) # richness = number of different species
398	159787x	p <- x / N
399	159787x	Hmax <- log(p, base)
400	159787x	Hmax[is.infinite(Hmax)] <- 0
401
402	159787x	H <- -sum(p * Hmax)
403	159787x	if (unbiased) {
404	!	if (ACE) S <- index_ace(x, ...)
405	1x	H <- H + (S - 1) / (2 * N)
406		}
407	386x	if (evenness) H <- H / log(S)
408	159787x	H
409		}
410		)
411
412		# Simpson ======================================================================
413		#' @export
414		#' @rdname index_simpson
415		#' @aliases index_simpson,numeric-method
416		setMethod(
417		f = "index_simpson",
418		signature = c(x = "numeric"),
419		definition = function(x, evenness = FALSE, unbiased = FALSE, na.rm = FALSE, ...) {
420		## Validation
421	16x	x <- x[x > 0] # Remove unobserved species
422	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
423	!	if (anyNA(x)) return(NA)
424
425	16x	N <- sum(x)
426	2x	if (unbiased) D <- sum(x * (x - 1)) / (N * (N - 1))
427	14x	else D <- sum((x / N)^2)
428
429	16x	if (evenness) {
430	5x	D <- 1 / D
431	5x	S <- length(x)
432	5x	D <- D / S
433		}
434	16x	D
435		}
436		)
437
438		# Squares ======================================================================
439		#' @export
440		#' @rdname index_squares
441		#' @aliases index_squares,numeric-method
442		setMethod(
443		f = "index_squares",
444		signature = c(x = "numeric"),
445		definition = function(x, na.rm = FALSE, ...) {
446		## Validation
447	1x	x <- x[x > 0] # Remove unobserved species
448	!	if (na.rm) x <- stats::na.omit(x) # Remove NAs
449	!	if (anyNA(x)) return(NA)
450
451	1x	S <- length(x) # Number of observed species
452	1x	N <- sum(x) # Number of individuals
453	1x	f1 <- sum(x == 1) # Number of singleton species
454
455	1x	Ssq <- S + f1^2 * sum(x^2) / (N^2 - f1 * S)
456	1x	Ssq
457		}
458		)

1		# RAREFACTION
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname rarefaction
7		#' @aliases rarefaction,matrix-method
8		setMethod(
9		f = "rarefaction",
10		signature = signature(object = "matrix"),
11		definition = function(object, sample = NULL, method = c("hurlbert", "baxter"),
12		step = 1) {
13		## Validation
14	2x	method <- match.arg(method, several.ok = FALSE)
15
16	2x	n <- nrow(object)
17	2x	if (is.null(sample)) {
18	!	sample <- rowSums(object)
19		}
20	2x	if (length(sample) == 1) {
21	2x	sample <- rep(sample, n)
22		}
23	2x	k <- seq(from = 1, to = max(sample), by = step)
24
25		## Matrix of results
26	2x	z <- matrix(data = NA_real_, nrow = n, ncol = length(k))
27	2x	row_names <- rownames(object) %\|\|% paste0("S", seq_len(n)) # Fix names
28	2x	dimnames(z) <- list(row_names, k)
29
30	2x	for (i in seq_len(n)) {
31	7x	spl <- k[k <= sample[[i]]]
32	7x	rare <- vapply(
33	7x	X = spl,
34	7x	FUN = function(s, x, f) f(x, s),
35	7x	FUN.VALUE = numeric(1),
36	7x	x = object[i, ],
37	7x	f = get_index(method) # Select method
38		)
39	7x	z[i, seq_along(rare)] <- rare
40		}
41
42	2x	.RarefactionIndex(
43	2x	z,
44	2x	labels = row_names,
45	2x	size = as.integer(k),
46	2x	method = method
47		)
48		}
49		)
50
51		#' @export
52		#' @rdname rarefaction
53		#' @aliases rarefaction,data.frame-method
54		setMethod(
55		f = "rarefaction",
56		signature = signature(object = "data.frame"),
57		definition = function(object, sample = NULL, method = c("hurlbert", "baxter"),
58		step = 1) {
59	1x	object <- data.matrix(object)
60	1x	methods::callGeneric(object, sample = sample, method = method, step = step)
61		}
62		)
63
64		# Plot =========================================================================
65		#' @export
66		#' @method plot RarefactionIndex
67		plot.RarefactionIndex <- function(x, color = NULL, symbol = FALSE,
68		xlab = NULL, ylab = NULL,
69		main = NULL, sub = NULL,
70		ann = graphics::par("ann"),
71		axes = TRUE, frame.plot = axes,
72		panel.first = NULL, panel.last = NULL,
73		legend = list(x = "topleft"), ...) {
74		## Prepare data
75	2x	n <- nrow(x)
76	2x	lab <- labels(x)
77
78		## Graphical parameters
79	2x	lwd <- list(...)$lwd %\|\|% graphics::par("lwd")
80	2x	if (length(lwd) == 1) lwd <- rep(lwd, length.out = n)
81
82	2x	lty <- list(...)$lty %\|\|% graphics::par("lty")
83	2x	if (length(lty) == 1) lty <- rep(lty, length.out = n)
84	!	if (!isFALSE(symbol)) lty <- khroma::palette_line(symbol)(lab)
85
86	2x	col <- list(...)$col %\|\|% graphics::par("col")
87	2x	if (length(col) == 1) col <- rep(col, length.out = n)
88	2x	if (!isFALSE(color)) col <- khroma::palette_color_discrete(color)(lab)
89
90		## Open new window
91	2x	grDevices::dev.hold()
92	2x	on.exit(grDevices::dev.flush(), add = TRUE)
93	2x	graphics::plot.new()
94
95		## Set plotting coordinates
96	2x	xlim <- range(x@size)
97	2x	ylim <- range(x, na.rm = TRUE)
98	2x	graphics::plot.window(xlim = xlim, ylim = ylim)
99
100		## Evaluate pre-plot expressions
101	2x	panel.first
102
103		## Plot
104	2x	for (i in seq_len(n)) {
105	10x	graphics::lines(x = x@size, y = x[i, ], col = col[i],
106	10x	lwd = lwd[i], lty = lty[i])
107		}
108
109		## Evaluate post-plot and pre-axis expressions
110	2x	panel.last
111
112		## Construct axis
113	2x	if (axes) {
114	2x	graphics::axis(side = 1, las = 1)
115	2x	graphics::axis(side = 2, las = 1)
116		}
117
118		## Plot frame
119	2x	if (frame.plot) {
120	2x	graphics::box()
121		}
122
123		## Add annotation
124	2x	if (ann) {
125	2x	xlab <- xlab %\|\|% tr_("Sample size")
126	2x	ylab <- ylab %\|\|% tr_("Expected species index")
127	2x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab)
128		}
129
130		## Legend
131	2x	if (is.list(legend) && length(legend) > 0) {
132	1x	args <- list(legend = lab, col = col, lty = lty, lwd = lwd, bty = "n")
133	1x	args <- utils::modifyList(args, legend)
134	1x	do.call(graphics::legend, args = args)
135		}
136
137	2x	invisible(x)
138		}
139
140		#' @export
141		#' @rdname plot_rarefaction
142		#' @aliases plot,RarefactionIndex,missing-method
143		setMethod("plot", c(x = "RarefactionIndex", y = "missing"), plot.RarefactionIndex)
144
145		# Index ========================================================================
146		#' @export
147		#' @rdname index_baxter
148		#' @aliases index_baxter,numeric-method
149		setMethod(
150		f = "index_baxter",
151		signature = signature(x = "numeric"),
152		definition = function(x, sample, ...) {
153	116x	x <- x[x > 0]
154	116x	N <- sum(x)
155
156	116x	E <- suppressWarnings(exp(ramanujan(N - x) + ramanujan(N - sample) -
157	116x	ramanujan(N - x - sample) - ramanujan(N)))
158	116x	sum(1 - E, na.rm = FALSE)
159		}
160		)
161
162		#' @export
163		#' @rdname index_hurlbert
164		#' @aliases index_hurlbert,numeric-method
165		setMethod(
166		f = "index_hurlbert",
167		signature = signature(x = "numeric"),
168		definition = function(x, sample, ...) {
169	27x	x <- x[x > 0]
170	27x	N <- sum(x)
171
172	27x	E <- vapply(
173	27x	X = x,
174	27x	FUN = function(x, N, sample) {
175	204x	if (N - x > sample) {
176	185x	combination(N - x, sample) / combination(N, sample)
177		} else {
178	19x	0
179		}
180		},
181	27x	FUN.VALUE = double(1),
182	27x	N, sample
183		)
184	27x	sum(1 - E, na.rm = FALSE)
185		}
186		)

1		# MUTATORS
2		#' @include AllClasses.R
3		NULL
4
5		# Getters ======================================================================
6		#' @export
7		#' @method labels DiversityIndex
8	2x	labels.DiversityIndex <- function(object, ...) object@labels
9
10		#' @rdname mutators
11		#' @aliases labels,DiversityIndex-method
12		setMethod("labels", "DiversityIndex", labels.DiversityIndex)
13
14		#' @export
15		#' @method labels RarefactionIndex
16	2x	labels.RarefactionIndex <- function(object, ...) object@labels
17
18		#' @rdname mutators
19		#' @aliases labels,RarefactionIndex-method
20		setMethod("labels", "RarefactionIndex", labels.RarefactionIndex)
21
22		#' @export
23		#' @rdname mutators
24		#' @aliases get_method,DiversityIndex-method
25		setMethod(
26		f = "get_method",
27		signature = "DiversityIndex",
28	14x	definition = function(x) x@method
29		)
30
31		# Setters ======================================================================

1		# HELPERS
2
3		## https://michaelchirico.github.io/potools/articles/developers.html
4		tr_ <- function(...) {
5	16x	enc2utf8(gettext(paste0(...), domain = "R-tabula"))
6		}
7
8		is_incidence <- function(x) {
9	1x	if (is.logical(x)) return(TRUE)
10	20x	x <- as.numeric(x)
11	20x	all(x == 0 \| x == 1, na.rm = TRUE)
12		}
13
14		#' Rolling Sum
15		#'
16		#' @param x A [`numeric`] vector.
17		#' @param n An [`integer`] giving the rolling window size.
18		#' @return A [`numeric`] vector.
19		#' @keywords internal
20		#' @noRd
21		roll_sum <- function(x, n = 2) {
22	3x	utils::tail(cumsum(x) - cumsum(c(rep(0, n), utils::head(x, -n))), -n + 1)
23		}
24
25		#' Plotting Dimensions of Character Strings
26		#'
27		#' Convert string length in inch to number of (margin) lines.
28		#' @param x A [`character`] vector of string whose length is to be calculated.
29		#' @param ... Further parameter to be passed to [graphics::strwidth()]`, such as
30		#' `cex`.
31		#' @return
32		#' A [`numeric`] vector (maximum string width in units of margin lines).
33		#' @note For internal use only.
34		#' @family graphic tools
35		#' @keywords internal
36		#' @noRd
37		inch2line <- function(x, ...) {
38	86x	(max(graphics::strwidth(x, units = "inch", ...)) /
39	86x	graphics::par("cin")[2] + graphics::par("mgp")[2]) * graphics::par("cex")
40		}
41
42		#' Circle
43		#'
44		#' Draws a circle.
45		#' @param x,y A length-one [`numeric`] vector giving the coordinates of the
46		#' center of the circle.
47		#' @param radius A length-one [`numeric`] vector giving the radius of the
48		#' circle.
49		#' @param n A length-on [`integer`] vector specifying the number of vertices to
50		#' draw the circle.
51		#' @param ... Further parameters to be passed to [graphics::polygon()].
52		#' @return
53		#' `circle()` is called it for its side-effects: it results in a graphic
54		#' being displayed.
55		#' @example inst/examples/ex-circle.R
56		#' @author N. Frerebeau
57		#' @keywords internal
58		#' @noRd
59		circle <- function(x, y, radius, ..., n = 100) {
60	308x	angle.inc <- 2 * pi / n
61	308x	angles <- seq(0, 2 * pi - angle.inc, by = angle.inc)
62
63	308x	xv <- cos(angles) * radius + x
64	308x	yv <- sin(angles) * radius + y
65	308x	graphics::polygon(xv, yv, ...)
66		}
67
68		#' Label Percentages
69		#'
70		#' @param x A [`numeric`] vector.
71		#' @param digits An [`integer`] indicating the number of decimal places.
72		#' If `NULL` (the default), breaks will have the minimum number of digits
73		#' needed to show the difference between adjacent values.
74		#' @param trim A [`logical`] scalar. If `FALSE` (the default), values are
75		#' right-justified to a common width (see [base::format()]).
76		#' @return A [`character`] vector.
77		#' @keywords internal
78		#' @noRd
79		label_percent <- function(x, digits = NULL, trim = FALSE) {
80	7x	i <- !is.na(x)
81	7x	y <- x[i]
82	7x	y <- abs(y) * 100
83	7x	y <- format(y, trim = trim, digits = digits)
84	7x	y <- paste0(y, "%")
85	7x	x[i] <- y
86	7x	x
87		}

1		# SIMILARITY
2		#' @include AllGenerics.R
3		NULL
4
5		# Similarity ===================================================================
6		#' @export
7		#' @rdname similarity
8		#' @aliases similarity,matrix-method
9		setMethod(
10		f = "similarity",
11		signature = c(object = "matrix"),
12		definition = function(object, method = c("brainerd", "bray", "jaccard",
13		"morisita", "sorenson")) {
14		## Validation
15	1x	if (any(method == "binomial")) {
16	!	stop("Use occurrence(x, method = \"bionomial\") instead.", call. = FALSE)
17		}
18	1x	method <- match.arg(method, several.ok = FALSE)
19		## Select method
20	1x	fun <- get_index(method)
21
22		## Pairwise comparison
23		## Sample/case comparisons
24	1x	m <- nrow(object)
25	1x	labels <- rownames(object)
26	1x	beta <- function(x, f) f(object[x[1], ], object[x[2], ])
27
28	1x	cbn <- utils::combn(seq_len(m), m = 2, simplify = TRUE)
29	1x	index <- apply(X = cbn, MARGIN = 2, FUN = beta, f = fun)
30
31		## Matrix of results
32	1x	diag_value <- ifelse(method == "brainerd", 200, 1)
33	1x	mtx <- matrix(data = diag_value, nrow = m, ncol = m,
34	1x	dimnames = list(labels, labels))
35	1x	mtx[lower.tri(mtx, diag = FALSE)] <- index
36	1x	mtx <- t(mtx)
37	1x	mtx[lower.tri(mtx, diag = FALSE)] <- index
38
39	1x	sim <- stats::as.dist(mtx)
40	1x	attr(sim, "method") <- method
41	1x	sim
42		}
43		)
44
45		#' @export
46		#' @rdname similarity
47		#' @aliases similarity,data.frame-method
48		setMethod(
49		f = "similarity",
50		signature = c(object = "data.frame"),
51		definition = function(object, method = c("brainerd", "bray", "jaccard",
52		"morisita", "sorenson")) {
53	1x	object <- data.matrix(object)
54	1x	methods::callGeneric(object, method = method)
55		}
56		)
57
58		# Index ========================================================================
59		## Qualitative index -----------------------------------------------------------
60		#' @export
61		#' @rdname index_jaccard
62		#' @aliases index_jaccard,character,character-method
63		setMethod(
64		f = "index_jaccard",
65		signature = c(x = "character", y = "character"),
66		definition = function(x, y) {
67	1x	inter <- length(intersect(x, y))
68	1x	union <- length(x) + length(y) - inter
69	1x	inter / union
70		}
71		)
72
73		#' @export
74		#' @rdname index_jaccard
75		#' @aliases index_jaccard,logical,logical-method
76		setMethod(
77		f = "index_jaccard",
78		signature = c(x = "logical", y = "logical"),
79		definition = function(x, y) {
80		## Validation
81	1x	arkhe::assert_length(y, length(x))
82
83	1x	a <- sum(x)
84	1x	b <- sum(y)
85	1x	j <- sum((x + y) == 2)
86	1x	Cj <- j / (a + b - j)
87	1x	Cj
88		}
89		)
90
91		#' @export
92		#' @rdname index_jaccard
93		#' @aliases index_jaccard,numeric,numeric-method
94		setMethod(
95		f = "index_jaccard",
96		signature = c(x = "numeric", y = "numeric"),
97		definition = function(x, y) {
98		## presence/absence
99	1x	x <- x > 0
100	1x	y <- y > 0
101	1x	methods::callGeneric(x, y)
102		}
103		)
104
105		#' @export
106		#' @rdname index_sorenson
107		#' @aliases index_sorenson,logical,logical-method
108		setMethod(
109		f = "index_sorenson",
110		signature = c(x = "logical", y = "logical"),
111		definition = function(x, y) {
112		## Validation
113	1x	arkhe::assert_length(y, length(x))
114
115	1x	a <- sum(x)
116	1x	b <- sum(y)
117	1x	j <- sum((x + y) == 2)
118	1x	Cs <- 2 * j / (a + b)
119	1x	Cs
120		}
121		)
122
123		#' @export
124		#' @rdname index_sorenson
125		#' @aliases index_sorenson,numeric,numeric-method
126		setMethod(
127		f = "index_sorenson",
128		signature = c(x = "numeric", y = "numeric"),
129		definition = function(x, y) {
130		## presence/absence
131	1x	x <- x > 0
132	1x	y <- y > 0
133	1x	methods::callGeneric(x, y)
134		}
135		)
136
137		## Quantitative index ----------------------------------------------------------
138		#' @export
139		#' @rdname index_brainerd
140		#' @aliases index_brainerd,numeric,numeric-method
141		setMethod(
142		f = "index_brainerd",
143		signature = c(x = "numeric", y = "numeric"),
144		definition = function(x, y) {
145		## Validation
146	37x	arkhe::assert_length(y, length(x))
147
148	37x	a <- x / sum(x)
149	37x	b <- y / sum(y)
150	37x	Cb <- 2 - sum(abs(a - b))
151	37x	Cb * 100
152		}
153		)
154
155		#' @export
156		#' @rdname index_bray
157		#' @aliases index_bray,numeric,numeric-method
158		setMethod(
159		f = "index_bray",
160		signature = c(x = "numeric", y = "numeric"),
161		definition = function(x, y) {
162		## Validation
163	1x	arkhe::assert_length(y, length(x))
164
165	1x	a <- sum(x) # Number of individuals in site A
166	1x	b <- sum(y) # Number of individuals in site B
167
168	1x	Cs <- 1 - sum(abs(y - x)) / sum(x + y)
169	1x	Cs
170		}
171		)
172
173		#' @export
174		#' @rdname index_morisita
175		#' @aliases index_morisita,numeric,numeric-method
176		setMethod(
177		f = "index_morisita",
178		signature = c(x = "numeric", y = "numeric"),
179		definition = function(x, y) {
180		## Validation
181	1x	arkhe::assert_length(y, length(x))
182
183	1x	a <- sum(x) # Number of individuals in site A
184	1x	b <- sum(y) # Number of individuals in site B
185	1x	da <- sum(x^2) / a^2
186	1x	db <- sum(y^2) / b^2
187
188	1x	Cm <- (2 * sum(x * y)) / ((da + db) * a * b)
189	1x	Cm
190		}
191		)

1		# DIVERSITY
2		#' @include AllGenerics.R
3		NULL
4
5		# Index ========================================================================
6		index_observed <- function(x, ...) {
7	159390x	sum(x > 0, ...) # Number of observed species
8		}
9		get_index <- function(x) {
10	319170x	match.fun(sprintf("index_%s", x))
11		}
12		do_index <- function(x, method, ...) {
13	319140x	f <- get_index(method)
14	319140x	f(x, ...)
15		}
16
17		#' Compute a Diversity Index
18		#'
19		#' @param x A [`numeric`] [`matrix`].
20		#' @param method A [`character`] string specifying the measure to be computed.
21		#' @param by_row A [`logical`] scalar: should `method` be computed for each row?
22		#' @param ... Further parameters to be passed to `method`.
23		#' @return A [DiversityIndex-class] object.
24		#' @author N. Frerebeau
25		#' @keywords internal
26		#' @noRd
27		index_diversity <- function(x, method, ..., by_row = TRUE) {
28	16x	fun <- get_index(method)
29	16x	if (by_row) {
30	14x	idx <- apply(X = x, MARGIN = 1, FUN = fun, ...)
31		} else {
32	2x	idx <- fun(x, ...)
33		}
34
35		## Fix names
36	16x	row_names <- rownames(x) %\|\|% paste0("S", seq_along(idx))
37
38	16x	.DiversityIndex(
39	16x	idx,
40	16x	labels = row_names,
41	16x	size = as.integer(rowSums(x)),
42	16x	data = x,
43	16x	method = method
44		)
45		}
46
47		# Heterogeneity ================================================================
48		#' @export
49		#' @rdname heterogeneity
50		#' @aliases heterogeneity,matrix-method
51		setMethod(
52		f = "heterogeneity",
53		signature = c(object = "matrix"),
54		definition = function(object, ...,
55		method = c("berger", "boone", "brillouin",
56		"mcintosh", "shannon", "simpson")) {
57	6x	method <- match.arg(method, several.ok = FALSE)
58	6x	by_row <- method != "boone"
59	6x	index <- index_diversity(object, method, ..., evenness = FALSE,
60	6x	by_row = by_row)
61	6x	.HeterogeneityIndex(index)
62		}
63		)
64
65		#' @export
66		#' @rdname heterogeneity
67		#' @aliases heterogeneity,data.frame-method
68		setMethod(
69		f = "heterogeneity",
70		signature = c(object = "data.frame"),
71		definition = function(object, ...,
72		method = c("berger", "boone", "brillouin",
73		"mcintosh", "shannon", "simpson")) {
74	6x	object <- data.matrix(object)
75	6x	methods::callGeneric(object, ..., method = method)
76		}
77		)
78
79		# Evenness =====================================================================
80		#' @export
81		#' @rdname heterogeneity
82		#' @aliases evenness,matrix-method
83		setMethod(
84		f = "evenness",
85		signature = c(object = "matrix"),
86		definition = function(object, ...,
87		method = c("shannon", "brillouin",
88		"mcintosh", "simpson")) {
89	4x	method <- match.arg(method, several.ok = FALSE)
90	4x	index <- index_diversity(object, method, ..., evenness = TRUE)
91	4x	.EvennessIndex(index)
92		}
93		)
94
95		#' @export
96		#' @rdname heterogeneity
97		#' @aliases evenness,data.frame-method
98		setMethod(
99		f = "evenness",
100		signature = c(object = "data.frame"),
101		definition = function(object, ...,
102		method = c("shannon", "brillouin",
103		"mcintosh", "simpson")) {
104	4x	object <- data.matrix(object)
105	4x	methods::callGeneric(object, ..., method = method)
106		}
107		)
108
109		# Richness =====================================================================
110		#' @export
111		#' @rdname richness
112		#' @aliases richness,matrix-method
113		setMethod(
114		f = "richness",
115		signature = c(object = "matrix"),
116		definition = function(object, ..., method = c("observed", "margalef", "menhinick")) {
117		## Backward compatibility
118	!	if (method == "count") method <- "observed"
119
120	4x	method <- match.arg(method, several.ok = FALSE)
121	4x	index <- index_diversity(object, method, ...)
122	4x	.RichnessIndex(index)
123		}
124		)
125
126		#' @export
127		#' @rdname richness
128		#' @aliases richness,data.frame-method
129		setMethod(
130		f = "richness",
131		signature = c(object = "data.frame"),
132		definition = function(object, ..., method = c("observed", "margalef", "menhinick")) {
133	4x	object <- data.matrix(object)
134	4x	methods::callGeneric(object, ..., method = method)
135		}
136		)
137
138		# Composition ==================================================================
139		#' @export
140		#' @rdname richness
141		#' @aliases composition,matrix-method
142		setMethod(
143		f = "composition",
144		signature = c(object = "matrix"),
145		definition = function(object, ...,
146		method = c("chao1", "ace", "squares", "chao2", "ice")) {
147	2x	method <- match.arg(method, several.ok = FALSE)
148	2x	by_row <- any(method == c("chao1", "ace", "squares"))
149	2x	index <- index_diversity(object, method, ..., by_row = by_row)
150	2x	.CompositionIndex(index)
151		}
152		)
153
154		#' @export
155		#' @rdname richness
156		#' @aliases composition,data.frame-method
157		setMethod(
158		f = "composition",
159		signature = c(object = "data.frame"),
160		definition = function(object, ...,
161		method = c("chao1", "ace", "squares", "chao2", "ice")) {
162	!	object <- data.matrix(object)
163	!	methods::callGeneric(object, ..., method = method)
164		}
165		)
166
167		# Turnover =====================================================================
168		#' @export
169		#' @rdname turnover
170		#' @aliases turnover,matrix-method
171		setMethod(
172		f = "turnover",
173		signature = c(object = "matrix"),
174		definition = function(object, ...,
175		method = c("whittaker", "cody", "routledge1",
176		"routledge2", "routledge3", "wilson")) {
177	6x	method <- match.arg(method, several.ok = FALSE)
178	6x	fun <- get_index(method)
179	6x	fun(object)
180		}
181		)
182
183		#' @export
184		#' @rdname turnover
185		#' @aliases turnover,data.frame-method
186		setMethod(
187		f = "turnover",
188		signature = c(object = "data.frame"),
189		definition = function(object, ...,
190		method = c("whittaker", "cody", "routledge1",
191		"routledge2", "routledge3", "wilson")) {
192	!	object <- data.matrix(object)
193	!	methods::callGeneric(object, ..., method = method)
194		}
195		)
196
197		# Resample =====================================================================
198		## Bootstrap -------------------------------------------------------------------
199		#' @export
200		#' @rdname bootstrap
201		#' @aliases bootstrap,DiversityIndex-method
202		setMethod(
203		f = "bootstrap",
204		signature = c(object = "DiversityIndex"),
205		definition = function(object, n = 1000, f = NULL) {
206
207	3x	w <- object@data
208	3x	m <- nrow(w)
209	3x	method <- object@method
210
211	3x	results <- vector(mode = "list", length = m)
212	3x	for (i in seq_len(m)) {
213	15x	results[[i]] <- arkhe::bootstrap(
214	15x	object = w[i, ],
215	15x	do = do_index,
216	15x	n = n,
217	15x	method = method,
218	15x	evenness = methods::is(object, "EvennessIndex"),
219	15x	f = f
220		)
221		}
222	3x	results <- do.call(rbind, results)
223	3x	rownames(results) <- rownames(w)
224	3x	as.data.frame(results)
225		}
226		)
227
228		## Jackknife -------------------------------------------------------------------
229		#' @export
230		#' @rdname jackknife
231		#' @aliases jackknife,DiversityIndex-method
232		setMethod(
233		f = "jackknife",
234		signature = c(object = "DiversityIndex"),
235		definition = function(object, f = NULL) {
236
237	3x	w <- object@data
238	3x	m <- nrow(w)
239	3x	method <- object@method
240
241	3x	results <- vector(mode = "list", length = m)
242	3x	for (i in seq_len(m)) {
243	15x	results[[i]] <- arkhe::jackknife(
244	15x	object = w[i, ],
245	15x	do = do_index,
246	15x	method = method,
247	15x	evenness = methods::is(object, "EvennessIndex"),
248	15x	f = f
249		)
250		}
251	3x	results <- do.call(rbind, results)
252	3x	rownames(results) <- rownames(w)
253	3x	as.data.frame(results)
254		}
255		)
256
257		## Simulate --------------------------------------------------------------------
258		#' @export
259		#' @rdname simulate
260		#' @aliases simulate,DiversityIndex-method
261		setMethod(
262		f = "simulate",
263		signature = c(object = "DiversityIndex"),
264		definition = function(object, n = 1000, step = 1,
265		interval = c("percentiles", "student", "normal"),
266		level = 0.80, progress = getOption("tabula.progress")) {
267		## Simulate
268		## Specify the probability for the classes
269	2x	data <- object@data
270	2x	method <- object@method # Select method
271
272		## Sample size
273	2x	size <- max(rowSums(data))
274	2x	sample_sizes <- seq(from = 1, to = size * 1.05, by = step)
275
276	2x	m <- length(sample_sizes)
277	2x	k <- seq_len(m)
278
279	2x	simulated <- vector(mode = "list", length = m)
280	2x	fun <- function(x) conf(x, level = level, type = interval)
281
282	2x	progress_bar <- interactive() && progress
283	!	if (progress_bar) pbar <- utils::txtProgressBar(max = m, style = 3)
284
285	2x	for (i in k) {
286	318x	simulated[[i]] <- resample(
287	318x	object = colSums(data),
288	318x	do = do_index,
289	318x	method = method,
290	318x	evenness = methods::is(object, "EvennessIndex"),
291	318x	n = n,
292	318x	size = sample_sizes[[i]],
293	318x	f = fun
294		)
295	!	if (progress_bar) utils::setTxtProgressBar(pbar, i)
296		}
297
298	!	if (progress_bar) close(pbar)
299
300	2x	simulated <- do.call(rbind, simulated)
301	2x	simulated <- cbind(size = sample_sizes, simulated)
302
303	2x	methods::initialize(object, simulation = simulated)
304		}
305		)
306
307		conf <- function(x, type = c("percentiles", "student", "normal"),
308		level = 0.80) {
309
310	318x	type <- match.arg(type, several.ok = FALSE)
311
312	318x	if (type == "percentiles") {
313		## Confidence interval as described in Kintigh 1989
314	318x	k <- (1 - level) / 2
315	318x	conf <- stats::quantile(x, probs = c(k, 1 - k), names = FALSE)
316		} else {
317		## Confidence interval
318	!	conf <- arkhe::confidence_mean(x, level = level, type = type)
319		}
320
321	318x	result <- c(mean(x), conf)
322	318x	names(result) <- c("mean", "lower", "upper")
323	318x	result
324		}

1		# PLOT BERTIN
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname plot_bertin
7		#' @aliases plot_bertin,matrix-method
8		setMethod(
9		f = "plot_bertin",
10		signature = signature(object = "matrix"),
11		definition = function(object, threshold = NULL, freq = FALSE, margin = 1,
12		color = c("white", "black"), flip = TRUE,
13		axes = TRUE, ...) {
14		## Conditional proportions
15	2x	if (freq) object <- prop.table(object, margin = margin)
16
17		## Compute threshold for each variable
18	5x	if (!freq && is.function(threshold)) {
19	1x	thr <- apply(X = object, MARGIN = 2, FUN = threshold)
20	1x	thr <- matrix(thr, nrow = nrow(object), ncol = ncol(object), byrow = TRUE)
21	1x	thr <- ifelse(object > thr, color[length(color)], color[1L])
22	1x	color <- as.vector(t(thr))
23		}
24
25		## /!\ Bertin plot flips x and y axis /!\
26	5x	data <- if (flip) t(object) else object
27
28	5x	plot_matrix(data, panel = panel_bertin, color = color,
29	5x	axes = axes, legend = FALSE, asp = NA)
30
31		## Legend
32		# TODO
33
34	5x	invisible(object)
35		}
36		)
37
38		#' @export
39		#' @rdname plot_bertin
40		#' @aliases plot_bertin,data.frame-method
41		setMethod(
42		f = "plot_bertin",
43		signature = signature(object = "data.frame"),
44		definition = function(object, threshold = NULL, freq = FALSE, margin = 1,
45		color = c("white", "black"), flip = TRUE,
46		axes = TRUE, ...) {
47	5x	object <- data.matrix(object)
48	5x	methods::callGeneric(object, threshold = threshold,
49	5x	freq = freq, margin = margin,
50	5x	color = color, flip = flip, axes = axes)
51		}
52		)

1		# PLOT RANK
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname plot_rank
7		#' @aliases plot_rank,matrix-method
8		setMethod(
9		f = "plot_rank",
10		signature = c(object = "matrix"),
11		definition = function(object, log = NULL,
12		color = NULL, symbol = FALSE,
13		xlab = NULL, ylab = NULL,
14		main = NULL, sub = NULL,
15		ann = graphics::par("ann"),
16		axes = TRUE, frame.plot = axes,
17		panel.first = NULL, panel.last = NULL,
18		legend = list(x = "topright"), ...) {
19		## Prepare data
20	2x	n <- nrow(object)
21	2x	lab <- rownames(object)
22	2x	object <- object / rowSums(object)
23	2x	rk <- t(apply(X = object, MARGIN = 1, FUN = function(x) rank(-x)))
24
25		## Graphical parameters
26	2x	lwd <- list(...)$lwd %\|\|% graphics::par("lwd")
27	2x	lty <- list(...)$lty %\|\|% graphics::par("lty")
28	2x	cex <- list(...)$cex %\|\|% graphics::par("cex")
29	2x	if (length(lwd) == 1) lwd <- rep(lwd, length.out = n)
30	2x	if (length(lty) == 1) lty <- rep(lty, length.out = n)
31	2x	if (length(cex) == 1) cex <- rep(cex, length.out = n)
32	1x	if (is.null(log)) log <- ""
33
34	2x	pch <- list(...)$pch %\|\|% 16
35	2x	if (length(pch) == 1) pch <- rep(pch, length.out = n)
36	!	if (!isFALSE(symbol)) pch <- khroma::palette_shape(symbol)(lab)
37
38	2x	col <- list(...)$col %\|\|% graphics::par("col")
39	2x	if (length(col) == 1) col <- rep(col, length.out = n)
40	2x	if (!isFALSE(color)) col <- khroma::palette_color_discrete(color)(lab)
41
42		## Open new window
43	2x	grDevices::dev.hold()
44	2x	on.exit(grDevices::dev.flush(), add = TRUE)
45	2x	graphics::plot.new()
46
47		## Set plotting coordinates
48	2x	xlim <- range(rk)
49	2x	ylim <- range(object)
50	2x	graphics::plot.window(xlim = xlim, ylim = ylim, log = log)
51
52		## Evaluate pre-plot expressions
53	2x	panel.first
54
55		## Plot
56	2x	for (i in seq_len(n)) {
57	10x	z <- order(rk[i, ])
58	10x	graphics::lines(x = rk[i, z], y = object[i, z],
59	10x	col = col[i], lty = lty[i], lwd = lwd[i])
60	10x	graphics::points(x = rk[i, z], y = object[i, z],
61	10x	col = col[i], pch = pch[i], cex = cex[i])
62		}
63
64		## Evaluate post-plot and pre-axis expressions
65	2x	panel.last
66
67		## Construct axis
68	2x	if (axes) {
69	2x	graphics::axis(side = 1, las = 1)
70	2x	graphics::axis(side = 2, las = 1)
71		}
72
73		## Plot frame
74	2x	if (frame.plot) {
75	2x	graphics::box()
76		}
77
78		## Add annotation
79	2x	if (ann) {
80	2x	xlab <- xlab %\|\|% tr_("Rank")
81	2x	ylab <- ylab %\|\|% tr_("Frequency")
82	2x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
83		}
84
85		## Legend
86	2x	if (is.list(legend) && length(legend) > 0) {
87	1x	args <- list(legend = lab, col = col, pch = pch,
88	1x	lty = lty, lwd = lwd, bty = "n")
89	1x	args <- utils::modifyList(args, legend)
90	1x	do.call(graphics::legend, args = args)
91		}
92
93	2x	invisible(object)
94		}
95		)
96
97		#' @export
98		#' @rdname plot_rank
99		#' @aliases plot_rank,data.frame-method
100		setMethod(
101		f = "plot_rank",
102		signature = signature(object = "data.frame"),
103		definition = function(object, log = NULL,
104		xlab = NULL, ylab = NULL,
105		main = NULL, sub = NULL,
106		ann = graphics::par("ann"),
107		axes = TRUE, frame.plot = axes,
108		panel.first = NULL, panel.last = NULL,
109		legend = list(x = "topright"), ...) {
110	2x	object <- data.matrix(object)
111	2x	methods::callGeneric(object, log = log,
112	2x	xlab = xlab, ylab = ylab,
113	2x	main = main, sub = sub, ann = ann, axes = axes,
114	2x	frame.plot = frame.plot, panel.first = panel.first,
115	2x	panel.last = panel.last, legend = legend, ...)
116		}
117		)

1		# PLOT DICE-LERASS
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname plot_diceleraas
7		#' @aliases plot_diceleraas,matrix-method
8		setMethod(
9		f = "plot_diceleraas",
10		signature = signature(object = "matrix"),
11		definition = function(object,
12		main = NULL, sub = NULL,
13		ann = graphics::par("ann"),
14		axes = TRUE, frame.plot = FALSE,
15		panel.first = NULL, panel.last = NULL, ...) {
16		## Prepare data
17	1x	object <- t(object)
18	1x	object[object == 0] <- NA
19	1x	n <- nrow(object)
20	1x	id <- seq_len(n)
21	1x	lab <- rownames(object) %\|\|% id
22
23		## Mean
24	1x	moy <- rowMeans(object, na.rm = TRUE)
25		## Standard deviation
26	1x	ec <- apply(X = object, MARGIN = 1, FUN = stats::sd, na.rm = TRUE)
27		## Standard error
28	1x	se <- ec / sqrt(rowSums(!is.na(object)))
29		## Range
30	1x	ran <- t(apply(X = object, MARGIN = 1, FUN = range, na.rm = TRUE))
31	1x	colnames(ran) <- c("min", "max")
32
33	1x	data <- data.frame(y = rev(id), mean = moy, sd = ec, se = se * 2, ran)
34
35		## Graphical parameters
36	1x	cex.axis <- list(...)$cex.axis %\|\|% graphics::par("cex.axis")
37	1x	col.axis <- list(...)$col.axis %\|\|% graphics::par("col.axis")
38	1x	font.axis <- list(...)$font.axis %\|\|% graphics::par("font.axis")
39	1x	lwd <- list(...)$lwd %\|\|% graphics::par("lwd")
40	1x	lty <- list(...)$lty %\|\|% graphics::par("lty")
41	1x	col <- list(...)$col %\|\|% c("black")
42	1x	if (length(lwd) < n) lwd <- rep(lwd, length.out = n)
43	1x	if (length(lty) < n) lty <- rep(lty, length.out = n)
44	1x	if (length(col) < n) col <- rep(col, length.out = n)
45
46		## Save and restore
47	1x	mar <- graphics::par("mar")
48	1x	mar[2] <- inch2line(lab, cex = cex.axis) + 0.5
49	1x	old_par <- graphics::par(mar = mar)
50	1x	on.exit(graphics::par(old_par))
51
52		## Open new window
53	1x	grDevices::dev.hold()
54	1x	on.exit(grDevices::dev.flush(), add = TRUE)
55	1x	graphics::plot.new()
56
57		## Set plotting coordinates
58	1x	xlim <- range(data$min, data$max, data$mean + data$sd,
59	1x	data$mean - data$sd, na.rm = TRUE)
60	1x	ylim <- c(1, n + 1)
61	1x	graphics::plot.window(xlim = xlim, ylim = ylim)
62
63		## Evaluate pre-plot expressions
64	1x	panel.first
65
66		## Plot
67	1x	for (i in id) {
68	5x	tmp <- data[i, ]
69	5x	graphics::polygon(
70	5x	x = tmp$sd * c(-1, -1, 1, 1) + tmp$mean,
71	5x	y = tmp$y + (1 / 3) * c(0, 1, 1, 0),
72	5x	col = "white", border = NA
73		)
74	5x	graphics::lines(
75	5x	x = tmp$sd * c(-1, -1, 1, 1) + tmp$mean,
76	5x	y = tmp$y + (1 / 3) * c(0, 1, 1, 0),
77	5x	col = col[i], lty = lty[i]
78		)
79	5x	graphics::polygon(
80	5x	x = tmp$se * c(-1, -1, 1, 1) + tmp$mean,
81	5x	y = tmp$y + (1 / 3) * c(0, 1, 1, 0),
82	5x	col = col[i], border = NA
83		)
84	5x	graphics::segments(
85	5x	x0 = tmp$mean, y0 = tmp$y,
86	5x	x1 = tmp$mean, y1 = tmp$y + 0.5,
87	5x	col = col[i], lty = lty[i]
88		)
89	5x	graphics::segments(
90	5x	x0 = tmp$min, y0 = tmp$y,
91	5x	x1 = tmp$max, y1 = tmp$y,
92	5x	col = col[i], lty = lty[i]
93		)
94		}
95
96		## Evaluate post-plot and pre-axis expressions
97	1x	panel.last
98
99		## Construct axis
100	1x	if (axes) {
101	1x	graphics::axis(side = 1, las = 1)
102	1x	graphics::mtext(rev(lab), side = 2, at = id, las = 1, padj = 0, line = 0,
103	1x	cex = cex.axis, col.axis = col.axis, font = font.axis)
104		}
105
106		## Plot frame
107	1x	if (frame.plot) {
108	!	graphics::box()
109		}
110
111		## Add annotation
112	1x	if (ann) {
113	1x	graphics::title(main = main, sub = sub,
114	1x	xlab = tr_("Absolute frequency"), ylab = NULL, ...)
115		}
116
117	1x	invisible(object)
118		}
119		)
120
121		#' @export
122		#' @rdname plot_diceleraas
123		#' @aliases plot_diceleraas,data.frame-method
124		setMethod(
125		f = "plot_diceleraas",
126		signature = signature(object = "data.frame"),
127		definition = function(object,
128		main = NULL, sub = NULL,
129		ann = graphics::par("ann"),
130		axes = TRUE, frame.plot = FALSE,
131		panel.first = NULL, panel.last = NULL, ...) {
132	1x	object <- data.matrix(object)
133	1x	methods::callGeneric(object,
134	1x	main = main, sub = sub, ann = ann,
135	1x	axes = axes, frame.plot = frame.plot,
136	1x	panel.first = panel.first, panel.last = panel.last, ...)
137		}
138		)

1		# MATRIGRAPH
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname matrigraph
7		#' @aliases pvi,matrix-method
8		setMethod(
9		f = "pvi",
10		signature = signature(object = "matrix"),
11		definition = function(object) {
12	3x	pvi <- object * 100 / expected(object)
13	3x	dimnames(pvi) <- dimnames(object)
14	3x	pvi
15		}
16		)
17
18		#' @export
19		#' @rdname matrigraph
20		#' @aliases pvi,data.frame-method
21		setMethod(
22		f = "pvi",
23		signature = signature(object = "data.frame"),
24		definition = function(object) {
25	1x	object <- data.matrix(object)
26	1x	methods::callGeneric(object)
27		}
28		)
29
30		#' @export
31		#' @rdname matrigraph
32		#' @aliases matrigraph,matrix-method
33		setMethod(
34		f = "matrigraph",
35		signature = signature(object = "matrix"),
36		definition = function(object, reverse = FALSE, axes = TRUE, ...) {
37	2x	plot_matrix(pvi(object), panel = panel_matrigraph, scale = FALSE,
38	2x	axes = axes, legend = FALSE, asp = 1, reverse = reverse)
39	2x	invisible(object)
40		}
41		)
42
43		#' @export
44		#' @rdname matrigraph
45		#' @aliases matrigraph,data.frame-method
46		setMethod(
47		f = "matrigraph",
48		signature = signature(object = "data.frame"),
49		definition = function(object, reverse = FALSE, ...) {
50	2x	object <- data.matrix(object)
51	2x	methods::callGeneric(object, reverse = reverse, ...)
52		}
53		)

1		# PLOT SPOT
2		#' @include AllClasses.R AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname plot_spot
7		#' @aliases plot_spot,matrix-method
8		setMethod(
9		f = "plot_spot",
10		signature = signature(object = "matrix"),
11		definition = function(object, type = c("ring", "plain"),
12		color = NULL,
13		diag = TRUE, upper = TRUE, lower = TRUE,
14		freq = FALSE, margin = 1,
15		axes = TRUE, legend = TRUE, ...) {
16		## Validation
17	7x	type <- match.arg(type, several.ok = FALSE)
18
19		## Backward compatibility
20	7x	col <- list(...)$col
21	!	if (!is.null(col) && is.null(color)) color <- col
22
23	7x	plot_matrix(object, panel = panel_spot, color = color,
24	7x	diag = diag, upper = upper, lower = lower,
25	7x	freq = freq, margin = margin, drop_zero = TRUE,
26	7x	axes = axes, legend = legend, type = type)
27
28	7x	invisible(object)
29		}
30		)
31
32		#' @export
33		#' @rdname plot_spot
34		#' @aliases plot_spot,data.frame-method
35		setMethod(
36		f = "plot_spot",
37		signature = signature(object = "data.frame"),
38		definition = function(object, type = c("ring", "plain"),
39		color = NULL,
40		diag = TRUE, upper = TRUE, lower = TRUE,
41		freq = FALSE, margin = 1,
42		axes = TRUE, legend = TRUE, ...) {
43	4x	object <- data.matrix(object)
44	4x	methods::callGeneric(object, type = type, color = color,
45	4x	diag = diag, upper = upper, lower = lower,
46	4x	freq = freq, margin = margin,
47	4x	axes = axes, legend = legend, ...)
48		}
49		)
50
51		#' @export
52		#' @rdname plot_spot
53		#' @aliases plot_spot,dist-method
54		setMethod(
55		f = "plot_spot",
56		signature = signature(object = "dist"),
57		definition = function(object, type = c("ring", "plain"),
58		color = NULL,
59		diag = FALSE, upper = FALSE, lower = !upper,
60		axes = TRUE, legend = TRUE, ...) {
61		# index_name <- attr(object, "method")
62	3x	object <- as.matrix(object)
63	3x	methods::callGeneric(object, type = type, color = color,
64	3x	diag = diag, upper = upper, lower = lower,
65	3x	axes = axes, legend = legend)
66		}
67		)

1		# PLOT HEATMAP
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname plot_heatmap
7		#' @aliases plot_heatmap,matrix-method
8		setMethod(
9		f = "plot_heatmap",
10		signature = signature(object = "matrix"),
11		definition = function(object, color = NULL,
12		diag = TRUE, upper = TRUE, lower = TRUE,
13		freq = FALSE, margin = 1, fixed_ratio = TRUE,
14		axes = TRUE, legend = TRUE, ...) {
15		## Backward compatibility
16	7x	col <- list(...)$col
17	!	if (!is.null(col) && is.null(color)) color <- col
18
19	7x	plot_matrix(object, panel = panel_tiles, color = color,
20	7x	diag = diag, upper = upper, lower = lower,
21	7x	freq = freq, margin = margin, drop_zero = FALSE,
22	7x	axes = axes, legend = legend, asp = fixed_ratio)
23
24	7x	invisible(object)
25		}
26		)
27
28		#' @export
29		#' @rdname plot_heatmap
30		#' @aliases plot_heatmap,data.frame-method
31		setMethod(
32		f = "plot_heatmap",
33		signature = signature(object = "data.frame"),
34		definition = function(object, color = NULL,
35		diag = TRUE, upper = TRUE, lower = TRUE,
36		freq = FALSE, margin = 1, fixed_ratio = TRUE,
37		axes = TRUE, legend = TRUE, ...) {
38	4x	object <- data.matrix(object)
39	4x	methods::callGeneric(object, color = color,
40	4x	diag = diag, upper = upper, lower = lower,
41	4x	freq = freq, margin = margin,
42	4x	fixed_ratio = fixed_ratio,
43	4x	axes = axes, legend = legend)
44		}
45		)
46
47		#' @export
48		#' @rdname plot_heatmap
49		#' @aliases plot_heatmap,dist-method
50		setMethod(
51		f = "plot_heatmap",
52		signature = signature(object = "dist"),
53		definition = function(object, color = NULL,
54		diag = FALSE, upper = FALSE, lower = !upper,
55		axes = TRUE, legend = TRUE, ...) {
56	2x	object <- as.matrix(object)
57	2x	methods::callGeneric(object, color = color,
58	2x	diag = diag, upper = upper, lower = lower,
59	2x	axes = axes, legend = legend)
60		}
61		)

1		# CO-OCCURRENCE
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname occurrence
7		#' @aliases occurrence,matrix-method
8		setMethod(
9		f = "occurrence",
10		signature = c(object = "matrix"),
11		definition = function(object, method = c("absolute", "relative", "binomial")) {
12		## Validation
13	3x	method <- match.arg(method, several.ok = FALSE)
14
15		## Pairwise comparison
16	3x	p <- ncol(object)
17	3x	labels <- colnames(object)
18
19	3x	if (method == "absolute" \|\| method == "relative") {
20	3x	incid <- object > 0
21	3x	fun <- function(x) sum(incid[, x[1]] + incid[, x[2]] == 2)
22		}
23	3x	if (method == "binomial") {
24	!	fun <- function(x) index_binomial(object[, x[1]], object[, x[2]])
25		}
26
27	3x	cbn <- utils::combn(p, m = 2, simplify = TRUE)
28	3x	index <- apply(X = cbn, MARGIN = 2, FUN = fun)
29
30		## Matrix of results
31	3x	mtx <- matrix(data = 0, nrow = p, ncol = p,
32	3x	dimnames = list(labels, labels))
33	3x	mtx[lower.tri(mtx, diag = FALSE)] <- index
34	3x	mtx <- t(mtx)
35	3x	mtx[lower.tri(mtx, diag = FALSE)] <- index
36
37	3x	if (method == "relative") {
38	!	mtx <- mtx / nrow(object)
39		}
40
41	3x	occ <- stats::as.dist(mtx)
42	3x	attr(occ, "total") <- nrow(object)
43	3x	occ
44		}
45		)
46
47		#' @export
48		#' @rdname occurrence
49		#' @aliases occurrence,data.frame-method
50		setMethod(
51		f = "occurrence",
52		signature = c(object = "data.frame"),
53		definition = function(object, method = c("absolute", "relative", "binomial")) {
54	1x	object <- data.matrix(object)
55	1x	methods::callGeneric(object, method = method)
56		}
57		)
58
59		## Binomial co-occurrence ------------------------------------------------------
60		#' @export
61		#' @rdname index_binomial
62		#' @aliases index_binomial,numeric,numeric-method
63		setMethod(
64		f = "index_binomial",
65		signature = c(x = "numeric", y = "numeric"),
66		definition = function(x, y) {
67		## Validation
68	1x	arkhe::assert_length(y, length(x))
69
70		## Total number of assemblages
71	1x	N <- length(x)
72		## Expected proportion of co-occurrences for artifact classes
73	1x	p <- sum(x > 0) * sum(y > 0) / N^2
74		## Number of observed co-occurence for artifact classes
75	1x	o <- sum((x > 0) + (y > 0) == 2)
76	1x	if (p == 1) {
77		## Avoid NaN generation
78		## TODO: warning ?
79	!	Cbi <- 0
80		} else {
81	1x	Cbi <- (o - N * p) / sqrt(N * p * (1 - p))
82		}
83	1x	Cbi
84		}
85		)

1		# COERCION
2		#' @include AllGenerics.R AllClasses.R
3		NULL
4
5		# To data.frame ================================================================
6		#' @method as.data.frame DiversityIndex
7		#' @export
8		as.data.frame.DiversityIndex <- function(x, ...) {
9	2x	data.frame(
10	2x	observed = apply(X = x@data, MARGIN = 1, FUN = observed),
11	2x	singleton = apply(X = x@data, MARGIN = 1, FUN = singleton),
12	2x	doubleton = apply(X = x@data, MARGIN = 1, FUN = doubleton),
13	2x	size = x@size,
14	2x	index = x@.Data,
15	2x	row.names = labels(x),
16	2x	stringsAsFactors = FALSE
17		)
18		}
19
20		#' @export
21		#' @rdname data.frame
22		#' @aliases as.data.frame,DiversityIndex-method
23		setMethod("as.data.frame", "DiversityIndex", as.data.frame.DiversityIndex)

1		# SHOW METHODS
2		#' @include AllClasses.R
3		NULL
4
5		setMethod(
6		f = "show",
7		signature = "DiversityIndex",
8		definition = function(object) {
9	!	methods::callGeneric(object@.Data)
10		}
11		)
12		setMethod(
13		f = "show",
14		signature = "RarefactionIndex",
15		definition = function(object) {
16	!	methods::callGeneric(object@.Data)
17		}
18		)