nexus coverage - 78.61%

Files
Source

# COERCION
#' @include AllGenerics.R
NULL

# To CompositionMatrix =========================================================
#' @export
#' @rdname as_composition
#' @aliases as_composition,numeric-method
setMethod(
  f = "as_composition",
  signature = c(from = "numeric"),
  definition = function(from) {
    from <- matrix(data = from, nrow = 1, ncol = length(from))
    methods::callGeneric(from)
  }
)

#' @export
#' @rdname as_composition
#' @aliases as_composition,matrix-method
setMethod(
  f = "as_composition",
  signature = c(from = "matrix"),
  definition = function(from) {
    ## Validation
    arkhe::assert_type(from, "numeric")

    ## Make row/column names
    lab <- make_names(x = NULL, n = nrow(from), prefix = "S")
    rownames(from) <- if (has_rownames(from)) rownames(from) else lab
    colnames(from) <- make_names(x = colnames(from), n = ncol(from), prefix = "V")

    ## Close
    totals <- rowSums(from, na.rm = TRUE)
    from <- from / totals

    .CompositionMatrix(from, totals = unname(totals))
  }
)

#' @export
#' @rdname as_composition
#' @aliases as_composition,data.frame-method
setMethod(
  f = "as_composition",
  signature = c(from = "data.frame"),
  definition = function(from, parts = NULL, groups = NULL, autodetect = TRUE,
                        verbose = getOption("nexus.verbose")) {
    ## Clean row/column names
    lab <- make_names(x = NULL, n = nrow(from), prefix = "S")
    rownames(from) <- if (has_rownames(from)) rownames(from) else lab
    colnames(from) <- make_names(x = colnames(from), n = ncol(from), prefix = "V")

    ## Remove non-numeric columns
    if (is.null(parts)) {
      if (isTRUE(autodetect)) {
        parts <- arkhe::detect(from, f = is.numeric, margin = 2)
        if (isTRUE(verbose)) {
          n <- sum(parts)
          what <- ngettext(n, "Found %g part (%s)", "Found %g parts (%s)")
          cols <- paste0(colnames(from)[parts], collapse = ", ")
          message(sprintf(what, n, cols))
        }
      } else {
        arkhe::assert_filled(parts)
      }
    } else {
      if (is.numeric(parts)) parts <- seq_len(ncol(from)) %in% parts
      if (is.character(parts)) parts <- colnames(from) %in% parts
    }
    coda <- from[, parts, drop = FALSE]
    arkhe::assert_filled(coda)

    ## Build matrix
    coda <- data.matrix(coda, rownames.force = NA)
    totals <- rowSums(coda, na.rm = TRUE)
    coda <- coda / totals

    z <- .CompositionMatrix(coda, totals = unname(totals))
    if (is.null(groups)) return(z)

    ## Group names
    grp <- from[, groups, drop = FALSE]
    group(z, by = grp, verbose = verbose)
  }
)

# To amounts ===================================================================
#' @export
#' @rdname as_amounts
#' @aliases as_amounts,CompositionMatrix-method
setMethod(
  f = "as_amounts",
  signature = c(from = "CompositionMatrix"),
  definition = function(from) {
    methods::as(from, "matrix") * totals(from)
  }
)

# To data.frame ================================================================
#' @method as.data.frame CompositionMatrix
#' @export
as.data.frame.CompositionMatrix <- function(x, row.names = rownames(x),
                                            optional = FALSE, ...) {
  as.data.frame(methods::as(x, "matrix"), row.names = row.names, optional = optional)
}

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

#' @method as.data.frame GroupedComposition
#' @export
as.data.frame.GroupedComposition <- function(x, row.names = rownames(x),
                                             optional = FALSE, ...,
                                             group_var = ".group",
                                             group_after = 0) {
  z <- data.frame(
    methods::as(x, "matrix"),
    row.names = row.names,
    check.names = !optional
  )
  arkhe::append_column(z, group_names(x), after = group_after, var = group_var)
}

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

#' @method as.data.frame LogRatio
#' @export
as.data.frame.LogRatio <- function(x, row.names = rownames(x),
                                   optional = FALSE, ...) {
  z <- as.data.frame(methods::as(x, "matrix"), row.names = row.names, optional = optional)

  lab <- labels(x)
  for (j in seq_len(ncol(z))) attr(z[[j]], "label") <- lab[[j]]
  z
}

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

#' @method as.data.frame GroupedLogRatio
#' @export
as.data.frame.GroupedLogRatio <- function(x, row.names = rownames(x),
                                          optional = FALSE, ...,
                                          group_var = ".group",
                                          group_after = 0) {
  z <- data.frame(
    methods::as(x, "matrix"),
    row.names = row.names,
    check.names = !optional
  )

  lab <- labels(x)
  for (j in seq_len(ncol(z))) attr(z[[j]], "label") <- lab[[j]]

  arkhe::append_column(z, group_names(x), after = group_after, var = group_var)
}

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

#' @method as.data.frame OutlierIndex
#' @export
as.data.frame.OutlierIndex <- function(x, row.names = rownames(x),
                                       optional = FALSE, ...) {
  as.data.frame(x@standard, row.names = row.names, optional = optional)
}

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

# OPERATORS
#' @include AllGenerics.R
NULL

# Closure ======================================================================
#' @export
#' @rdname closure
#' @aliases closure,numeric-method
setMethod(
  f = "closure",
  signature = c(x = "numeric"),
  definition = function(x, total = 1, na.rm = TRUE) {
    x * total / sum(x, na.rm = na.rm)
  }
)

#' @export
#' @rdname closure
#' @aliases closure,matrix-method
setMethod(
  f = "closure",
  signature = c(x = "matrix"),
  definition = function(x, total = 1, na.rm = TRUE) {
    x * total / rowSums(x, na.rm = na.rm)
  }
)

# Perturbation =================================================================
#' @export
#' @rdname arithmetic
#' @aliases `%perturbe%`,CompositionMatrix,CompositionMatrix-method
setMethod(
  f = "%perturbe%",
  signature = c(x = "CompositionMatrix", y = "CompositionMatrix"),
  definition = function(x, y) {
    arkhe::assert_dimensions(y, dim(x))

    if (all(x <= 0)) x <- 1 / x
    if (all(y <= 0)) y <- 1 / y

    z <- x * y
    z <- as_composition(z)
    rownames(z) <- rownames(x)

    if (is_grouped(x)) {
      x@group_indices <- group_indices(x)
      x@group_levels <- group_levels(x)
    }

    z
  }
)

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

#' @export
#' @rdname perturbation
#' @aliases perturbation,CompositionMatrix,numeric-method
setMethod(
  f = "perturbation",
  signature = c(x = "CompositionMatrix", y = "numeric"),
  definition = function(x, y) {
    y <- matrix(data = y, nrow = nrow(x), ncol = length(y), byrow = TRUE)
    x %perturbe% as_composition(y)
  }
)

#' @export
#' @rdname perturbation
#' @aliases perturbation,CompositionMatrix,CompositionMatrix-method
setMethod(
  f = "perturbation",
  signature = c(x = "CompositionMatrix", y = "matrix"),
  definition = function(x, y) {
    x %perturbe% as_composition(y)
  }
)

# Powering =====================================================================
#' @export
#' @rdname arithmetic
#' @aliases `%power%`,CompositionMatrix,numeric-method
setMethod(
  f = "%power%",
  signature = c(x = "CompositionMatrix", y = "numeric"),
  definition = function(x, y) {
    arkhe::assert_length(y, 1L)

    z <- x ^ y
    z <- as_composition(z)
    rownames(z) <- rownames(x)

    if (is_grouped(x)) {
      x@group_indices <- group_indices(x)
      x@group_levels <- group_levels(x)
    }

    z
  }
)

#' @export
#' @rdname arithmetic
#' @aliases `%power%`,numeric,CompositionMatrix-method
setMethod(
  f = "%power%",
  signature = c(x = "numeric", y = "CompositionMatrix"),
  definition = function(x, y) {
    methods::callGeneric(x = y, y = x)
  }
)

#' @export
#' @rdname powering
#' @aliases powering,numeric,numeric-method
setMethod(
  f = "powering",
  signature = c(x = "numeric", a = "numeric"),
  definition = function(x, a) {
    arkhe::assert_length(a, 1L)
    closure(x ^ a)
  }
)

#' @export
#' @rdname powering
#' @aliases powering,numeric,numeric-method
setMethod(
  f = "powering",
  signature = c(x = "CompositionMatrix", a = "numeric"),
  definition = function(x, a) {
    x %power% a
  }
)

# Scalar product ===============================================================
#' @export
#' @rdname scalar
#' @aliases scalar,numeric,numeric-method
setMethod(
  f = "scalar",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    n <- length(x)
    arkhe::assert_length(y, n)

    D <- seq_len(n)
    z <- 0
    for (i in D) {
      j <- utils::tail(D, -i)
      z <- z + sum(log(x[i] / x[j]) * log(y[i] / y[j]))
    }
    (1 / n) * z
  }
)

#' @export
#' @rdname scalar
#' @aliases scalar,CompositionMatrix,CompositionMatrix-method
setMethod(
  f = "scalar",
  signature = c(x = "CompositionMatrix", y = "CompositionMatrix"),
  definition = function(x, y) {
    arkhe::assert_dimensions(y, dim(x))
    m <- nrow(x)

    z <- numeric(m)
    for (i in seq_len(m)) {
      z[i] <- scalar(x[i, , drop = TRUE], y[i, , drop = TRUE])
    }
    z
  }
)

#' Norm of a Composition
#'
#' @param x A [`CompositionMatrix-class`] object.
#' @return A [`numeric`] vector.
#' @keywords internal
#' @noRd
norm <- function(x) {
  sqrt(scalar(x, x))
}

#' Aitchison Distance
#'
#' @param x A [`CompositionMatrix-class`] object.
#' @param diag A [`logical`] scalar: should the diagonal of the distance matrix
#'  be printed?
#' @param upper A [`logical`] scalar: should the upper triangle of the distance
#'  matrix be printed?
#' @return A [`dist`] object.
#' @keywords internal
#' @noRd
aitchison <- function(x, diag = FALSE, upper = FALSE) {
  m <- nrow(x)
  spl <- rownames(x)

  d <- utils::combn(
    x = seq_len(m),
    m = 2,
    FUN = function(i, coda) {
      x <- coda[i[1], , drop = TRUE]
      y <- coda[i[2], , drop = TRUE]
      norm(x / y)
    },
    coda = x
  )

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

  stats::as.dist(mtx, diag = diag, upper = upper)
}

# STATISTICS: MEAN
#' @include AllGenerics.R
NULL

# Mean =========================================================================
#' @export
#' @method mean CompositionMatrix
mean.CompositionMatrix <- function(x, ..., ignore_na = FALSE, ignore_zero = TRUE) {
  m <- apply(
    X = x,
    MARGIN = 2,
    FUN = gmean,
    ignore_na = ignore_na,
    ignore_zero = ignore_zero,
    simplify = TRUE
  )
  m <- closure(m)
  names(m) <- colnames(x)
  m
}

#' @export
#' @rdname mean
#' @aliases mean,CompositionMatrix-method
setMethod("mean", "CompositionMatrix", mean.CompositionMatrix)

#' Geometric Mean
#'
#' @param x A [`numeric`] vector.
#' @param trim A length-one [`numeric`] vector specifying the fraction (0 to 0.5)
#'  of observations to be trimmed from each end of `x` before the mean is
#'  computed.
#' @param ignore_na A [`logical`] scalar: should [missing values][NA] be
#'  stripped before the computation proceeds?
#' @param ignore_zero A [`logical`] scalar: should zeros be stripped before the
#'  computation proceeds?
#' @return A [`numeric`] vector.
#' @keywords internal
gmean <- function(x, trim = 0, ignore_na = FALSE, ignore_zero = TRUE) {
  if (ignore_na) x <- x[is.finite(x)]
  if (ignore_zero) x <- x[x > 0]
  exp(mean(log(unclass(x)), trim = trim))
}

# GROUPS
#' @include AllGenerics.R
NULL

.group_extract <- function(object, which) {
  ## Validation
  arkhe::assert_type(which, "character")
  if (!any_assigned(object)) stop("No group is defined.", call. = FALSE)

  ok <- group_names(object) %in% which
  if (!any(ok)) {
    msg <- ngettext(length(which), "No sample belongs to the group: %s.",
                    "No sample belongs to the groups: %s.")
    message(sprintf(msg, paste0(dQuote(which), collapse = ", ")))
    return(object)
  }

  object[ok, , drop = FALSE]
}

#' @export
#' @rdname group_subset
#' @aliases group_subset,GroupedComposition-method
setMethod(
  f = "group_subset",
  signature = c("GroupedComposition"),
  definition = .group_extract
)

#' @export
#' @rdname group_subset
#' @aliases group_subset,GroupedLogRatio-method
setMethod(
  f = "group_subset",
  signature = c("GroupedLogRatio"),
  definition = .group_extract
)

# Groups =======================================================================
#' Compute Groups
#'
#' @param x A ([`list`] of) [`factors`] for which interaction is to be computed.
#' @param drop_levels A [`logical`] scalar: should unused factor levels be
#'  dropped?
#' @param allow_na A [`logical`] scalar: should `NA` be considered an extra
#'  level?
#' @return A [`factor`] vector.
#' @keywords internal
#' @noRd
compute_groups <- function(x, drop_levels = TRUE, allow_na = TRUE) {
  if (!is.list(x)) x <- list(x)
  x <- interaction(x, sep = "_")
  if (drop_levels) x <- droplevels(x)
  if (allow_na) x <- addNA(x, ifany = TRUE)

  x
}

#' @export
#' @rdname group
#' @aliases group,CompositionMatrix-method
setMethod(
  f = "group",
  signature = "CompositionMatrix",
  definition = function(object, by, verbose = getOption("nexus.verbose"), ...) {
    ## Compute groups
    by <- compute_groups(by, ...)

    ## Validation
    arkhe::assert_length(by, nrow(object))
    if (nlevels(by) == 0) {
      stop(tr_("Nothing to group by."), call. = FALSE)
    }
    if (isTRUE(verbose)) {
      if (nlevels(by) == nrow(object)) {
        message(tr_("As many groups as individuals."))
      }

      n <- nlevels(by)
      what <- ngettext(n, "Found %g group (%s)", "Found %g groups (%s)")
      grp <- paste0(levels(by), collapse = ", ")
      message(sprintf(what, n, grp))
    }

    .GroupedComposition(
      object,
      group_indices = as.integer(by),
      group_levels = levels(by)
    )
  }
)

#' @export
#' @rdname group
#' @aliases group,GroupedComposition-method
setMethod(
  f = "group",
  signature = "GroupedComposition",
  definition = function(object, by, add = FALSE,
                        verbose = getOption("nexus.verbose"), ...) {
    ## Compute groups
    if (isTRUE(add)) {
      if (!is.list(by)) by <- list(by)
      by <- c(list(group_factor(object)), by)
    }
    methods::callNextMethod(object, by = by, verbose = verbose, ...)
  }
)

#' @export
#' @rdname group
#' @aliases ungroup,GroupedComposition-method
setMethod(
  f = "ungroup",
  signature = "GroupedComposition",
  definition = function(object) {
    methods::as(object, "CompositionMatrix", strict = TRUE)
  }
)

#' @export
#' @rdname group
#' @aliases ungroup,GroupedLR-method
setMethod(
  f = "ungroup",
  signature = "GroupedLR",
  definition = function(object) {
    methods::as(object, "LR", strict = TRUE)
  }
)

#' @export
#' @rdname group
#' @aliases ungroup,GroupedCLR-method
setMethod(
  f = "ungroup",
  signature = "GroupedCLR",
  definition = function(object) {
    methods::as(object, "CLR", strict = TRUE)
  }
)

#' @export
#' @rdname group
#' @aliases ungroup,GroupedALR-method
setMethod(
  f = "ungroup",
  signature = "GroupedALR",
  definition = function(object) {
    methods::as(object, "ALR", strict = TRUE)
  }
)

#' @export
#' @rdname group
#' @aliases ungroup,GroupedILR-method
setMethod(
  f = "ungroup",
  signature = "GroupedILR",
  definition = function(object) {
    methods::as(object, "ILR", strict = TRUE)
  }
)

#' @export
#' @rdname group
#' @aliases ungroup,GroupedPLR-method
setMethod(
  f = "ungroup",
  signature = "GroupedPLR",
  definition = function(object) {
    methods::as(object, "PLR", strict = TRUE)
  }
)

# Metadata =====================================================================
#' @export
#' @describeIn group_metadata returns a [`character`] vector giving the group
#'  names.
#' @aliases group_levels,ReferenceGroups-method
setMethod(
  f = "group_levels",
  signature = "ReferenceGroups",
  definition = function(object) object@group_levels
)

#' @export
#' @describeIn group_metadata returns a [`character`] vector giving the name of
#'  the group that each observation belongs to.
#' @aliases group_names,ReferenceGroups-method
setMethod(
  f = "group_names",
  signature = "ReferenceGroups",
  definition = function(object) group_levels(object)[group_indices(object)]
)

group_factor <- function(object) {
  factor(
    x = group_names(object),
    levels = group_levels(object),
    exclude = NULL
  )
}

#' @export
#' @describeIn group_metadata returns an [`integer`] vector giving the group
#'  that each value belongs to.
#' @aliases group_indices,ReferenceGroups-method
setMethod(
  f = "group_indices",
  signature = "ReferenceGroups",
  definition = function(object) object@group_indices
)

#' @export
#' @describeIn group_metadata returns a `list` of [`integer`] vectors giving the
#'  observation that each group contains.
#' @aliases group_rows,ReferenceGroups-method
setMethod(
  f = "group_rows",
  signature = "ReferenceGroups",
  definition = function(object) {
    i <- group_factor(object)
    split(seq_along(i), f = i)
  }
)

#' @export
#' @describeIn group_metadata gives the total number of groups.
#' @aliases group_length,ReferenceGroups-method
setMethod(
  f = "group_length",
  signature = "ReferenceGroups",
  definition = function(object) length(group_levels(object))
)

#' @export
#' @describeIn group_metadata gives the size of each group.
#' @aliases group_size,ReferenceGroups-method
setMethod(
  f = "group_size",
  signature = "ReferenceGroups",
  definition = function(object) lengths(group_rows(object))
)

#' @export
#' @describeIn group_metadata returns a [`logical`] vector specifying whether or
#'  not an observation belongs to a group.
#' @aliases is_assigned,ReferenceGroups-method
setMethod(
  f = "is_assigned",
  signature = "ReferenceGroups",
  definition = function(object) !is.na(group_names(object))
)

#' @export
#' @describeIn group_metadata returns an [`logical`] scalar specifying if any
#'  observation belongs to a group.
#' @aliases any_assigned,ReferenceGroups-method
setMethod(
  f = "any_assigned",
  signature = "ReferenceGroups",
  definition = function(object) any(is_assigned(object))
)

#' @export
#' @describeIn group_metadata returns an [`logical`] scalar specifying if all
#'  observations belong to a group.
#' @aliases all_assigned,ReferenceGroups-method
setMethod(
  f = "all_assigned",
  signature = "ReferenceGroups",
  definition = function(object) all(is_assigned(object))
)

# HISTOGRAM
#' @include AllGenerics.R
NULL

# CompositionMatrix ============================================================
#' @export
#' @method hist CompositionMatrix
hist.CompositionMatrix <- function(x, ..., select = 1,
                                   breaks = "Sturges",
                                   freq = FALSE, labels = FALSE,
                                   main = NULL, sub = NULL,
                                   ann = graphics::par("ann"),
                                   axes = TRUE, frame.plot = axes) {
  ## Validation
  if (is.character(select)) select <- match(select, labels(x))
  arkhe::assert_length(select, 1)

  ## Graphical parameters
  cex.lab <- list(...)$cex.lab %||% graphics::par("cex.lab")
  col.lab <- list(...)$col.lab %||% graphics::par("col.lab")
  font.lab <- list(...)$font.lab %||% graphics::par("font.lab")
  cex.main <- list(...)$cex.main %||% graphics::par("cex.main")
  col.main <- list(...)$col.main %||% graphics::par("col.main")
  font.main <- list(...)$font.main %||% graphics::par("font.main")

  ## Compute univariate ilr transformation
  z <- univariate_ilr(x)

  ## Select one compositonal part
  xi <- x[, select, drop = TRUE]
  zi <- z[, select, drop = TRUE]

  ## Compute axis in percent
  lab_i <- pretty(xi, n = 6)
  lab_i <- lab_i[lab_i > 0]
  at_i <- univariate_ilr(lab_i)

  ## Plot histogram
  h <- graphics::hist(x = zi, breaks = breaks, plot = FALSE)
  xlim <- range(at_i, h$breaks, finite = TRUE)
  plot(h, freq = freq, xlim = xlim, labels = labels, ...,
       main = main, sub = sub, xlab = NULL, ylab = NULL, axes = FALSE)

  ## Construct axis
  if (axes) {
    graphics::axis(side = 1, xpd = NA, las = 1)
    graphics::axis(side = 3, at = at_i, labels = label_percent(lab_i),
                   xpd = NA, las = 1)
    graphics::axis(side = 2, xpd = NA, las = 1)
  }

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

  ## Add annotation
  if (ann) {
    xlab <- labels(x)[select]
    ylab <- tr_("Frequency")
    graphics::mtext(sprintf("ilr(%s)", xlab), side = 1, line = 3,
                    cex = cex.lab, col = col.lab, font = font.lab)
    graphics::mtext(sprintf("%s %%", xlab), side = 3, line = 3,
                    cex = cex.lab, col = col.lab, font = font.lab)
    graphics::mtext(ylab, side = 2, line = 3,
                    cex = cex.lab, col = col.lab, font = font.lab)
  }

  invisible(x)
}

#' @export
#' @rdname hist
#' @aliases hist,CompositionMatrix-method
setMethod("hist", c(x = "CompositionMatrix"), hist.CompositionMatrix)

# ACCESSORS
#' @include AllGenerics.R
NULL

# Extract ======================================================================
.subscript1 <- function(x, i) {
  x@.Data[i]
}

.subscript2 <- function(x, i, j, drop) {
  ## Rows
  if (missing(i)) i <- seq_len(nrow(x))
  if (is.character(i)) i <- match(i, dimnames(x)[1L])
  totals <- totals(x)[i]

  ## Columns
  if (missing(j)) j <- seq_len(ncol(x))

  ## Subset
  z <- x@.Data[i, j, drop = drop]
  if (drop) return(z)

  ## /!\ Subcomposition /!\
  # if (ncol(z) < ncol(x)) {
  #   tot <- rowSums(z, na.rm = TRUE)
  #   totals <- totals * tot
  #   z <- z / tot
  # }

  if (is_grouped(x)) {
    g <- droplevels(group_factor(x)[i])
    methods::initialize(
      x, z,
      totals = totals,
      group_indices = as.integer(g),
      group_levels = levels(g)
    )
  } else {
    methods::initialize(x, z, totals = totals)
  }
}

wrong_dimensions <- function(i, j) {
  msg <- sprintf("M[%s%s]: incorrect number of dimensions.", i, j)
  stop(msg, call. = FALSE)
}

## CompositionMatrix -----------------------------------------------------------
#' @export
#' @rdname subset
#' @aliases [,CompositionMatrix,missing,missing,missing-method
setMethod(
  f = "[",
  signature = c(x = "CompositionMatrix", i = "missing", j = "missing", drop = "missing"),
  definition = function(x, i, j, ..., drop) {
    na <- nargs()
    if (na == 2L) return(x) # x[]
    if (na == 3L) return(x) # x[, ]
    wrong_dimensions(".", ".") # x[, , ], etc.
  }
)

#' @export
#' @rdname subset
#' @aliases [,CompositionMatrix,missing,missing,logical-method
setMethod(
  f = "[",
  signature = c(x = "CompositionMatrix", i = "missing", j = "missing", drop = "logical"),
  definition = function(x, i, j, ..., drop) {
    na <- nargs()
    if (na < 4L) return(x) # x[drop=], x[, drop=], x[drop=, ]
    if (na == 4L) {
      x <- if (drop) x@.Data else x
      return(x) # x[, , drop=], x[, drop=, ], x[drop=, , ]
    }
    wrong_dimensions(".", ".") # x[, , , drop=], etc.
  }
)

#' @export
#' @rdname subset
#' @aliases [,CompositionMatrix,index,missing,missing-method
setMethod(
  f = "[",
  signature = c(x = "CompositionMatrix", i = "index", j = "missing", drop = "missing"),
  definition = function(x, i, j, ..., drop) {
    na <- nargs()
    if (na == 2L) {
      x <- .subscript1(x, i)
      return(x) # x[i=]
    }
    if (na == 3L) {
      #/!\ DROP /!\
      x <- .subscript2(x, i, , drop = FALSE)
      return(x) # x[i=, ], x[, i=]
    }
    wrong_dimensions("i", ".") # x[i=, , ], etc.
  }
)

#' @export
#' @rdname subset
#' @aliases [,CompositionMatrix,index,missing,logical-method
setMethod(
  f = "[",
  signature = c(x = "CompositionMatrix", i = "index", j = "missing", drop = "logical"),
  definition = function(x, i, j, ..., drop) {
    na <- nargs()
    if (na == 3L) {
      x <- .subscript1(x, i)
      return(x) # x[i=, drop=]
    }
    if (na == 4L) {
      x <- .subscript2(x, i, , drop = drop)
      return(x) # x[i=, , drop=], x[, i=, drop=]
    }
    wrong_dimensions("i", ".") # x[i=, , , drop=], etc.
  }
)

#' @export
#' @rdname subset
#' @aliases [,CompositionMatrix,missing,index,missing-method
setMethod(
  f = "[",
  signature = c(x = "CompositionMatrix", i = "missing", j = "index", drop = "missing"),
  definition = function(x, i, j, ..., drop) {
    na <- nargs()
    if (na == 2L) {
      x <- .subscript1(x, j) # x[j=]
      return(x)
    }
    if (na == 3L) {
      # /!\ DROP /!\
      x <- .subscript2(x, , j, drop = FALSE) # x[j=, ], x[, j=]
      return(x)
    }
    wrong_dimensions(".", "j") # x[, j=, ], etc.
  }
)

#' @export
#' @rdname subset
#' @aliases [,CompositionMatrix,missing,index,logical-method
setMethod(
  f = "[",
  signature = c(x = "CompositionMatrix", i = "missing", j = "index", drop = "logical"),
  definition = function(x, i, j, ..., drop) {
    na <- nargs()
    if (na == 3L) {
      x <- .subscript1(x, j) # x[j=, drop=]
      return(x)
    }
    if (na == 4L) {
      x <- .subscript2(x, , j, drop = drop) # x[j=, , drop=], x[, j=, drop=]
      return(x)
    }
    wrong_dimensions(".", "j") # x[, j=, , drop=], etc.
  }
)

#' @export
#' @rdname subset
#' @aliases [,CompositionMatrix,index,index,missing-method
setMethod(
  f = "[",
  signature = c(x = "CompositionMatrix", i = "index", j = "index", drop = "missing"),
  definition = function(x, i, j, ..., drop) {
    na <- nargs()
    if (na == 3L) {
      # /!\ DROP /!\
      x <- .subscript2(x, i, j, drop = FALSE) # x[i=, j=], x[j=, i=]
      return(x)
    }
    wrong_dimensions("i", "j") # x[i=, j=, ], etc.
  }
)

#' @export
#' @rdname subset
#' @aliases [,CompositionMatrix,index,index,logical-method
setMethod(
  f = "[",
  signature = c(x = "CompositionMatrix", i = "index", j = "index", drop = "logical"),
  definition = function(x, i, j, ..., drop) {
    na <- nargs()
    if (na == 4L) {
      x <- .subscript2(x, i, j, drop = drop) # x[i=, j=, drop=], x[j=, i=, drop=]
      return(x)
    }
    wrong_dimensions("i", "j") # x[i=, j=, , drop=], etc.
  }
)

# Replace ======================================================================
## [<- -------------------------------------------------------------------------
#' @export
#' @rdname subset
#' @aliases [<-,CompositionMatrix-method
setMethod(
  f = "[<-",
  signature = c(x = "CompositionMatrix"),
  function(x, i, j, ..., value) {
    z <- methods::callNextMethod()
    methods::validObject(z)
    z
  }
)

## [[<- ------------------------------------------------------------------------
#' @export
#' @rdname subset
#' @aliases [[<-,CompositionMatrix-method
setMethod(
  f = "[[<-",
  signature = c(x = "CompositionMatrix"),
  function(x, i, j, ..., value) {
    z <- methods::callNextMethod()
    methods::validObject(z)
    z
  }
)

# Transpose ====================================================================
#' @export
#' @rdname t
#' @aliases t,CompositionMatrix-method
setMethod(
  f = "t",
  signature = c(x = "CompositionMatrix"),
  function(x) {
    t(x@.Data)
  }
)

#' @export
#' @rdname t
#' @aliases t,LogRatio-method
setMethod(
  f = "t",
  signature = c(x = "LogRatio"),
  function(x) {
    t(x@.Data)
  }
)

# VARIATION MATRIX
#' @include AllGenerics.R
NULL

# Variation matrix =============================================================
#' @export
#' @rdname variation
#' @aliases variation,CompositionMatrix-method
setMethod(
  f = "variation",
  signature = c(x = "CompositionMatrix"),
  definition = function(x) {
    J <- ncol(x)
    parts <- colnames(x)

    varia <- utils::combn(
      x = seq_len(J),
      m = 2,
      FUN = function(i, coda) {
        r <- coda[, i[[1]], drop = TRUE] / coda[, i[[2]], drop = TRUE]
        z <- log(r, base = exp(1))
        stats::var(z)
      },
      coda = x
    )

    mtx <- matrix(data = 0, nrow = J, ncol = J)
    mtx[lower.tri(mtx, diag = FALSE)] <- varia
    mtx <- t(mtx)
    mtx[lower.tri(mtx, diag = FALSE)] <- varia

    dimnames(mtx) <- list(parts, parts)
    mtx
  }
)

#' @export
#' @rdname pip
#' @aliases pip,CompositionMatrix-method
setMethod(
  f = "pip",
  signature = c(x = "CompositionMatrix"),
  definition = function(x) {
    v <- variation(x)
    1 / (1 + sqrt(v))
  }
)

# Variation array ==============================================================
# @export
# @rdname variation_array
# @aliases variation_array,CompositionMatrix-method
# setMethod(
#   f = "variation_array",
#   signature = c(object = "CompositionMatrix"),
#   definition = function(object) {
#     J <- ncol(object)
#     cbn <- utils::combn(seq_len(J), 2)
#     varia <- apply(
#       X = cbn,
#       MARGIN = 2,
#       FUN = function(j, x) {
#         mean(log(x[, j[1]] / x[, j[2]]))
#       },
#       x = object
#     )
#
#     mtx <- variation(object)
#     mtx[lower.tri(mtx, diag = FALSE)] <- varia
#     mtx
#   }
# )

# DATA TRANSFORMATION: INVERSE
#' @include AllGenerics.R
NULL

# Backtransform ================================================================
## CLR -------------------------------------------------------------------------
#' @export
#' @rdname transform_inverse
#' @aliases transform_inverse,CLR,missing-method
setMethod(
  f = "transform_inverse",
  signature = c(object = "CLR", origin = "missing"),
  definition = function(object) {
    y <- methods::as(object, "matrix") # Drop slots
    y <- exp(y)
    y <- y / rowSums(y)

    dimnames(y) <- list(rownames(object), object@parts)
    .CompositionMatrix(y, totals = totals(object))
  }
)

#' @export
#' @rdname transform_inverse
#' @aliases transform_inverse,GroupedCLR,missing-method
setMethod(
  f = "transform_inverse",
  signature = c(object = "GroupedCLR", origin = "missing"),
  definition = function(object) {
    z <- methods::callNextMethod()
    .GroupedComposition(z, group_indices = group_indices(object),
                        group_levels = group_levels(object))
  }
)

## ALR -------------------------------------------------------------------------
#' @export
#' @rdname transform_inverse
#' @aliases transform_inverse,ALR,missing-method
setMethod(
  f = "transform_inverse",
  signature = c(object = "ALR", origin = "missing"),
  definition = function(object) {
    y <- exp(object)
    y <- y / (1 + rowSums(y))
    z <- 1 - rowSums(y)

    y <- cbind(y, z)
    dimnames(y) <- list(rownames(object), object@parts)
    y <- y[, object@order]

    .CompositionMatrix(y, totals = totals(object))
  }
)

#' @export
#' @rdname transform_inverse
#' @aliases transform_inverse,GroupedALR,missing-method
setMethod(
  f = "transform_inverse",
  signature = c(object = "GroupedALR", origin = "missing"),
  definition = function(object) {
    z <- methods::callNextMethod()
    .GroupedComposition(z, group_indices = group_indices(object),
                        group_levels = group_levels(object))
  }
)

## ILR -------------------------------------------------------------------------
#' @export
#' @rdname transform_inverse
#' @aliases transform_inverse,ILR,missing-method
setMethod(
  f = "transform_inverse",
  signature = c(object = "ILR", origin = "missing"),
  definition = function(object) {
    y <- tcrossprod(object@.Data, object@base)
    y <- exp(y)
    y <- y / rowSums(y)

    dimnames(y) <- list(rownames(object), object@parts)
    y <- y[, object@order]

    .CompositionMatrix(y, totals = totals(object))
  }
)

#' @export
#' @rdname transform_inverse
#' @aliases transform_inverse,GroupedILR,missing-method
setMethod(
  f = "transform_inverse",
  signature = c(object = "GroupedILR", origin = "missing"),
  definition = function(object) {
    z <- methods::callNextMethod()
    .GroupedComposition(z, group_indices = group_indices(object),
                        group_levels = group_levels(object))
  }
)

#' @export
#' @rdname transform_inverse
#' @aliases transform_inverse,matrix,ILR-method
setMethod(
  f = "transform_inverse",
  signature = c(object = "matrix", origin = "ILR"),
  definition = function(object, origin) {
    y <- tcrossprod(object, origin@base)
    y <- exp(y)
    y <- y / rowSums(y)

    dimnames(y) <- list(rownames(object), origin@parts)
    y <- y[, origin@order]

    y
  }
)

# MUTATORS
#' @include AllGenerics.R
NULL

# Not exported
get_transformation <- function(x) {
  if (methods::is(x, "LR")) return("Pairwise Log-Ratio")
  if (methods::is(x, "CLR")) return("Centered Log-Ratio")
  if (methods::is(x, "ALR")) return("Additive Log-Ratio")
  if (methods::is(x, "ILR")) return("Isometric Log-Ratio")
  if (methods::is(x, "PLR")) return("Pivot Log-Ratio")
}

# Predicates ===================================================================
#' @export
#' @rdname CompositionMatrix-class
is_composition <- function(object) {
  methods::is(object, "CompositionMatrix")
}

#' @export
#' @rdname LogRatio-class
is_logratio <- function(object) {
  methods::is(object, "LogRatio")
}

#' @export
#' @rdname ReferenceGroups-class
is_grouped <- function(object) {
  methods::is(object, "ReferenceGroups")
}

# Getter =======================================================================
#' @export
#' @method labels CompositionMatrix
labels.CompositionMatrix <- function(object, ...) {
  colnames(object)
}

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

#' @export
#' @method labels LogRatio
labels.LogRatio <- function(object, ...) {
  object@ratio
}

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

# Weights ======================================================================
#' @export
#' @method weights ALR
weights.ALR <- function(object, ...) {
  w <- object@weights
  w[-1] * w[1]
}

#' @export
#' @rdname mutators
#' @aliases weights,ALR-method
setMethod("weights", "ALR", weights.ALR)

#' @export
#' @method weights LR
weights.LR <- function(object, ...) {
  w <- object@weights
  w <- utils::combn(
    x = w,
    m = 2,
    FUN = function(x) Reduce(`*`, x),
    simplify = FALSE
  )
  unlist(w)
}

#' @export
#' @rdname mutators
#' @aliases weights,LR-method
setMethod("weights", "LR", weights.LR)

#' @export
#' @method weights LogRatio
weights.LogRatio <- function(object, ...) {
  object@weights
}

#' @export
#' @rdname mutators
#' @aliases weights,LogRatio-method
setMethod("weights", "LogRatio", weights.LogRatio)

# Totals =======================================================================
#' @export
#' @rdname totals
#' @aliases totals,CompositionMatrix-method
setMethod("totals", "CompositionMatrix", function(object) object@totals)

#' @export
#' @rdname totals
#' @aliases totals,LogRatio-method
setMethod("totals", "LogRatio", function(object) object@totals)

#' @export
#' @rdname totals
setMethod(
  f = "totals<-",
  signature = "CompositionMatrix",
  definition = function(object, value) {
    object@totals <- if (is.null(value)) rowSums(object) else as.numeric(value)
    methods::validObject(object)
    object
  }
)

# DATA SUMMARY: DESCRIBE
#' @include AllGenerics.R
NULL

#' @export
#' @rdname describe
#' @aliases describe,CompositionMatrix-method
setMethod(
  f = "describe",
  signature = c(x = "CompositionMatrix"),
  definition = function(x) {
    ## Dimensions
    m <- nrow(x)
    p <- ncol(x)
    rows <- sprintf(ngettext(m, "%d composition", "%d compositions"), m)
    cols <- sprintf(ngettext(p, "with %d part", "with %d parts"), p)
    msg_tbl <- sprintf("%s %s:", rows, cols)

    ## Message
    cat(msg_tbl)
    .describe(x)

    invisible(x)
  }
)

#' @export
#' @rdname describe
#' @aliases describe,GroupedComposition-method
setMethod(
  f = "describe",
  signature = c(x = "GroupedComposition"),
  definition = function(x) {
    ## Dimensions
    m <- nrow(x)
    p <- ncol(x)

    rows <- sprintf(ngettext(m, "%d composition", "%d compositions"), m)
    cols <- sprintf(ngettext(p, "with %d part", "with %d parts"), p)
    msg_tbl <- sprintf("%s %s:", rows, cols)

    ## Groups
    i <- group_length(x)
    ls_grp <- paste0(dQuote(group_levels(x)), collapse = ", ")
    msg_grp <- sprintf(ngettext(i, "%d group", "%d groups"), i)
    msg_grp <- sprintf("%s: %s.", msg_grp, ls_grp)

    j <- sum(!is_assigned(x))
    msg_ung <- sprintf(ngettext(j, "%d unassigned sample.", "%d unassigned samples."), j)

    cat(msg_tbl, msg_grp, msg_ung, sep = "\n* ")
    .describe(x)

    invisible(x)
  }
)

.describe <- function(x) {
  m <- nrow(x)
  p <- ncol(x)

  ## Missing values
  n_NA <- sum(count(x, f = is.na))
  m_NA <- sum(detect(x, f = is.na, margin = 1))
  p_NA <- sum(detect(x, f = is.na, margin = 2))
  pc <- label_percent(c(m_NA / m, p_NA / p), digits = 1, trim = TRUE)

  msg_NA <- sprintf(ngettext(n_NA, "\n%d missing value:", "\n%d missing values:"), n_NA)

  rows_NA <- ngettext(m_NA, "%d observation (%s) contains missing values.",
                      "%d observations (%s) contain missing values.")
  msg_row_NA <- sprintf(rows_NA, m_NA, pc[[1]])

  cols_NA <- ngettext(p_NA, "%d variable (%s) contains missing values.",
                      "%d variables (%s) contain missing values.")
  msg_col_NA <- sprintf(cols_NA, p_NA, pc[[2]])

  ## Constant columns
  p_var <- sum(detect(x, f = function(x) is_unique(x), margin = 2))
  cols_var <- ngettext(p_var, "%d part with no variance.",
                       "%d parts with no variance.")
  msg_col_var <- sprintf(cols_var, p_var)

  ## Sparsity
  spa <- sparsity(x, count = FALSE)
  msg_spa <- sprintf(tr_("%s of values are zero."), label_percent(spa, digits = 1))

  cat(msg_NA, msg_row_NA, msg_col_NA, sep = "\n* ")
  cat(tr_("\nData checking:"), msg_spa, msg_col_var, sep = "\n* ")
}

# DISTANCES
#' @include AllGenerics.R
NULL

# Distances ====================================================================
#' @export
#' @method dist CompositionMatrix
dist.CompositionMatrix <- function(x, method = "euclidean",
                                   diag = FALSE, upper = FALSE, p = 2) {
  x <- transform_clr(x)
  stats::dist(x, method = method, diag = diag, upper = upper, p = p)
}

#' @export
#' @rdname dist
#' @aliases dist,CompositionMatrix-method
setMethod("dist", "CompositionMatrix", dist.CompositionMatrix)

# Mahalanobis ==================================================================
#' @export
#' @method mahalanobis CompositionMatrix
mahalanobis.CompositionMatrix <- function(x, center, cov, ..., robust = TRUE,
                                          method = c("mve", "mcd")) {
  ## Transformation
  x <- transform_ilr(x)
  mahalanobis(x, center, cov, ..., robust = robust, method = method)
}

#' @export
#' @rdname mahalanobis
#' @aliases mahalanobis,CompositionMatrix-method
setMethod("mahalanobis", "CompositionMatrix", mahalanobis.CompositionMatrix)

#' @export
#' @method mahalanobis ILR
mahalanobis.ILR <- function(x, center, cov, ..., robust = TRUE,
                            method = c("mve", "mcd")) {

  if (missingORnull(center) | missingORnull(cov)) {
    if (!robust) method <- "classical" # Standard estimators
    else method <- match.arg(method, several.ok = FALSE) # Robust estimators
    v <- MASS::cov.rob(x, method = method, ...)
  }

  est <- list(center = NULL, cov = NULL)
  est$center <- if (missingORnull(center)) v$center else center
  est$cov <- if (missingORnull(cov)) v$cov else cov

  message(v$sing)

  stats::mahalanobis(x, center = est$center, cov = est$cov)
}

#' @export
#' @rdname mahalanobis
#' @aliases mahalanobis,ILR-method
setMethod("mahalanobis", "ILR", mahalanobis.ILR)

# BARPLOT
#' @include AllGenerics.R
NULL

# CompositionMatrix ============================================================
#' @export
#' @method barplot CompositionMatrix
barplot.CompositionMatrix <- function(height, ...,
                                      order_columns = FALSE, order_rows = NULL,
                                      decreasing = TRUE, names = TRUE,
                                      space = 0.2, offset = 0.025,
                                      palette_color = palette_color_discrete(),
                                      border = NA, axes = TRUE, legend = TRUE) {
  ## Validation
  if (ncol(height) < 2) {
    stop(tr_("At least two compositional parts are needed."), call. = FALSE)
  }

  ## Prepare data
  xy <- prepare_barplot(height, order_columns = order_columns,
                        order_rows = order_rows, decreasing = decreasing,
                        offset = offset)
  parts <- factor(xy$data$column, levels = colnames(height))
  col <- palette_color(parts)
  n <- nrow(height)

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

  ## Save and restore
  mfrow <- graphics::par("mfrow")
  mar <- graphics::par("mar")
  nlines <- height2line("M", cex = cex.axis)
  mar[1] <- 3
  mar[2] <- if (names) width2line(rownames(height), cex = cex.axis) + 0.1 else nlines
  mar[3] <- nlines
  mar[4] <- nlines

  old_par <- graphics::par(mfrow = mfrow, 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
  space <- 1 / n * 0.5 - space / n * 0.5
  xlim <- range(0, 1)
  ylim <- range(xy$data$y)
  ylim <- ylim + c(0, 2 * (offset + space) * legend) + c(-1, 1) * space
  graphics::plot.window(xlim = xlim, ylim = ylim, xaxs = "i", yaxs = "i")

  ## Plot
  graphics::rect(
    xleft = xy$data$xmin,
    ybottom = xy$data$y - space,
    xright = xy$data$xmax,
    ytop = xy$data$y + space,
    col = col,
    border = border
  )

  ## Construct axis
  if (axes) {
    at <- graphics::axTicks(side = 1)
    graphics::axis(side = 1, at = at, labels = label_percent(at),
                   xpd = NA, las = 1, cex.axis = cex.axis, col.axis = col.axis,
                   font.axis = font.axis)
    graphics::mtext(text = names(xy$groups), side = 4, line = 0, at = xy$groups,
                    cex = cex.axis, col = col.axis, font = font.axis)
    if (names) {
      graphics::axis(side = 2, at = unique(xy$data$y), labels = unique(xy$data$row),
                     las = 2, lty = 0, cex.axis = cex.axis, col.axis = col.axis,
                     font.axis = font.axis)
    }
  }

  ## Add legend
  if (legend) {
    graphics::rect(
      xleft = cumsum(xy$mean) - xy$mean,
      ybottom = max(xy$data$y) + offset + space,
      xright = cumsum(xy$mean),
      ytop = max(xy$data$y) + 2 * (offset + space),
      col = unique(col),
      border = border
    )
    lab <- label_chemical(names(xy$mean))
    graphics::mtext(text = parse(text = lab), side = 3, line = 0,
                    at = cumsum(xy$mean) - xy$mean / 2,
                    cex = cex.axis, col = unique(col), font = font.axis)
  }

  invisible(height)
}

#' @export
#' @rdname barplot
#' @aliases barplot,CompositionMatrix-method
setMethod("barplot", c(height = "CompositionMatrix"), barplot.CompositionMatrix)

prepare_barplot <- function(x, order_rows = NULL, order_columns = FALSE,
                            decreasing = TRUE, offset = 0.025) {
  ## Relative frequencies
  n <- nrow(x)
  x <- x / rowSums(x)

  ## Validation
  stopifnot(is_composition(x))

  ## Row order
  if (!is.null(order_rows)) {
    j <- x[, order_rows, drop = TRUE]
    i <- order(j, decreasing = decreasing)
    x <- x[i, , drop = FALSE]
  }

  ## Columns order
  center <- mean(x)
  if (order_columns) {
    col_order <- order(center, decreasing = FALSE)
    center <- center[col_order]
    x <- x[, col_order, drop = FALSE]
  }

  ## Grouping
  if (!is_grouped(x)) {
    x <- group(x, by = rep(NA, n), verbose = FALSE)
  }
  spl <- group_split(x)
  z <- do.call(rbind, spl)

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

  xmax <- t(apply(X = z, MARGIN = 1, FUN = cumsum))
  xmin <- xmax - z
  data$xmin <- as.vector(xmin)
  data$xmax <- as.vector(xmax)
  data$y <- as.vector(n + 1 - as.numeric(row)) / n # Reverse levels order

  ## Offset
  n_grp <- group_length(x)
  n_spl <- group_size(x)
  offset <- rev(seq_len(n_grp)) * offset - offset
  data$y <- data$y + rep(offset, n_spl)[as.numeric(row)]

  list(
    data = data,
    mean = center,
    groups = 1 - cumsum(n_spl) / n + n_spl / n * 0.5 + offset
  )
}

# OUTLIERS
#' @include AllGenerics.R
NULL

# Find =========================================================================
#' @export
#' @rdname detect_outlier
#' @aliases detect_outlier,CompositionMatrix,missing-method
setMethod(
  f = "detect_outlier",
  signature = c(object = "CompositionMatrix", reference = "missing"),
  definition = function(object, ..., robust = TRUE, method = c("mve", "mcd"),
                        quantile = 0.975) {
    methods::callGeneric(object, reference = object, robust = robust,
                         method = method, quantile = quantile)
  }
)

#' @export
#' @rdname detect_outlier
#' @aliases detect_outlier,CompositionMatrix,CompositionMatrix-method
setMethod(
  f = "detect_outlier",
  signature = c(object = "CompositionMatrix", reference = "CompositionMatrix"),
  definition = function(object, reference, ..., quantile = 0.975, robust = TRUE,
                        method = c("mve", "mcd")) {
    ## Validation
    if (!identical(colnames(object), colnames(reference))) {
      stop(tr_("Column names do not match!"), call. = FALSE)
    }

    ## Transformation
    obj <- transform_ilr(object)
    ref <- transform_ilr(reference)

    ## Clean
    n <- nrow(ref)
    p <- ncol(ref)
    if (n < (p + 1)) {
      msg <- tr_("Sample size is too small (%d).")
      stop(sprintf(msg, n), call. = FALSE)
    }
    if (n < (2 * p)) {
      msg <- tr_("Possibly too small sample size (%d).")
      warning(sprintf(msg, n), call. = FALSE)
    }

    ## Compute center and spread + Mahalanobis distance
    ## Standard estimators
    estc <- list(center = colMeans(ref, na.rm = TRUE), cov = cov(ref))
    dc <- stats::mahalanobis(obj, center = estc$center, cov = estc$cov)

    ## Robust estimators
    dr <- rep(NA_real_, nrow(obj))
    if (robust) {
      method <- match.arg(method, several.ok = FALSE)
      estr <- MASS::cov.rob(ref, method = method, ...)
      dr <- stats::mahalanobis(obj, center = estr$center, cov = estr$cov)
    }

    ## Threshold
    limit <- stats::qchisq(p = quantile, df = p)

    .OutlierIndex(
      samples = rownames(obj),
      standard = sqrt(dc),
      robust = sqrt(dr),
      limit = sqrt(limit),
      dof = p
    )
  }
)

#' @export
#' @rdname detect_outlier
#' @aliases is_outlier,OutlierIndex-method
setMethod(
  f = "is_outlier",
  signature = c("OutlierIndex"),
  definition = function(object, robust = TRUE) {
    d <- if (robust) object@robust else object@standard
    out <- d > object@limit
    names(out) <- object@samples
    out
  }
)

# Plot =========================================================================
#' @export
#' @method plot OutlierIndex
plot.OutlierIndex <- function(x, ...,
                              type = c("dotchart", "distance"),
                              robust = TRUE,
                              symbols = c(16, 1, 3),
                              xlim = NULL, ylim = 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 = "topleft")) {
  ## Get data
  dc <- x@standard
  dr <- x@robust
  dof <- x@dof
  limit <- x@limit
  n <- length(dc)

  ## Validation
  if (all(is.na(dr))) {
    robust <- FALSE
    type <- "dotchart"
  }
  type <- match.arg(type, several.ok = FALSE)

  ## Graphical parameters
  shape <- rep(symbols[[1L]], n)
  if (robust || type == "distance") shape[dr > limit] <- symbols[[3L]]
  if (!robust || type == "distance") shape[dc > limit] <- symbols[[2L]]

  cy <- if (robust) dr else dc
  dlab <- ifelse(robust, tr_("Robust Mahalanobis distance"),
                 tr_("Standard Mahalanobis distance"))
  ylab <- ylab %||% dlab

  if (type == "dotchart") {
    asp <- NA
    cx <- seq_along(dc)
    xlab <- xlab %||% tr_("Index")
    panel <- function() {
      graphics::points(x = cx, y = cy, pch = shape, ...)
      graphics::abline(h = limit, lty = 1)
    }
  }
  if (type == "distance") {
    asp <- 1
    cx <- dc
    cy <- dr
    xlab <- xlab %||% tr_("Standard Mahalanobis distance")
    ylab <- ylab %||% tr_("Robust Mahalanobis distance")
    panel <- function() {
      graphics::points(x = cx, y = cy, pch = shape, ...)
      graphics::abline(h = limit, lty = 1)
      graphics::abline(v = limit, lty = 1)
      graphics::abline(a = 0, b = 1, lty = 2, col = "darkgrey")
    }
  }

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

  ## Set plotting coordinates
  xlim <- xlim %||% range(cx, finite = TRUE)
  ylim <- ylim %||% range(cy, finite = TRUE)
  graphics::plot.window(xlim = xlim, ylim = ylim, asp = asp)

  ## Evaluate pre-plot expressions
  panel.first

  ## Plot
  panel()

  ## 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) {
    graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab)
  }

  ## Add legend
  if (is.list(legend)) {
    if (type == "distance") {
      lab <- c(tr_("No outlier"), tr_("Robust only"), tr_("Both"))
      pch <- symbols
    } else {
      lab <- c(tr_("No outlier"), tr_("Outlier"))
      pch <- symbols[-2 - !robust]
    }
    args <- list(x = "topleft", legend = lab, pch = pch, bty = "n", xpd = NA)
    args <- utils::modifyList(args, legend)
    do.call(graphics::legend, args = args)
  }

  invisible(x)
}

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

# BIND
#' @include AllGenerics.R
NULL

#' @export
#' @rdname bind
#' @aliases rbind2,CompositionMatrix,CompositionMatrix-method
setMethod(
  f = "rbind2",
  signature = c(x = "CompositionMatrix", y = "CompositionMatrix"),
  definition = function(x, y) {
    mtx_x <- methods::as(x, "matrix")
    mtx_y <- methods::as(y, "matrix")

    spl <- c(rownames(x), rownames(y))
    if (any(duplicated(spl))) {
      warning(tr_("Duplicated rownames!"), call. = FALSE)
      spl <- make.unique(spl, sep = "_")
    }

    z <- rbind(mtx_x, mtx_y)
    rownames(z) <- spl

    .CompositionMatrix(z, totals = c(totals(x), totals(y)))
  }
)

#' @export
#' @rdname bind
#' @aliases rbind2,GroupedComposition,GroupedComposition-method
setMethod(
  f = "rbind2",
  signature = c(x = "GroupedComposition", y = "GroupedComposition"),
  definition = function(x, y) {
    stop(tr_("Combining grouped matrices by rows is not currently supported."),
         call. = FALSE)
  }
)

# REPLACE
#' @include AllGenerics.R
NULL

# Zeros ========================================================================
#' @export
#' @rdname replace_zero
#' @aliases replace_zero,CompositionMatrix-method
setMethod(
  f = "replace_zero",
  signature = c(x = "CompositionMatrix"),
  definition = function(x, value, delta = 2/3) {
    ## Validation
    D <- ncol(x)
    if (is.null(value)) return(x)
    if (length(value) == 1) value <- rep(value, D)
    if (length(value) > 1) arkhe::assert_length(value, D)
    if (length(delta) > 1) arkhe::assert_length(delta, D)

    sigma <- value * delta
    r <- apply(X = x, MARGIN = 1, FUN = zero_multiplicative, sigma = sigma)

    methods::initialize(x, t(r))
  }
)

# zero_additive <- function(x, sigma) {
#   D <- length(x)
#
#   is_zero <- x == 0
#   Z <- sum(is_zero)
#
#   x[is_zero] <- (sigma * (Z + 1) * (D - Z)) / D^2
#   x[!is_zero] <- x[!is_zero] - (sigma * (Z + 1) * Z) / D^2
#
#   x
# }
zero_multiplicative <- function(x, sigma) {
  is_zero <- x == 0 & !is.na(x)
  x[is_zero] <- sigma[is_zero]
  x[!is_zero] <- x[!is_zero] * (1 - sum(sigma[is_zero]) / 1)
  x
}

# Missing values ===============================================================
#' @export
#' @rdname replace_NA
#' @aliases replace_NA,CompositionMatrix-method
setMethod(
  f = "replace_NA",
  signature = c(x = "CompositionMatrix"),
  definition = function(x, value) {
    ## Validation
    D <- ncol(x)
    if (is.null(value)) return(x)
    if (length(value) == 1) value <- rep(value, D)
    if (length(value) > 1) arkhe::assert_length(value, D)

    r <- apply(X = x, MARGIN = 1, FUN = missing_multiplicative, sigma = value)

    methods::initialize(x, t(r))
  }
)

missing_multiplicative <- function(x, sigma) {
  is_missing <- is.na(x)
  x[is_missing] <- sigma[is_missing]
  x[!is_missing] <- x[!is_missing] * (1 - sum(sigma[is_missing])) / sum(x[!is_missing])
  x
}

# MIXED-MODE ANALYSIS
#' @include AllGenerics.R
NULL

# First approach ===============================================================
#' @export
#' @describeIn mix First approach of mixed-mode analysis.
#' @aliases mix,matrix,matrix-method
setMethod(
  f = "mix",
  signature = c(x = "matrix", y = "matrix"),
  definition = function(x, y, lambda = 1, ...) {
    ## Validation
    arkhe::assert_package("cluster")
    stopifnot(nrow(y) == nrow(x))

    X <- vector(mode = "list", length = lambda + 1)
    X[[1]] <- y
    lambda <- as.integer(lambda[[1L]])
    for (i in seq_len(lambda)) X[[i + 1]] <- x
    X <- do.call(cbind, X)

    d <- cluster::daisy(X, metric = "gower", ...)
    as.dist(d)
  }
)

# Second approach ==============================================================
#' @export
#' @describeIn mix Second approach of mixed-mode analysis.
#' @aliases mix,dist,dist-method
setMethod(
  f = "mix",
  signature = c(x = "dist", y = "dist"),
  definition = function(x, y, mu = 0.5) {
    ## Validation
    stopifnot(mu >= 0 & mu <= 1)

    x <- as.matrix(x)
    y <- as.matrix(y)

    d <- mu * x + (1 - mu) * y
    as.dist(d)
  }
)

# DATA TRANSFORMATION: ADDITIVE LOG RATIO
#' @include AllGenerics.R
NULL

# ALR ==========================================================================
alr_base <- function(D) {
  V <- diag(1, nrow = D, ncol = D - 1)
  V[D, ] <- -1

  V
}

#' @export
#' @rdname transform_alr
#' @aliases transform_alr,CompositionMatrix-method
setMethod(
  f = "transform_alr",
  signature = c(object = "CompositionMatrix"),
  definition = function(object, j = ncol(object), weights = FALSE) {
    D <- ncol(object)
    parts <- colnames(object)

    ## Reorder (move denominator)
    j <- if (is.character(j)) which(parts == j) else as.integer(j)
    ordering <- c(which(j != seq_len(D)), j)
    parts <- parts[ordering]
    z <- object[, ordering, drop = FALSE]

    ## Compute ratios
    base <- alr_base(D)
    alr <- log(z, base = exp(1)) %*% base
    rownames(alr) <- rownames(object)
    colnames(alr) <- paste(parts[-D], parts[D], sep = "_")

    ## Compute weights
    weights <- make_weights(object, weights = weights)

    .ALR(
      alr,
      parts = parts,
      ratio = paste(parts[-D], parts[D], sep = "/"),
      order = order(ordering),
      base = base,
      weights = weights,
      totals = totals(object)
    )
  }
)

#' @export
#' @rdname transform_alr
#' @aliases transform_alr,GroupedComposition-method
setMethod(
  f = "transform_alr",
  signature = c(object = "GroupedComposition"),
  definition = function(object, j = ncol(object), weights = FALSE) {
    z <- methods::callNextMethod()
    .GroupedALR(z, group_indices = group_indices(object),
                group_levels = group_levels(object))
  }
)

#' @export
#' @rdname transform_alr
#' @aliases transform_alr,CLR-method
setMethod(
  f = "transform_alr",
  signature = c(object = "CLR"),
  definition = function(object, j = ncol(object)) {
    D <- ncol(object)
    parts <- object@parts

    ## Reorder (move denominator)
    j <- if (is.character(j)) which(parts == j) else as.integer(j)
    ordering <- c(which(j != seq_len(D)), j)
    parts <- parts[ordering]
    z <- object[, ordering, drop = FALSE]

    ## Compute ratios
    base <- alr_base(D)
    alr <- z %*% base
    rownames(alr) <- rownames(object)
    colnames(alr) <- paste(parts[-D], parts[D], sep = "_")

    .ALR(
      alr,
      parts = parts,
      ratio = paste(parts[-D], parts[D], sep = "/"),
      order = order(ordering),
      base = base,
      weights = object@weights,
      totals = totals(object)
    )
  }
)

#' @export
#' @rdname transform_alr
#' @aliases transform_alr,GroupedCLR-method
setMethod(
  f = "transform_alr",
  signature = c(object = "GroupedCLR"),
  definition = function(object, j = ncol(object), weights = FALSE) {
    z <- methods::callNextMethod()
    .GroupedALR(z, group_indices = group_indices(object),
                group_levels = group_levels(object))
  }
)

# DATA TRANSFORMATION: CENTERED LOG RATIO
#' @include AllGenerics.R
NULL

# CLR ==========================================================================
clr_base <- function(D, weights = rep(1 / D, D)) {
  diag(D) - matrix(data = weights, nrow = D, ncol = D)
}

#' @export
#' @rdname transform_clr
#' @aliases transform_clr,CompositionMatrix-method
setMethod(
  f = "transform_clr",
  signature = c(object = "CompositionMatrix"),
  definition = function(object, weights = FALSE) {
    D <- ncol(object)
    parts <- colnames(object)

    weights <- make_weights(object, weights = weights)
    base <- clr_base(D, weights = weights)
    clr <- log(object, base = exp(1)) %*% base
    dimnames(clr) <- dimnames(object)

    .CLR(
      clr,
      parts = parts,
      ratio = parts,
      order = seq_len(D),
      base = base,
      weights = weights,
      totals = totals(object)
    )
  }
)

#' @export
#' @rdname transform_clr
#' @aliases transform_clr,GroupedComposition-method
setMethod(
  f = "transform_clr",
  signature = c(object = "GroupedComposition"),
  definition = function(object, weights = FALSE) {
    z <- methods::callNextMethod()
    .GroupedCLR(z, group_indices = group_indices(object),
                group_levels = group_levels(object))
  }
)

#' @export
#' @rdname transform_clr
#' @aliases transform_clr,ALR-method
setMethod(
  f = "transform_clr",
  signature = c(object = "ALR"),
  definition = function(object) {
    D <- ncol(object) + 1
    w <- object@weights

    base <- clr_base(D, weights = w)
    clr <- object %*% base[-D, ]
    dimnames(clr) <- list(rownames(object), object@parts)

    .CLR(
      clr,
      parts = object@parts,
      ratio = object@parts,
      order = seq_len(D),
      base = base,
      weights = w,
      totals = totals(object)
    )
  }
)

#' @export
#' @rdname transform_clr
#' @aliases transform_clr,GroupedALR-method
setMethod(
  f = "transform_clr",
  signature = c(object = "GroupedALR"),
  definition = function(object) {
    z <- methods::callNextMethod()
    .GroupedCLR(z, group_indices = group_indices(object),
                group_levels = group_levels(object))
  }
)

# DATA TRANSFORMATION: LOG RATIO
#' @include AllGenerics.R
NULL

# LR ===========================================================================
#' @export
#' @rdname transform_lr
#' @aliases transform_lr,CompositionMatrix-method
setMethod(
  f = "transform_lr",
  signature = c(object = "CompositionMatrix"),
  definition = function(object, weights = FALSE) {
    J <- ncol(object)
    parts <- colnames(object)

    ## Compute weights
    weights <- make_weights(object, weights = weights)

    ## Computes ratios
    jj <- utils::combn(seq_len(J), 2, simplify = FALSE)
    lr <- matrix(data = 0, nrow = nrow(object), ncol = length(jj))
    for (i in seq_along(jj)) {
      a <- jj[[i]][[1]]
      b <- jj[[i]][[2]]
      r <- object[, a, drop = TRUE] / object[, b, drop = TRUE]
      lr[, i] <- log(r, base = exp(1))
    }

    ## Make names
    ratio <- unlist(utils::combn(parts, 2, FUN = paste, collapse = "/", simplify = FALSE))
    rownames(lr) <- rownames(object)
    colnames(lr) <- ratio

    .LR(
      lr,
      parts = parts,
      ratio = ratio,
      order = seq_len(J),
      weights = weights,
      totals = totals(object)
    )
  }
)

#' @export
#' @rdname transform_lr
#' @aliases transform_lr,GroupedComposition-method
setMethod(
  f = "transform_lr",
  signature = c(object = "GroupedComposition"),
  definition = function(object, weights = FALSE) {
    z <- methods::callNextMethod()
    .GroupedLR(z, group_indices = group_indices(object),
               group_levels = group_levels(object))
  }
)

# COMPOSITIONAL MEAN
#' @include AllGenerics.R
NULL

#' @export
#' @rdname condense
#' @aliases condense,CompositionMatrix-method
setMethod(
  f = "condense",
  signature = "CompositionMatrix",
  definition = function(x, by, ignore_na = FALSE, ignore_zero = TRUE,
                        verbose = FALSE, ...) {
    x <- group(x, by = by)
    y <- methods::callGeneric(x = x, ignore_na = ignore_na,
                              ignore_zero = ignore_zero, verbose = verbose, ...)
    ungroup(y)
  }
)

#' @export
#' @rdname condense
#' @aliases condense,GroupedComposition-method
setMethod(
  f = "condense",
  signature = "GroupedComposition",
  definition = function(x, by = NULL, ignore_na = FALSE, ignore_zero = TRUE,
                        verbose = FALSE, ...) {
    ## Grouping
    grp <- group_factor(x)
    if (!is.null(by)) x <- group(x, by = by, verbose = verbose)

    ## Compute mean
    z <- aggregate(
      x = x,
      FUN = mean,
      ignore_na = ignore_na,
      ignore_zero = ignore_zero,
      simplify = TRUE
    )
    tot <- tapply(X = totals(x), INDEX = group_factor(x), FUN = mean)

    z <- .CompositionMatrix(z, totals = as.numeric(tot))
    group(z, by = flatten_chr(x = grp, by = group_factor(x)), verbose = verbose)
  }
)

flatten_chr <- function(x, by) {
  x <- as.character(x)
  z <- tapply(X = x, INDEX = by, FUN = unique, simplify = FALSE)
  z <- vapply(X = z, FUN = paste0, FUN.VALUE = character(1), collapse = ":")
  z
}

# HELPERS

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

missingORnull <- function(x) {
  missing(x) || is.null(x)
}

has_rownames <- function(x) {
  .row_names_info(x, type = 1L) > 0L &&
    !is.na(.row_names_info(x, type = 0L)[[1L]])
}

make_names <- function(x, n = length(x), prefix = "X") {
  x <- if (n > 0) x %||% paste0(prefix, seq_len(n)) else character(0)
  x <- make.unique(x, sep = "_")
  x
}

#' 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
width2line <- function(x, ...) {
  (max(graphics::strwidth(x, units = "inch", ...)) /
     graphics::par("cin")[2] + graphics::par("mgp")[2]) * graphics::par("cex")
}
height2line <- function(x, ...) {
  (max(graphics::strheight(x, units = "inch", ...)) /
     graphics::par("cin")[2] + graphics::par("mgp")[2]) * graphics::par("cex")
}

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

#' Label Chemical Formula
#'
#' @param x A [`character`] vector.
#' @return A [`character`] vector.
#' @keywords internal
#' @noRd
label_chemical <- function(x, digits = NULL, trim = FALSE) {
  if (!all(is_chemical(x))) return(x)
  x <- gsub(pattern = "([[:digit:]]+)", replacement = "[\\1]", x = x, fixed = FALSE)
  x <- gsub(pattern = "([[:alpha:]]+)", replacement = "\"\\1\"", x = x, fixed = FALSE)
  gsub(pattern = "]\"", replacement = "]*\"", x = x, fixed = FALSE)
}

#' Column Weights
#'
#' Computes column weights.
#' @param x A `numeric` [`matrix`].
#' @param weights A [`logical`] scalar: should varying weights (column means)
#'  be computed? If `FALSE` (the default), equally-weighted parts are used.
#'  Alternatively, a positive [`numeric`] vector of weights can be specified
#'  (will be rescaled to sum to \eqn{1}).
#' @return A [`numeric`] vector.
#' @keywords internal
#' @noRd
make_weights <- function(x, weights = FALSE) {
  D <- ncol(x)

  w <- if (isTRUE(weights)) colMeans(x) else rep(1 / D, D)
  if (is.numeric(weights)) {
    arkhe::assert_length(weights, D)
    arkhe::assert_positive(weights, strict = TRUE)
    w <- weights / sum(weights) # Sum up to 1
  }

  unname(w)
}

# PCA
#' @include AllGenerics.R
NULL

#' @export
#' @describeIn pca PCA of centered log-ratio, i.e. log-ratio analysis (LRA).
#' @aliases pca,CompositionMatrix-method
setMethod(
  f = "pca",
  signature = c("CompositionMatrix"),
  definition = function(object, center = TRUE, scale = FALSE, rank = NULL,
                        sup_row = NULL, sup_col = NULL,
                        weight_row = NULL, weight_col = NULL) {
    message(tr_("PCA of centered log-ratio (CLR)."))
    x <- transform_clr(object)
    methods::callGeneric(object = x, center = center, scale = scale,
                         rank = rank, sup_row = sup_row, sup_col = sup_col,
                         weight_row = weight_row, weight_col = weight_col)
  }
)

#' @export
#' @rdname pca
#' @aliases pca,LogRatio-method
setMethod(
  f = "pca",
  signature = c("LogRatio"),
  definition = function(object, center = TRUE, scale = FALSE, rank = NULL,
                        sup_row = NULL, sup_col = NULL,
                        weight_row = NULL, weight_col = NULL) {
    methods::callNextMethod(object, center = center, scale = scale,
                            rank = rank, sup_row = sup_row,
                            sup_col = sup_col, weight_row = weight_row,
                            weight_col = weight_col)
  }
)

# MARGIN
#' @include AllGenerics.R
NULL

# Margin =======================================================================
#' @export
#' @rdname margin
#' @aliases margin,CompositionMatrix-method
setMethod(
  f = "margin",
  signature = c("CompositionMatrix"),
  definition = function(x, parts = c(1, 2), name = "*") {
    ## Validation
    p <- NCOL(x)
    parts <- unique(parts)
    if (is.character(parts)) parts <- match(parts, colnames(x))
    if (p <= length(parts)) return(x)
    if (p == length(parts) + 1) {
      d <- seq_len(p)
      d <- c(d[parts], d[-parts])
      return(x[, d, drop = FALSE])
    }

    rest <- x[, -parts, drop = FALSE]
    star <- apply(X = rest, MARGIN = 1, FUN = gmean)
    mar <- cbind(x[, parts, drop = FALSE], star)
    colnames(mar) <- c(colnames(x)[parts], name[[1L]])

    clo <- closure(mar)
    methods::initialize(x, clo)
  }
)

# GENERIC METHODS
#' @include AllClasses.R
NULL

# S4 dispatch to base S3 generic ===============================================
setGeneric("dist", package = "stats")
setGeneric("mahalanobis", package = "stats")

# Import S4 generics ===========================================================
#' @importMethodsFrom arkhe describe
#' @importMethodsFrom arkhe replace_NA
#' @importMethodsFrom arkhe replace_zero
#' @importMethodsFrom dimensio pca
NULL

# CoDa =========================================================================
#' Coerce to a Closed Compositional Matrix
#'
#' Coerces an object to a `CompositionMatrix` object.
#' @param from A [`matrix`] or [`data.frame`] to be coerced.
#' @param parts A `vector` giving the index of the column to be used a
#'  compositional parts. If `NULL` and `autodetect` is `TRUE` (the default),
#'  all `numeric` columns will be used.
#' @param groups An [`integer`] giving the index of the column to be used to
#'  group the samples. If `NULL` (the default), no grouping is stored.
#' @param autodetect A [`logical`] scalar: should `numeric` variables be
#'  automatically used as compositional parts?
#' @param verbose A [`logical`] scalar: should \R report extra information
#'  on progress?
#' @param ... Currently not used.
#' @details
#'  See `vignette("nexus")`.
#' @return A [`CompositionMatrix-class`] object.
#' @example inst/examples/ex-coerce.R
#' @author N. Frerebeau
#' @docType methods
#' @family compositional data tools
#' @aliases as_composition-method
setGeneric(
  name = "as_composition",
  def = function(from, ...) standardGeneric("as_composition"),
  valueClass = "CompositionMatrix"
)

#' Coerce to Amounts
#'
#' @param from A [`CompositionMatrix-class`] object.
#' @param ... Currently not used.
#' @return A [`numeric`] [`matrix`].
#' @example inst/examples/ex-coerce.R
#' @author N. Frerebeau
#' @docType methods
#' @family compositional data tools
#' @aliases as_amounts-method
setGeneric(
  name = "as_amounts",
  def = function(from, ...) standardGeneric("as_amounts"),
  valueClass = "matrix"
)

#' Data Description
#'
#' Describes an object.
#' @param x A [`CompositionMatrix-class`] object.
#' @return
#'  `describe()` is called for its side-effects. Invisibly returns `x`.
#' @example inst/examples/ex-describe.R
#' @author N. Frerebeau
#' @docType methods
#' @family data summaries
#' @name describe
#' @rdname describe
NULL

# Simplex ======================================================================
#' Operations in the Simplex
#'
#' Operators performing operations in the simplex.
#' @param x A [`CompositionMatrix-class`] object.
#' @param y A [`CompositionMatrix-class`] object or a [`numeric`] vector.
#' @details
#'  \describe{
#'   \item{`%perturbe%`}{[Perturbation operation][perturbation()].}
#'   \item{`%power%`}{[Powering operation][powering()].}
#'  }
#' @return
#'  A [`CompositionMatrix-class`] object or a [`numeric`] vector (same as `x`).
#' @example inst/examples/ex-arith.R
#' @author N. Frerebeau
#' @docType methods
#' @family operations in the simplex
#' @name arithmetic
#' @rdname arithmetic
NULL

#' @rdname arithmetic
#' @aliases `%perturbe%`-method
setGeneric(
  name = "%perturbe%",
  def = function(x, y) standardGeneric("%perturbe%")
)

#' @rdname arithmetic
#' @aliases `%power%`-method
setGeneric(
  name = "%power%",
  def = function(x, y) standardGeneric("%power%")
)

#' Closure Operation
#'
#' Closes compositions to sum to 1.
#' @param x A [`numeric`] vector or matrix.
#' @param total A [numeric] vector specifying the total amount to which the
#'  compositions should be closed (defaults to 1).
#' @param na.rm A [`logical`] scalar: should missing values be removed?
#' @param ... Currently not used.
#' @return A [`numeric`] vector or matrix (same as `x`).
#' @example inst/examples/ex-arith.R
#' @author N. Frerebeau
#' @docType methods
#' @family operations in the simplex
#' @aliases closure-method
setGeneric(
  name = "closure",
  def = function(x, ...) standardGeneric("closure")
)

#' Perturbation Operation
#'
#' Perturbation of two compositions.
#' @param x,y A [`numeric`] vector of compositional data or a
#'  [`CompositionMatrix-class`] object.
#' @param ... Currently not used.
#' @details
#'  In compositional geometry, perturbation plays the role of sum (translation).
#'  It is the closed component-wise product of two compositions.
#' @return A [`numeric`] vector.
#' @example inst/examples/ex-arith.R
#' @author N. Frerebeau
#' @docType methods
#' @family operations in the simplex
#' @aliases perturbation-method
setGeneric(
  name = "perturbation",
  def = function(x, y, ...) standardGeneric("perturbation")
)

#' Powering Operation
#'
#' Perturbation of two compositions.
#' @param x A [`numeric`] vector of compositional data or a
#'  [`CompositionMatrix-class`] object.
#' @param a A [`numeric`] constant.
#' @param ... Currently not used.
#' @details
#'  In compositional geometry, powering replaces the product of a vector by a
#'  scalar (scaling) and is defined as the closed powering of the components by
#'  a given scalar.
#' @return A [`numeric`] vector.
#' @example inst/examples/ex-arith.R
#' @author N. Frerebeau
#' @docType methods
#' @family operations in the simplex
#' @aliases powering-method
setGeneric(
  name = "powering",
  def = function(x, a, ...) standardGeneric("powering")
)

#' Scalar Product
#'
#' Computes the Aitchison scalar product of two compositions.
#' @param x,y A [`CompositionMatrix-class`] object.
#' @param ... Currently not used.
#' @return A [`numeric`] vector.
#' @example inst/examples/ex-arith.R
#' @author N. Frerebeau
#' @docType methods
#' @family operations in the simplex
#' @aliases scalar-method
setGeneric(
  name = "scalar",
  def = function(x, y, ...) standardGeneric("scalar")
)

# Extract ======================================================================
## Subset ----------------------------------------------------------------------
#' Extract or Replace Parts of an Object
#'
#' Operators acting on objects to extract or replace parts.
#' @param x An object from which to extract element(s) or in which to replace
#'  element(s).
#' @param i,j Indices specifying elements to extract or replace. Indices are
#'  [`numeric`], [`integer`] or [`character`] vectors or empty (missing) or
#'  `NULL`. Numeric values are coerced to [`integer`] as by [as.integer()].
#'  Character vectors will be matched to the name of the elements.
#'  An empty index (a comma separated blank) indicates that all entries in that
#'  dimension are selected.
#' @param value A possible value for the element(s) of `x`.
#' @param drop A [`logical`] scalar: should the result be coerced to
#'  the lowest possible dimension? This only works for extracting elements,
#'  not for the replacement. Defaults to `FALSE`.
#' @param ... Currently not used.
# @section Subcomposition:
#  If `drop` is `FALSE`, subsetting some of the possible components of a
#  [`CompositionMatrix-class`] object will produce a closed *subcomposition*
#  (see examples).
#' @return
#'  A subsetted object of the same sort as `x`.
#' @example inst/examples/ex-subset.R
#' @author N. Frerebeau
#' @docType methods
#' @family subsetting methods
#' @name subset
#' @rdname subset
NULL

#' Combine Two Composition Matrices
#'
#' @param x,y A [`CompositionMatrix-class`] object.
#' @details
#'  `rbind2()` combine by rows.
#' @return
#'  A [`CompositionMatrix-class`] objects.
#' @example inst/examples/ex-split.R
#' @author N. Frerebeau
#' @docType methods
#' @family subsetting methods
#' @name bind
#' @rdname bind
NULL

## Mutators --------------------------------------------------------------------
#' Get or Set Parts of an Object
#'
#' Getters and setters to retrieve or set parts of an object.
#' @param object An object from which to get or set element(s).
# @param value A possible value for the element(s) of `x`.
#' @param ... Currently not used.
# @example inst/examples/ex-mutators.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name mutators
#' @rdname mutators
#' @aliases get set
NULL

#' Row Sums
#'
#' Retrieves or defines the row sums (before [closure][closure()]).
#' @param object An object from which to get or set `totals`.
#' @param value A possible value for the `totals` of `x`.
#' @return
#'  * `totals() <- value` returns an object of the same sort as `x` with the new
#'    row sums assigned.
#'  * `totals()` returns the row sums of `x`.
#' @example inst/examples/ex-coerce.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @aliases totals-method
setGeneric(
  name = "totals",
  def = function(object) standardGeneric("totals")
)

#' @rdname totals
setGeneric(
  name = "totals<-",
  def = function(object, value) standardGeneric("totals<-")
)

#' Coerce to a Data Frame
#'
#' @param x An \R object (typically, a [`CompositionMatrix-class`] object).
#' @param row.names A [`character`] vector giving the row names for the data
#'  frame, or `NULL`.
#' @param optional A [`logical`] scalar: should the names of the variables in
#'  the data frame be checked? If `FALSE` then the names of the variables in the
#'  data frame are checked to ensure that they are syntactically valid variable
#'  names and are not duplicated.
#' @param group_var A [`character`] string specifying the name of the column to
#'  create for group attribution.
#' @param group_after An [`integer`] specifying a subscript, after which the new
#'  `group_var` column is to be appended.
#' @param ... Currently not used.
#' @return
#'  A [`data.frame`].
#' @example inst/examples/ex-coerce.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name as.data.frame
#' @rdname as.data.frame
NULL

#' Matrix Transpose
#'
#' @param x A [`CompositionMatrix-class`] object.
#' @return
#'  A `matrix`, with dim and dimnames constructed appropriately from those of `x`.
#' @note
#'  Implemented only to ensure that `t()` always returns a base `matrix`.
#' @example inst/examples/ex-subset.R
#' @author N. Frerebeau
#' @docType methods
# @family mutators
#' @keywords internal
#' @name t
#' @rdname t
NULL

# Groups =======================================================================
#' Working With Groups
#'
#' Define or remove the (reference) groups to which the observations belong.
#' @param object An \R object (typically, a [`CompositionMatrix-class`] object).
#' @param by A possible value for the groups of `object` (typically, a
#'  [`character`] vector). If `value` is a [`list`],
#'  [`interaction(by)`][interaction()] defines the grouping.
#' @param add A [`logical`] scalar. If `TRUE`, add to existing groups.
#' @param verbose A [`logical`] scalar: should \R report extra information
#'  on progress?
#' @param ... Further parameters to be passed to internal methods.
#' @details
#'  Missing values (`NA`) can be used to specify that a sample does not belong
#'  to any group.
#' @return
#'  * `group()` returns a grouped object of the same sort as `object`.
#'  * `ungroup()` returns an ungrouped object of the same sort as `object`.
#' @example inst/examples/ex-group.R
#' @author N. Frerebeau
#' @docType methods
#' @family grouping methods
#' @aliases group-method
setGeneric(
  name = "group",
  def = function(object, ...) standardGeneric("group")
)

#' @rdname group
setGeneric(
  name = "ungroup",
  def = function(object, ...) standardGeneric("ungroup")
)

#' Grouping Metadata
#'
#' Retrieve the (reference) groups to which the observations belong.
#' @param object A [grouped][group()] \R object.
#' @example inst/examples/ex-group.R
#' @author N. Frerebeau
#' @docType methods
#' @family grouping methods
#' @name group_metadata
#' @rdname group_metadata
NULL

#' @rdname group_metadata
#' @aliases group_levels-method
setGeneric(
  name = "group_levels",
  def = function(object) standardGeneric("group_levels")
)

#' @rdname group_metadata
#' @aliases group_names-method
setGeneric(
  name = "group_names",
  def = function(object) standardGeneric("group_names")
)

#' @rdname group_metadata
#' @aliases group_rows-method
setGeneric(
  name = "group_rows",
  def = function(object) standardGeneric("group_rows")
)

#' @rdname group_metadata
#' @aliases group_length-method
setGeneric(
  name = "group_length",
  def = function(object) standardGeneric("group_length")
)

#' @rdname group_metadata
#' @aliases group_size-method
setGeneric(
  name = "group_size",
  def = function(object) standardGeneric("group_size")
)

#' @rdname group_metadata
#' @aliases group_indices-method
setGeneric(
  name = "group_indices",
  def = function(object) standardGeneric("group_indices")
)

#' @rdname group_metadata
#' @aliases is_assigned-method
setGeneric(
  name = "is_assigned",
  def = function(object) standardGeneric("is_assigned")
)

#' @rdname group_metadata
#' @aliases any_assigned-method
setGeneric(
  name = "any_assigned",
  def = function(object) standardGeneric("any_assigned")
)

#' @rdname group_metadata
#' @aliases any_assigned-method
setGeneric(
  name = "all_assigned",
  def = function(object) standardGeneric("all_assigned")
)

#' Divide into Groups
#'
#' Divides a compositional matrix by groups.
#' @param object,x A [`CompositionMatrix-class`] object.
#' @param f A `factor` in the sense that [`as.factor(f)`][as.factor()] defines
#'  the grouping, or a list of such factors in which case their interaction is
#'  used for the grouping (see [base::split()]).
#' @param drop A [`logical`] scalar: should levels that do not occur be dropped?
#' @param ... Currently not used.
#' @return
#'  A `list` of [`CompositionMatrix-class`] objects.
#' @example inst/examples/ex-split.R
#' @author N. Frerebeau
#' @docType methods
#' @family grouping methods
#' @aliases group_split-method
setGeneric(
  name = "group_split",
  def = function(object, ...) standardGeneric("group_split")
)

#' Group-based Subset
#'
#' @param object A [`GroupedComposition-class`] object.
#' @param which A [`character`] vector specifying the [groups][group()] of
#'  `object` to extract.
#' @param ... Currently not used.
#' @return
#'  A [`CompositionMatrix-class`] object.
#' @example inst/examples/ex-group.R
#' @author N. Frerebeau
#' @docType methods
#' @family grouping methods
#' @aliases group_subset-method
setGeneric(
  name = "group_subset",
  def = function(object, ...) standardGeneric("group_subset")
)

# Tools ========================================================================
#' Chemical Elements and Oxides
#'
#' Identify oxides and major, minor and traces elements in a compositional data
#' matrix.
#' @param object A [`CompositionMatrix-class`] object.
#' @param min A length-one [`numeric`] vector specifying the lower bound for
#'  element identification.
#' @param max A length-one [`numeric`] vector specifying the upper bound for
#'  element identification.
#' @param ... Currently not used.
#' @details
#'  There is no definite classification of what are the major, minor and trace
#'  elements are. By default, the following rule of thumb is used:
#'  \describe{
#'   \item{major elements}{The major elements are those that define the material
#'   under study. Major elements usually have concentrations of above 1%.}
#'   \item{minor elements}{Minor elements usually have concentrations between
#'   1% and 0.1%}
#'   \item{trace elements}{Trace elements usually have concentrations of less
#'   than 0.1%.}
#'  }
#' @note
#'  `is_oxide()` uses a regular expression (it does not check if elements exist
#'  or if stoichiometry is valid).
#' @return A [`logical`] vector.
#' @example inst/examples/ex-chemistry.R
#' @author N. Frerebeau
#' @docType methods
#' @family tools
#' @name chemistry
#' @rdname chemistry
NULL

#' @rdname chemistry
#' @aliases is_element_major-method
setGeneric(
  name = "is_element_major",
  def = function(object, ...) standardGeneric("is_element_major"),
  valueClass = "logical"
)

#' @rdname chemistry
#' @aliases is_element_minor-method
setGeneric(
  name = "is_element_minor",
  def = function(object, ...) standardGeneric("is_element_minor"),
  valueClass = "logical"
)

#' @rdname chemistry
#' @aliases is_element_trace-method
setGeneric(
  name = "is_element_trace",
  def = function(object, ...) standardGeneric("is_element_trace"),
  valueClass = "logical"
)

#' @rdname chemistry
#' @aliases is_oxide-method
setGeneric(
  name = "is_oxide",
  def = function(object, ...) standardGeneric("is_oxide"),
  valueClass = "logical"
)

# Log-Ratio ====================================================================
## LR --------------------------------------------------------------------------
#' Pairwise Log-Ratios (LR)
#'
#' Computes all pairwise log-ratio transformation.
#' @param object A [`CompositionMatrix-class`] object.
#' @param weights A [`logical`] scalar: should varying weights (column means)
#'  be computed? If `FALSE` (the default), equally-weighted parts are used.
#'  Alternatively, a positive [`numeric`] vector of weights can be specified
#'  (will be rescaled to sum to \eqn{1}). Weights will be used internally by
#'  other methods (e.g. [variance()]).
#' @param ... Currently not used.
#' @return
#'  A [`LR-class`] object.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#'
#'  Greenacre, M. J. (2021). Compositional Data Analysis. *Annual Review of
#'  Statistics and Its Application*, 8(1): 271-299.
#'  \doi{10.1146/annurev-statistics-042720-124436}.
#' @example inst/examples/ex-transform-lr.R
#' @author N. Frerebeau
#' @docType methods
#' @family log-ratio transformations
#' @aliases transform_lr-method
setGeneric(
  name = "transform_lr",
  def = function(object, ...) standardGeneric("transform_lr"),
  valueClass = "LR"
)

## CLR -------------------------------------------------------------------------
#' Centered Log-Ratios (CLR)
#'
#' Computes CLR transformation.
#' @param object A [`CompositionMatrix-class`] object.
#' @param weights A [`logical`] scalar: should varying weights (column means)
#'  be used? If `FALSE` (the default), equally-weighted parts are used.
#'  Alternatively, a positive [`numeric`] vector of weights can be specified
#'  (will be rescaled to sum to \eqn{1}).
#' @param ... Currently not used.
#' @details
#'  The CLR transformation computes the log of each part relative to the
#'  geometric mean of all parts.
#' @return
#'  A [`CLR-class`] object.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#'
#'  Greenacre, M. J. (2021). Compositional Data Analysis. *Annual Review of
#'  Statistics and Its Application*, 8(1): 271-299.
#'  \doi{10.1146/annurev-statistics-042720-124436}.
#' @example inst/examples/ex-transform-clr.R
#' @author N. Frerebeau
#' @docType methods
#' @family log-ratio transformations
#' @aliases transform_clr-method
setGeneric(
  name = "transform_clr",
  def = function(object, ...) standardGeneric("transform_clr"),
  valueClass = "CLR"
)

## ALR -------------------------------------------------------------------------
#' Additive Log-Ratios (ALR)
#'
#' Computes ALR transformation.
#' @param object A [`CompositionMatrix-class`] object.
#' @param j An [`integer`] giving the index of the rationing part (denominator).
#' @param weights A [`logical`] scalar: should varying weights (column means)
#'  be computed? If `FALSE` (the default), equally-weighted parts are used.
#'  Alternatively, a positive [`numeric`] vector of weights can be specified
#'  (will be rescaled to sum to \eqn{1}). Weights will be used internally by
#'  other methods (e.g. [variance()]).
#' @param ... Currently not used.
#' @details
#'  The ALR transformation is the logratio of a pair of parts with respect to a
#'  fixed part.
#' @return
#'  An [`ALR-class`] object.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#'
#'  Greenacre, M. J. (2021). Compositional Data Analysis. *Annual Review of
#'  Statistics and Its Application*, 8(1): 271-299.
#'  \doi{10.1146/annurev-statistics-042720-124436}.
#' @example inst/examples/ex-transform-alr.R
#' @author N. Frerebeau
#' @docType methods
#' @family log-ratio transformations
#' @aliases transform_alr-method
setGeneric(
  name = "transform_alr",
  def = function(object, ...) standardGeneric("transform_alr"),
  valueClass = "ALR"
)

## ILR -------------------------------------------------------------------------
#' Isometric Log-Ratios (ILR)
#'
#' Computes ILR transformations.
#' @param object A [`CompositionMatrix-class`] object.
#' @param ... Currently not used.
#' @details
#'  The ILR transformation provides the coordinates of any composition with
#'  respect to a given orthonormal basis. `transform_ilr()` uses the orthonormal
#'  basis (Helmert matrix) originally defined by Egozcue *et al.* (2003).
#' @return
#'  An [`ILR-class`] object.
#' @references
#'  Egozcue, J. J., Pawlowsky-Glahn, V., Mateu-Figueras, G. & Barceló-Vidal, C.
#'  (2003). Isometric Logratio Transformations for Compositional Data Analysis.
#'  *Mathematical Geology*, 35(3), 279-300. \doi{10.1023/A:1023818214614}.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#'
#'  Greenacre, M. J. (2021). Compositional Data Analysis. *Annual Review of
#'  Statistics and Its Application*, 8(1): 271-299.
#'  \doi{10.1146/annurev-statistics-042720-124436}.
#' @example inst/examples/ex-transform-ilr.R
#' @author N. Frerebeau
#' @docType methods
#' @family log-ratio transformations
#' @aliases transform_ilr-method
setGeneric(
  name = "transform_ilr",
  def = function(object, ...) standardGeneric("transform_ilr"),
  valueClass = "ILR"
)

#' Univariate Isometric Log-Ratios (ILR)
#'
#' Computes univariate ILR coordinates.
#' @param object A [`CompositionMatrix-class`] object.
#' @param ... Currently not used.
#' @details
#'  The ILR transformation provides the coordinates of any composition with
#'  respect to a given orthonormal basis. `transform_ilr()` uses the orthonormal
#'  basis (Helmert matrix) originally defined by Egozcue *et al.* (2003).
#' @return
#'  Same as `object`.
#' @references
#'  Filzmoser, P., Hron, K. & Reimann, C. (2009). Univariate Statistical
#'  Analysis of Environmental (Compositional) Data: Problems and Possibilities.
#'  *Science of The Total Environment*, 407(23), 6100-6108.
#'  \doi{10.1016/j.scitotenv.2009.08.008}.
#'
#'  Filzmoser, P., Hron, K. & Reimann, C. (2010). The Bivariate Statistical
#'  Analysis of Environmental (Compositional) Data. *Science of The Total
#'  Environment*, 408(19), 4230-4238. \doi{10.1016/j.scitotenv.2010.05.011}.
#' @example inst/examples/ex-transform-ilr.R
#' @author N. Frerebeau
#' @docType methods
# @family log-ratio transformations
#' @aliases univariate_ilr-method
#' @keywords internal
setGeneric(
  name = "univariate_ilr",
  def = function(object, ...) standardGeneric("univariate_ilr")
)

## PLR -------------------------------------------------------------------------
#' Pivot Log-Ratios (PLR)
#'
#' Computes PLR transformations.
#' @param object A [`CompositionMatrix-class`] object.
#' @param pivot An [`integer`] giving the index of the pivotal variable.
#' @param ... Currently not used.
#' @return
#'  A [`PLR-class`] object.
#' @references
#'  Fišerová, E. & Hron, K. (2011). On the Interpretation of Orthonormal
#'  Coordinates for Compositional Data. *Mathematical Geosciences*, 43(4),
#'  455‑468. \doi{10.1007/s11004-011-9333-x}.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#'
#'  Greenacre, M. J. (2021). Compositional Data Analysis. *Annual Review of
#'  Statistics and Its Application*, 8(1): 271-299.
#'  \doi{10.1146/annurev-statistics-042720-124436}.
#'
#'  Hron, K., Filzmoser, P., de Caritat, P., Fišerová, E. & Gardlo, A. (2017).
#'  Weighted Pivot Coordinates for Compositional Data and Their Application to
#'  Geochemical Mapping. *Mathematical Geosciences*, 49(6), 797-814.
#'  \doi{10.1007/s11004-017-9684-z}.
#' @example inst/examples/ex-transform-ilr.R
#' @author N. Frerebeau
#' @docType methods
#' @family log-ratio transformations
#' @aliases transform_plr-method
setGeneric(
  name = "transform_plr",
  def = function(object, ...) standardGeneric("transform_plr"),
  valueClass = "PLR"
)

## Inverse ---------------------------------------------------------------------
#' Inverse Log-Ratio Transformation
#'
#' Computes inverse log-ratio transformations.
#' @param object A [`LogRatio-class`] object.
#' @param origin A [`LogRatio-class`] object to be used for the inverse
#'  transformation.
#' @param ... Currently not used.
#' @return
#'  A [`CompositionMatrix-class`] object.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall.
#'
#'  Egozcue, J. J., Pawlowsky-Glahn, V., Mateu-Figueras, G. & Barceló-Vidal, C.
#'  (2003). Isometric Logratio Transformations for Compositional Data Analysis.
#'  *Mathematical Geology*, 35(3), 279-300. \doi{10.1023/A:1023818214614}.
#'
#'  Fišerová, E. & Hron, K. (2011). On the Interpretation of Orthonormal
#'  Coordinates for Compositional Data. *Mathematical Geosciences*, 43(4),
#'  455‑468. \doi{10.1007/s11004-011-9333-x}.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#' @example inst/examples/ex-transform-clr.R
#' @author N. Frerebeau
#' @docType methods
#' @family log-ratio transformations
#' @aliases transform_inverse-method
setGeneric(
  name = "transform_inverse",
  def = function(object, origin, ...) standardGeneric("transform_inverse"),
  valueClass = "matrix"
)

# Statistics ===================================================================
#' Compute Summary Statistics of Data Subsets
#'
#' Splits the data into subsets, computes summary statistics for each, and
#' returns the result.
#' @param x A [`CompositionMatrix-class`] object.
#' @param by A `vector` or a list of grouping elements, each as long as the
#'  variables in `x` (see [group()]).
#' @param FUN A [`function`] to compute the summary statistics.
#' @param simplify A [`logical`] scalar: should the results be simplified to a
#'  matrix if possible?
#' @param drop A [`logical`] scalar indicating whether to drop unused
#'  combinations of grouping values.
#' @param ... Further arguments to be passed to `FUN`.
#' @return A [`matrix`].
#' @example inst/examples/ex-aggregate.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name aggregate
#' @rdname aggregate
NULL

#' Compositional Mean
#'
#' @param x A [`CompositionMatrix-class`] object.
#' @inheritParams gmean
#' @param ... Currently not used.
#' @details
#'  Closed vector of the columns geometric means.
#' @return A [`numeric`] vector.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall, p. 64-91.
#' @example inst/examples/ex-mean.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name mean
#' @rdname mean
NULL

#' Sample Quantiles
#'
#' @param x A [`CompositionMatrix-class`] object.
#' @param probs A [`numeric`] vector of probabilities with values in \eqn{[0,1]}.
#' @param na.rm A [`logical`] scalar: should missing values be removed?
#' @param names A [`logical`] scalar: should results be named?
#' @param ... Currently not used.
#' @return A [`numeric`] matrix.
#' @references
#'  Filzmoser, P., Hron, K. & Reimann, C. (2009). Univariate Statistical
#'  Analysis of Environmental (Compositional) Data: Problems and Possibilities.
#'  *Science of The Total Environment*, 407(23): 6100-6108.
#'  \doi{10.1016/j.scitotenv.2009.08.008}.
#' @example inst/examples/ex-mean.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name quantile
#' @rdname quantile
NULL

#' Compositional Mean of Data Subsets
#'
#' Splits the data into subsets and computes compositional mean for each.
#' @param x A [`CompositionMatrix-class`] object.
#' @param by A `vector` or a list of grouping elements, each as long as the
#'  variables in `x` (see [group()]).
#' @inheritParams gmean
#' @param verbose A [`logical`] scalar: should \R report extra information
#'  on progress?
#' @param ... Currently not used.
#' @return A [`CompositionMatrix-class`] object.
#' @seealso [mean()], [aggregate()]
#' @example inst/examples/ex-condense.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases condense-method
setGeneric(
  name = "condense",
  def = function(x, ...) standardGeneric("condense")
  # valueClass = "CompositionMatrix"
)

#' Marginal Compositions
#'
#' @param x A [`CompositionMatrix-class`] object.
#' @param parts An [`integer`] or a [`character`] vector specifying the columns
#'  to be selected.
#' @param name A [`character`] string giving the name of the amalgamation
#'  column.
#' @param ... Currently not used.
#' @return A [`CompositionMatrix-class`] object.
#' @example inst/examples/ex-margin.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases margin-method
setGeneric(
  name = "margin",
  def = function(x, ...) standardGeneric("margin"),
  valueClass = "CompositionMatrix"
)

#' Covariance Matrix
#'
#' Computes the (centered) log-ratio covariance matrix (see below).
#' @param x A [`CompositionMatrix-class`] object.
#' @param center A [`logical`] scalar: should the *centered* log-ratio
#'  covariance matrix be computed?
#' @param method A [`character`] string indicating which covariance is to be
#'  computed (see [stats::cov()]).
#' @param ... Currently not used.
#' @return A [`matrix`].
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall, p. 64-91.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#' @example inst/examples/ex-covariance.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases covariance-method
setGeneric(
  name = "covariance",
  def = function(x, ...) standardGeneric("covariance"),
  valueClass = "matrix"
)

#' Log-Ratios Variances
#'
#' Computes log-ratio (weighted) variances.
#' @param x A [`CompositionMatrix-class`] object.
#' @param row_weights A [`numeric`] vector of row weights. If `NULL` (the
#'  default), equal weights are used.
#' @param column_weights A [`logical`] scalar: should the weights of the
#'  log-ratio be used? If `FALSE`, equally-weighted parts are used.
#'  Alternatively, a positive [`numeric`] vector of weights can be specified.
#' @param ... Currently not used.
#' @return A [`numeric`] vector of individual variances.
#' @references
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#' @example inst/examples/ex-variance.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases variance-method
setGeneric(
  name = "variance",
  def = function(x, ...) standardGeneric("variance"),
  valueClass = "numeric"
)

#' Total Variance
#'
#' Computes the total (or metric) variance, a global measure of spread.
#' @inheritParams variance
#' @param sd A [`logical`] scalar: should the metric standard deviation be
#'  returned instead of the metric variance?
#' @return A [`numeric`] vector.
#' @details
#'  Two methods are available, see below.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall, p. 64-91.
#'
#'  Aitchison, J. (1997). The One-Hour Course in Compositional Data Analysis or
#'  Compositional Data Analysis Is Simple. In V. Pawlowsky-Glahn (ed.),
#'  *IAMG'97*. Barcelona: International Center for Numerical Methods in
#'  Engineering (CIMNE), p. 3-35.
#'
#'  Boogaart, K. G. van den & Tolosana-Delgado, R. (2013). *Analyzing
#'  Compositional Data with R*. Berlin Heidelberg: Springer-Verlag.
#'  \doi{10.1007/978-3-642-36809-7}.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#'
#'  Hron, K. & Kubáček. L. (2011). Statistical Properties of the Total Variation
#'  Estimator for Compositional Data. *Metrika*, 74 (2): 221-230.
#'  \doi{10.1007/s00184-010-0299-3}.
#'
#'  Pawlowsky-Glahn, V. & Egozcue, J. J. (2001). Geometric Approach to
#'  Statistical Analysis on the Simplex. *Stochastic Environmental Research and
#'  Risk Assessment*, 15(5): 384-398. \doi{10.1007/s004770100077}.
#' @example inst/examples/ex-variance.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases variance_total-method
setGeneric(
  name = "variance_total",
  def = function(x, ...) standardGeneric("variance_total"),
  valueClass = "numeric"
)

#' Variation Matrix
#'
#' Computes the variation matrix (Aitchison 1986, definition 4.4).
#' @param x A [`CompositionMatrix-class`] object.
#' @param ... Currently not used.
#' @return A [`matrix`].
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall, p. 64-91.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#' @example inst/examples/ex-variation.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases variation-method
setGeneric(
  name = "variation",
  def = function(x, ...) standardGeneric("variation"),
  valueClass = "matrix"
)

#' Proportionality Index of Parts (PIP)
#'
#' Computes an index of association between parts.
#' @param x A [`CompositionMatrix-class`] object.
#' @param ... Currently not used.
#' @return A [`matrix`].
#' @details
#'  The proportionality index of parts (PIP) is based on the
#'  [variation matrix][variation()], but maintains the range of values whithin
#'  \eqn{(0,1)}.
#' @references
#'  Egozcue, J. J.. & Pawlowsky-Glahn, V. (2023). Subcompositional Coherence
#'  and and a Novel Proportionality Index of Parts. *SORT*, 47(2): 229-244.
#'  \doi{10.57645/20.8080.02.7}.
#' @example inst/examples/ex-variation.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases pip-method
setGeneric(
  name = "pip",
  def = function(x, ...) standardGeneric("pip"),
  valueClass = "matrix"
)

# Variation Array
#
# Computes the compositional variation array.
# @param object A [`CompositionMatrix-class`] object.
# @param ... Currently not used.
# @return A [`matrix`].
# @details
#  The compositional variation array is a square matrix where the upper
#  triangular part displays the pairwise log-ratio variances and the lower
#  triangular part displays the pairwise log-ratio means.
# @references
#  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#  London: Chapman and Hall, p. 64-91.
# @example inst/examples/ex-variation_array.R
# @author N. Frerebeau
# @docType methods
# @family statistics
# @aliases variation_array-method
# setGeneric(
#   name = "variation_array",
#   def = function(object, ...) standardGeneric("variation_array"),
#   valueClass = "matrix"
# )

#' Scaling and Centering of Compositional Data
#'
#' @param x A [`CompositionMatrix-class`] object.
#' @param center A [`logical`] scalar or a [`numeric`] vector giving the center
#'  to be substracted.
#' @param scale A [`logical`] scalar or a length-one [`numeric`] vector giving a
#'  scaling factor for multiplication.
#' @return A [`CompositionMatrix-class`] object.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall, p. 64-91.
#'
#'  Boogaart, K. G. van den & Tolosana-Delgado, R. (2013). *Analyzing
#'  Compositional Data with R*. Berlin Heidelberg: Springer-Verlag.
#'  \doi{10.1007/978-3-642-36809-7}.
#' @example inst/examples/ex-scale.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name scale
#' @rdname scale
NULL

# Distances ====================================================================
#' Distances
#'
#' Computes the distances between all rows of in `x`.
#' @param x A [`CompositionMatrix-class`] object.
#' @param method A [`character`] string specifying the distance measure to be
#'  used. See [stats::dist()] for the available distances.
#' @param diag A [`logical`] scalar indicating whether the diagonal of the
#'  distance matrix should be printed.
#' @param upper A [`logical`] scalar indicating whether the upper triangle of
#'  the distance matrix should be printed.
#' @param p An [`integer`] giving the power of the Minkowski distance.
#' @details
#'  Distances are computed on [CLR-transformed][transform_clr] data.
#' @return A [`stats::dist`] object.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall, p. 64-91.
#'
#'  Greenacre, M. J. (2019). *Compositional Data Analysis in Practice*.
#'  Boca Raton: CRC Press.
#' @example inst/examples/ex-dist.R
#' @seealso [stats::dist()]
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name dist
#' @rdname dist
NULL

#' Mahalanobis Distance
#'
#' Computes the squared Mahalanobis distance of all rows in `x`.
#' @param x A [`CompositionMatrix-class`] or an [`ILR-class`] object.
#' @param center A [`numeric`] vector giving the mean vector of the
#'  distribution. If missing, will be estimated from `x`.
#' @param cov A [`numeric`] matrix giving the covariance of the
#'  distribution. If missing, will be estimated from `x`.
#' @param robust A [`logical`] scalar: should robust location and scatter
#'  estimation be used?
#' @param method A [`character`] string specifying the method to be used.
#'  It must be one of "`mve`" (minimum volume ellipsoid) or "`mcd`" (minimum
#'  covariance determinant). Only used if `robust` is `TRUE`.
#' @param ... Extra parameters to be passed to [MASS::cov.rob()].
#'  Only used if `robust` is `TRUE`.
#' @return A [`numeric`] vector.
#' @example inst/examples/ex-mahalanobis.R
#' @seealso [stats::mahalanobis()]
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @name mahalanobis
#' @rdname mahalanobis
NULL

# Plot =========================================================================
## Barplot ---------------------------------------------------------------------
#' Barplot of Compositional Data
#'
#' Displays a compositional bar chart.
#' @param height A [`CompositionMatrix-class`] object.
#' @param order_columns A [`logical`] scalar: should should columns be reorderd?
#' @param order_rows An [`integer`] vector giving the index of the column to be
#'  used for the ordering of the data.
#' @param decreasing A [`logical`] scalar: should the sort order of rows be
#'  increasing or decreasing?
#' @param names A [`logical`] scalar: should row names be displayed?
#' @param space A length-one [`numeric`] vector giving the the amount of space
#'  (as a fraction of the width of a bar) left between each bar
#'  (defaults to \eqn{0.2}).
#' @param offset A length-one [`numeric`] vector giving the the amount of space
#'  (as a fraction) left between groups (defaults to \eqn{0.025}). Only used if
#'  `groups` is not `NULL`.
#' @param palette_color A palette [`function`] that when called with a single
#'  argument returns a `character` vector of colors.
#' @param border The color to draw the borders.
#' @param axes A [`logical`] scalar: should axes be drawn on the plot?
#' @param legend A [`logical`] scalar: should the legend be displayed?
#' @param ... Further graphical parameters.
#' @return
#'  `barplot()` is called for its side-effects: is results in a graphic being
#'  displayed (invisibly return `height`).
#' @example inst/examples/ex-barplot.R
#' @author N. Frerebeau
#' @docType methods
#' @family plot methods
#' @name barplot
#' @rdname barplot
NULL

## Histogram -------------------------------------------------------------------
#' Histogram of Compositional Data
#'
#' Produces an histogram of univariate ILR data (see Filzmoser *et al.*, 2009).
#' @param x A [`CompositionMatrix-class`] object.
#' @param select A length-one `vector` of column indices.
#' @param breaks An object specifying how to compute the breakpoints
#'  (see [graphics::hist()]).
#' @param freq A [`logical`] scalar: should absolute frequencies (counts) be
#'  displayed? If `FALSE` (the default), relative frequencies (probabilities)
#'  are displayed (see [graphics::hist()]).
#' @param labels A [`logical`] scalar: should labels be drawn on top of bars?
#'  If `TRUE`, draw the counts or rounded densities; if `labels` is a
#'  `character` vector, draw itself.
#' @param main A [`character`] string giving a main title for the plot.
#' @param sub A [`character`] string giving a subtitle for the plot.
#' @param ann A [`logical`] scalar: should the default annotation (title and x
#'  and y 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 ... Further graphical parameters.
#' @return
#'  `hist()` is called for its side-effects: is results in a graphic being
#'  displayed (invisibly return `x`).
#' @references
#'  Filzmoser, P., Hron, K. & Reimann, C. (2009). Univariate Statistical
#'  Analysis of Environmental (Compositional) Data: Problems and Possibilities.
#'  *Science of The Total Environment*, 407(23): 6100-6108.
#'  \doi{10.1016/j.scitotenv.2009.08.008}.
#' @example inst/examples/ex-hist.R
#' @author N. Frerebeau
#' @docType methods
#' @family plot methods
#' @name hist
#' @rdname hist
NULL

## Ternary ---------------------------------------------------------------------
#' Plot Compositional Data
#'
#' Displays a matrix of ternary plots.
#' @param x A [`CompositionMatrix-class`] object.
#' @param palette_color A palette [`function`] that when called with a single
#'  argument returns a `character` vector of colors.
#' @param palette_symbol A palette [`function`] that when called with a single
#'  argument returns a vector of symbols.
#' @inheritParams isopleuros::ternary_pairs
#' @return
#'  `plot()` is called for its side-effects: is results in a graphic being
#'  displayed (invisibly return `x`).
#' @seealso [isopleuros::ternary_pairs()], [isopleuros::ternary_plot()]
#' @example inst/examples/ex-pairs.R
#' @author N. Frerebeau
#' @docType methods
#' @family plot methods
#' @name pairs
#' @rdname pairs
NULL

## Density ---------------------------------------------------------------------
#' Plot Log-Ratios
#'
#' Displays a scatter plot.
#' @param x A [`LogRatio-class`] object.
#' @param factor,amount A length-one [`numeric`] vector specifying the amount of
#'  jitter (see [jitter()]).
#' @param palette_color A palette [`function`] that when called with a single
#'  argument returns a `character` vector of colors (only used if `x`
#'  [is grouped][group()]).
#' @param palette_symbol A palette [`function`] that when called with a single
#'  argument returns a vector of symbols (only used if `x`
#'  [is grouped][group()]).
#' @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
#'  and y 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 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.
#' @return
#'  `plot()` is called for its side-effects: is results in a graphic being
#'  displayed (invisibly return `x`).
#' @example inst/examples/ex-plot.R
#' @author N. Frerebeau
#' @docType methods
#' @family plot methods
#' @name plot
#' @rdname plot
NULL

## Graph -----------------------------------------------------------------------
#' Graph of Log-ratios
#'
#' Produces a graph of log-ratios.
#' @param object A [`LogRatio-class`] object.
#' @param ... Currently not used.
#' @return
#'  An \pkg{igraph} graph object.
#' @example inst/examples/ex-graph.R
#' @author N. Frerebeau
#' @docType methods
#' @family plot methods
#' @aliases as_graph-method
setGeneric(
  name = "as_graph",
  def = function(object, ...) standardGeneric("as_graph")
)

# Multivariate analysis ========================================================
#' Principal Components Analysis
#'
#' Computes a principal components analysis based on the singular value
#' decomposition.
#' @param object A [`CompositionMatrix-class`] or [`LogRatio-class`] object.
#' @inheritParams dimensio::pca
#' @return
#'  A [`dimensio::PCA-class`] object. See [dimensio::pca()] for details.
#' @references
#'  Aitchison, J. and Greenacre, M. (2002). Biplots of compositional data.
#'  *Journal of the Royal Statistical Society: Series C (Applied Statistics)*,
#'  51: 375-392. \doi{10.1111/1467-9876.00275}.
#'
#'  Filzmoser, P., Hron, K. and Reimann, C. (2009). Principal component analysis
#'  for compositional data with outliers. *Environmetrics*, 20: 621-632.
#'  \doi{10.1002/env.966}.
#' @example inst/examples/ex-pca.R
#' @seealso [dimensio::pca()], [dimensio::biplot()], [dimensio::screeplot()],
#'  [dimensio::viz_individuals()], [dimensio::viz_variables()]
#' @author N. Frerebeau
#' @docType methods
#' @family multivariate analysis
#' @name pca
#' @rdname pca
#' @aliases lra
NULL

# Missign Values ===============================================================
#' Missing Values Policy
#'
#' @details
#'  Compositional data are quantitative positive descriptions of the parts
#'  of some whole, carrying relative, rather than absolute, information
#'  (ie. only relative changes are relevant; Aitchison 1986).
#'
#'  Basically, three situations can be outlined regarding missing values in
#'  compositions:
#'
#'  * Unobserved quantities.
#'  * Amounts observed, but which happen to be below the detection limit
#'    (thus interpreted as small unknown values).
#'  * Absolutely zero quantities.
#'
#' These situations can be represented in several ways:
#'
#'  * The presence of zeros.
#'  * The presence of missing values (`NA`).
#'
#'  When creating a [`CompositionMatrix-class`] object, the presence of zero
#'  and [`NA`] values is allowed: this makes it possible to explore and
#'  visualize the data while preserving the missing structure. However, **the
#'  user must deal with these missing values before proceeding further** (e.g.
#'  by removing incomplete cases or replacing the values concerned): log-ratio
#'  transformations cannot be computed in the presence of zeros or missing
#'  values.
#' @note
#'  If you need more advanced features (e.g. imputation of missing values),
#'  you should consider the \pkg{compositions} or \pkg{robCompositions} package.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall.
#' @family imputation methods
#' @name missing
#' @rdname missing
NULL

#' Zero-Replacement
#'
#' Multiplicative replacement of zeros.
#' @param x A [`CompositionMatrix-class`] object.
#' @param value A [`numeric`] vector giving the detection limits of each part
#'  (in \eqn{(0,1)}).
#' @param delta A [`numeric`] vector specifying the fraction of the detection
#'  limit to be used in replacement.
#' @return
#'  An [`CompositionMatrix-class`] object, where all zero values have been
#'  replaced.
#' @references
#'  Aitchison, J. (1986). *The Statistical Analysis of Compositional Data*.
#'  London: Chapman and Hall.
#'
#'  Martín-Fernández, J. A., Barceló-Vidal, C. & Pawlowsky-Glahn, V. (2003).
#'  Dealing with Zeros and Missing Values in Compositional Data Sets Using
#'  Nonparametric Imputation. *Mathematical Geology*, 35(3): 253-278.
#'  \doi{10.1023/A:1023866030544}.
#' @example inst/examples/ex-zero.R
#' @author N. Frerebeau
#' @docType methods
#' @family imputation methods
#' @name replace_zero
#' @rdname replace_zero
NULL

#' Missing Values Replacement
#'
#' Multiplicative replacement of missing values.
#' @param x A [`CompositionMatrix-class`] object.
#' @param value A [`numeric`] vector giving the replacement values.
#' @return
#'  An [`CompositionMatrix-class`] object, where all missing values have been
#'  replaced.
#' @references
#'  Martín-Fernández, J. A., Barceló-Vidal, C. & Pawlowsky-Glahn, V. (2003).
#'  Dealing with Zeros and Missing Values in Compositional Data Sets Using
#'  Nonparametric Imputation. *Mathematical Geology*, 35(3): 253-278.
#'  \doi{10.1023/A:1023866030544}.
#' @example inst/examples/ex-missing.R
#' @author N. Frerebeau
#' @docType methods
#' @family imputation methods
#' @name replace_NA
#' @rdname replace_NA
NULL

# Outliers =====================================================================
#' Outlier Detection
#'
#' @param object A [`CompositionMatrix-class`].
#' @param reference A [`CompositionMatrix-class`]. If missing, `object` is used.
#' @param robust A [`logical`] scalar: should robust estimators be used?
#' @param method A [`character`] string specifying the method to be used.
#'  It must be one of "`mve`" (minimum volume ellipsoid) or "`mcd`" (minimum
#'  covariance determinant; see [MASS::cov.rob()]).
#'  Only used if `robust` is `TRUE`.
#' @param quantile A length-one [`numeric`] vector giving the significance level.
#'  `quantile` is used as a cut-off value for outlier detection: observations
#'  with larger (squared) Mahalanobis distance are considered as potential
#'  outliers.
#' @param ... Further parameters to be passed to [MASS::cov.rob()].
#' @details
#'  An outlier can be defined as having a very large Mahalanobis distance from
#'  all observations. In this way, a certain proportion of the observations can
#'  be identified, e.g. the top 2% of values (i.e. values above the 0.98th
#'  percentile of the Chi-2 distribution).
#'
#'  On the one hand, the Mahalanobis distance is likely to be strongly
#'  affected by the presence of outliers. Rousseeuw and van Zomeren (1990) thus
#'  recommend using robust methods (which are not excessively affected by the
#'  presence of outliers).
#'
#'  On the other hand, the choice of the threshold for classifying an
#'  observation as an outlier should be discussed. There is no apparent reason
#'  why a particular threshold should be applicable to all data sets
#'  (Filzmoser, Garrett, and Reimann 2005).
#' @return
#'  * `detect_outlier()` returns an [`OutlierIndex-class`] object.
#'  * `is_outlier()` returns a [`logical`] vector.
#' @references
#'  Filzmoser, P., Garrett, R. G. & Reimann, C. (2005). Multivariate outlier
#'  detection in exploration geochemistry. *Computers & Geosciences*,
#'  31(5), 579-587. \doi{10.1016/j.cageo.2004.11.013}.
#'
#'  Filzmoser, P. & Hron, K. (2008). Outlier Detection for Compositional Data
#'  Using Robust Methods. *Mathematical Geosciences*, 40(3), 233-248.
#'  \doi{10.1007/s11004-007-9141-5}.
#'
#'  Filzmoser, P., Hron, K. & Reimann, C. (2012). Interpretation of multivariate
#'  outliers for compositional data. *Computers & Geosciences*, 39, 77-85.
#'  \doi{10.1016/j.cageo.2011.06.014}.
#'
#'  Rousseeuw, P. J. & van Zomeren, B. C. (1990). Unmasking Multivariate Outliers
#'  and Leverage Points. *Journal of the American Statistical Association*,
#'  85(411): 633-639. \doi{10.1080/01621459.1990.10474920}.
#'
#'  Santos, F. (2020). Modern methods for old data: An overview of some robust
#'  methods for outliers detection with applications in osteology. *Journal of
#'  Archaeological Science: Reports*, 32, 102423.
#'  \doi{10.1016/j.jasrep.2020.102423}.
#' @example inst/examples/ex-outliers.R
#' @author N. Frerebeau
#' @docType methods
#' @family outlier detection methods
#' @aliases detect_outlier-method
setGeneric(
  name = "detect_outlier",
  def = function(object, reference, ...) standardGeneric("detect_outlier")
)

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

#' Plot Outliers
#'
#' @param x An [`OutlierIndex-class`] object.
#' @param type A [`character`] string specifying the type of plot that should be
#'  made. It must be one of "`dotchart`" or "`distance`".
#'  Any unambiguous substring can be given.
#' @param robust A [`logical`] scalar: should robust Mahalanobis distances be
#'  displayed? Only used if `type` is "`dotchart`".
#' @param symbols A lenth-three vector of symbol specification for non-outliers
#'  and outliers (resp.).
#' @param xlim A length-two [`numeric`] vector giving the x limits of the plot.
#'  The default value, `NULL`, indicates that the range of the
#'  [finite][is.finite()] values to be plotted should be used.
#' @param ylim A length-two [`numeric`] vector giving the y limits of the plot.
#'  The default value, `NULL`, indicates that the range of the
#'  [finite][is.finite()] values to be plotted should be used.
#' @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
#'  and y axis labels) appear on the plot?
#' @param axes A [`logical`] scalar: should axes be drawn on 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 frame.plot A [`logical`] scalar: should a box be drawn around the
#'  plot?
#' @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 parameters to be passed to [graphics::points()].
#' @return
#'  `plot()` is called for its side-effects: is results in a graphic being
#'  displayed (invisibly return `x`).
#' @references
#'  Filzmoser, P., Garrett, R. G. & Reimann, C. (2005). Multivariate outlier
#'  detection in exploration geochemistry. *Computers & Geosciences*,
#'  31(5), 579-587. \doi{10.1016/j.cageo.2004.11.013}.
#'
#'  Filzmoser, P. & Hron, K. (2008). Outlier Detection for Compositional Data
#'  Using Robust Methods. *Mathematical Geosciences*, 40(3), 233-248.
#'  \doi{10.1007/s11004-007-9141-5}.
#'
#'  Filzmoser, P., Hron, K. & Reimann, C. (2012). Interpretation of multivariate
#'  outliers for compositional data. *Computers & Geosciences*, 39, 77-85.
#'  \doi{10.1016/j.cageo.2011.06.014}.
#' @example inst/examples/ex-outliers.R
#' @author N. Frerebeau
#' @docType methods
#' @family outlier detection methods
#' @name plot_outlier
#' @rdname plot_outlier
NULL

# Sourcing =====================================================================
#' Mixed-Mode Analysis
#'
#' Mixes chemical and petrographic matrices.
#' @param x A [`matrix`] of chemical compositional data or a
#'  [dissimilarity matrix][stats::dist] for these chemical compositional data.
#' @param y A [`matrix`] of coded mineralogical binary data or a
#'  [dissimilarity matrix][stats::dist] for these mineralogical data.
#' @param lambda A length-one [`numeric`] vector giving a weighting factor.
#' @param mu A length-one [`numeric`] vector that lies between 0 and 1 giving
#'  the mixing parameter.
#' @param ... Extra parameters to be passed to [cluster::daisy()].
#' @return
#'  A [stats::dist] object.
#' @references
#'  Baxter, M. J., Beardah, C. C., Papageorgiou, I., Cau, M. A., Day, P. M. &
#'  Kilikoglou, V. (2008). On Statistical Approaches to the Study of Ceramic
#'  Artefacts Using Geochemical and Petrographic Data. *Archaeometry*, 50(1):
#'  142-157. \doi{10.1111/j.1475-4754.2007.00359.x}.
#'
#'  Beardah, C. C., Baxter, M. J., Papageorgiou, I. & Cau, M. A. (2003).
#'  "Mixed-Mode" Approaches to the Grouping of Ceramic Artefacts Using S-Plus.
#'  In M. Doerr and A. Sarris, *The Digital Heritage of Archaeology*, p. 261-266.
#'  Athens: Archive of Monuments and Publications, Hellenic Ministry of Culture.
#'
#'  Gower, J. C. (1971). A general coefficient of similarity and some of its
#'  properties. *Biometrics*, 27(4):857-874. \doi{10.2307/2528823}.
#' @note
#'  **Experimental.**
#' @example inst/examples/ex-mix.R
#' @author N. Frerebeau
#' @docType methods
#' @family sourcing methods
#' @aliases mix-method
setGeneric(
  name = "mix",
  def = function(x, y, ...) standardGeneric("mix"),
  valueClass = "dist"
)

# PLOT
#' @include AllGenerics.R
NULL

# LogRatio =====================================================================
#' @export
#' @method plot LogRatio
plot.LogRatio <- function(x, ..., factor = 1, amount = NULL,
                          palette_color = palette_color_discrete(),
                          palette_symbol = palette_shape(),
                          xlab = NULL, ylab = NULL,
                          main = NULL, sub = NULL, ann = graphics::par("ann"),
                          axes = TRUE, frame.plot = axes,
                          legend = list(x = "topright")) {
  ## Get data
  xy <- data.frame(
    x = jitter(as.vector(col(x)), factor = factor, amount = amount),
    y = as.vector(x)
  )

  ## Graphical parameters
  col <- list(...)$col %||% graphics::par("col")
  bg <- list(...)$bg %||% graphics::par("bg")
  pch <- list(...)$pch %||% graphics::par("pch")
  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")

  ## Grouping
  if (is_grouped(x)) {
    lvl <- group_names(x)
    col <- palette_color(lvl)
    bg <- grDevices::adjustcolor(col, alpha.f = 0.5)
    pch <- palette_symbol(lvl)
  }

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

  ## Set plotting coordinates
  xlim <- range(xy$x)
  ylim <- range(xy$y)
  graphics::plot.window(xlim = xlim, ylim = ylim)

  ## Plot
  graphics::points(x = xy$x, y = xy$y, pch = pch, col = col, bg = bg)

  ## Construct axis
  if (axes) {
    graphics::axis(side = 1, at = seq_len(ncol(x)),
                   labels = labels(x), las = 1,
                   cex.axis = cex.axis, col.axis = col.axis,
                   font.axis = font.axis)
    graphics::axis(side = 2, las = 1)
  }

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

  ## Add annotation
  if (ann) {
    ylab <- ylab %||% get_transformation(x)
    graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
  }

  ## Add legend
  if (is.list(legend) && is_grouped(x)) {
    ## Compute legend position
    args <- list(x = "topright", legend = unique(lvl), pch = unique(pch),
                 col = unique(col), bg = unique(bg), bty = "n")
    args <- utils::modifyList(args, legend)

    ## Plot legend
    do.call(graphics::legend, args = args)
  }

  invisible(x)
}

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

# GRAPH
#' @include AllGenerics.R
NULL

#' @export
#' @rdname as_graph
#' @aliases as_graph,LR-method
setMethod(
  f = "as_graph",
  signature = c(object = "LR"),
  definition = function(object) {
    ## Validation
    arkhe::assert_package("igraph")

    ratio <- labels(object)
    edges <- do.call(rbind, strsplit(ratio, "/"))
    edges <- edges[, c(2, 1)]
    igraph::graph_from_edgelist(edges, directed = FALSE)
  }
)

#' @export
#' @rdname as_graph
#' @aliases as_graph,ALR-method
setMethod(
  f = "as_graph",
  signature = c(object = "ALR"),
  definition = function(object) {
    ## Validation
    arkhe::assert_package("igraph")

    ratio <- labels(object)
    edges <- do.call(rbind, strsplit(ratio, "/"))
    edges <- edges[, c(2, 1)]
    igraph::graph_from_edgelist(edges, directed = TRUE)
  }
)

#' @export
#' @rdname as_graph
#' @aliases as_graph,ILR-method
setMethod(
  f = "as_graph",
  signature = c(object = "ILR"),
  definition = function(object) {
    ## Validation
    arkhe::assert_package("igraph")

    ratio <- labels(object)
    ratio <- gsub(pattern = "[\\(\\)]", replacement =  "", x = ratio)
    edges <- lapply(
      X = strsplit(ratio, "/"),
      FUN = function(x) {
        a <- unlist(strsplit(x[[1]], ","))
        b <- unlist(strsplit(x[[2]], ","))
        expand.grid(b, a)
      }
    )
    edges <- do.call(rbind, edges)
    igraph::graph_from_data_frame(edges, directed = TRUE)
  }
)

# DATA TRANSFORMATION: ISOMETRIC LOG RATIO
#' @include AllGenerics.R
NULL

# ILR ==========================================================================
ilr_base <- function(D, method = "basic") {
  ## Validation
  method <- match.arg(method, several.ok = FALSE)

  seq_parts <- seq_len(D - 1)

  ## Original ILR transformation defined by Egozcue et al. 2003
  if (method == "basic") {
    ## Helmert matrix (rotation matrix)
    H <- stats::contr.helmert(D)                  # D x D-1
    H <- t(H) / sqrt((seq_parts + 1) * seq_parts) # D-1 x D

    ## Center
    M <- diag(x = 1, nrow = D) - matrix(data = 1 / D, nrow = D, ncol = D)
    V <- tcrossprod(M, H)
  }

  V
}

#' @export
#' @rdname transform_ilr
#' @aliases transform_ilr,CompositionMatrix,missing-method
setMethod(
  f = "transform_ilr",
  signature = c(object = "CompositionMatrix"),
  definition = function(object) {
    weights <- rep(1 / ncol(object), ncol(object))
    base <- ilr_base(D = ncol(object), method = "basic")
    .transform_ilr(object, base, weights)
  }
)

#' @export
#' @rdname transform_ilr
#' @aliases transform_ilr,GroupedComposition,missing-method
setMethod(
  f = "transform_ilr",
  signature = c(object = "GroupedComposition"),
  definition = function(object) {
    z <- methods::callNextMethod()
    .GroupedILR(z, group_indices = group_indices(object),
                group_levels = group_levels(object))
  }
)

#' @export
#' @rdname transform_ilr
#' @aliases transform_ilr,CLR,missing-method
setMethod(
  f = "transform_ilr",
  signature = c(object = "CLR"),
  definition = function(object) {
    weights <- object@weights
    base <- ilr_base(D = ncol(object), method = "basic")
    object@.Data <- exp(object@.Data)
    .transform_ilr(object, base, weights)
  }
)

#' @export
#' @rdname transform_ilr
#' @aliases transform_ilr,GroupedCLR,missing-method
setMethod(
  f = "transform_ilr",
  signature = c(object = "GroupedCLR"),
  definition = function(object) {
    z <- methods::callNextMethod()
    .GroupedILR(z, group_indices = group_indices(object),
                group_levels = group_levels(object))
  }
)

#' @export
#' @rdname transform_ilr
#' @aliases transform_ilr,ALR,missing-method
setMethod(
  f = "transform_ilr",
  signature = c(object = "ALR"),
  definition = function(object) {
    object <- transform_clr(object)
    methods::callGeneric(object)
  }
)

#' @export
#' @rdname transform_ilr
#' @aliases transform_ilr,GroupedALR,missing-method
setMethod(
  f = "transform_ilr",
  signature = c(object = "GroupedALR"),
  definition = function(object) {
    z <- methods::callNextMethod()
    .GroupedILR(z, group_indices = group_indices(object),
                group_levels = group_levels(object))
  }
)

.transform_ilr <- function(object, base, weights) {
  D <- ncol(object)
  seq_parts <- seq_len(D - 1)
  parts <- colnames(object)

  ## Rotated and centered values
  y <- log(object, base = exp(1))
  ilr <- y %*% base

  ratio <- vapply(
    X = seq_parts,
    FUN = function(i, k) {
      sprintf("(%s)/%s", paste0(k[seq_len(i)], collapse = ","), k[i + 1])
    },
    FUN.VALUE = character(1),
    k = parts
  )
  colnames(ilr) <- paste0("Z", seq_parts)
  rownames(ilr) <- rownames(object)

  .ILR(
    ilr,
    parts = parts,
    ratio = ratio,
    order = seq_len(D),
    base = base,
    weights = weights,
    totals = totals(object)
  )
}

# Univariate ILR ===============================================================
#' @export
#' @rdname univariate_ilr
#' @aliases univariate_ilr,numeric-method
setMethod(
  f = "univariate_ilr",
  signature = c(object = "numeric"),
  definition = function(object) {
    sqrt(1 / 2) * log(object / (1 - object))
  }
)

#' @export
#' @rdname univariate_ilr
#' @aliases univariate_ilr,matrix-method
setMethod(
  f = "univariate_ilr",
  signature = c(object = "matrix"),
  definition = function(object) {
    apply(X = object, MARGIN = 2, FUN = univariate_ilr)
  }
)

# DATA TRANSFORMATION: PIVOT LOG RATIO
#' @include AllGenerics.R
NULL

# Pivot ========================================================================
#' @export
#' @rdname transform_plr
#' @aliases transform_plr,CompositionMatrix-method
setMethod(
  f = "transform_plr",
  signature = c(object = "CompositionMatrix"),
  definition = function(object, pivot = 1) {
    J <- ncol(object)
    parts <- colnames(object)

    ## Reorder
    pivot <- if (is.character(pivot)) which(parts == pivot) else as.integer(pivot)
    ordering <- c(pivot, which(pivot != seq_len(J)))
    parts <- parts[ordering]
    obj <- object[, ordering, drop = FALSE]

    x <- seq_len(J - 1)
    balances <- diag(sqrt((J - x) / (J - x + 1)))
    z <- 1 / matrix(data = seq_len(J) - J, nrow = J, ncol = J)
    z[lower.tri(z)] <- 0
    diag(z) <- 1
    z <- z[-nrow(z), ]

    H <- t(balances %*% z)
    plr <- log(obj, base = exp(1)) %*% H

    ratio <- vapply(
      X = seq_len(J - 1),
      FUN = function(i, parts) {
        j <- length(parts)
        sprintf("%s/(%s)", parts[1], paste0(parts[(i+1):j], collapse = ","))
      },
      FUN.VALUE = character(1),
      parts = parts
    )
    colnames(plr) <- paste0("Z", seq_len(J - 1))
    rownames(plr) <- rownames(object)

    .PLR(
      plr,
      parts = parts,
      ratio = ratio,
      order = order(ordering),
      base = H,
      weights = rep(1 / J, J),
      totals = totals(object)
    )
  }
)

#' @export
#' @rdname transform_plr
#' @aliases transform_plr,GroupedComposition-method
setMethod(
  f = "transform_plr",
  signature = c(object = "GroupedComposition"),
  definition = function(object, pivot = 1) {
    z <- methods::callNextMethod()
    .GroupedPLR(z, group_indices = group_indices(object),
                group_levels = group_levels(object))
  }
)

# CHEMISTRY
#' @include AllGenerics.R
NULL

is_chemical <- function(object) {
  pattern <- "^([A-Z]{1}[a-z]?[1-9]*)+$"
  z <- grepl(pattern, x = object)
  names(z) <- object
  z
}

#' @export
#' @rdname chemistry
#' @aliases is_oxide,character-method
setMethod(
  f = "is_oxide",
  signature = c(object = "character"),
  definition = function(object) {
    pattern <- "^[A-Z]{1}[a-z]?[1-9]*[O]{1}[1-9]*$"
    z <- grepl(pattern, x = object)
    names(z) <- object
    z
  }
)

#' @export
#' @rdname chemistry
#' @aliases is_oxide,CompositionMatrix-method
setMethod(
  f = "is_oxide",
  signature = c(object = "CompositionMatrix"),
  definition = function(object) {
    methods::callGeneric(colnames(object))
  }
)

#' @export
#' @rdname chemistry
#' @aliases is_element_major,CompositionMatrix-method
setMethod(
  f = "is_element_major",
  signature = c(object = "CompositionMatrix"),
  definition = function(object, min = 1 / 100, max = Inf) {
    .element_threshold(object, min = min, max = max)
  }
)

#' @export
#' @rdname chemistry
#' @aliases is_element_minor,CompositionMatrix-method
setMethod(
  f = "is_element_minor",
  signature = c(object = "CompositionMatrix"),
  definition = function(object, min = 0.1 / 100, max = 1 / 100) {
    .element_threshold(object, min = min, max = max)
  }
)

#' @export
#' @rdname chemistry
#' @aliases is_element_trace,CompositionMatrix-method
setMethod(
  f = "is_element_trace",
  signature = c(object = "CompositionMatrix"),
  definition = function(object, min = -Inf, max = 0.1 / 100) {
    .element_threshold(object, min = min, max = max)
  }
)

.element_threshold <- function(x, min = -Inf, max = Inf) {
  x <- mean(x)
  x >= min & x < max
}

# PAIRS
#' @include AllGenerics.R
NULL

# CompositionMatrix ============================================================
#' @export
#' @method pairs CompositionMatrix
pairs.CompositionMatrix <- function(x, margin = NULL, ...) {
  isopleuros::ternary_pairs(x, margin = margin, ...)
  invisible(x)
}

#' @export
#' @rdname pairs
#' @aliases pairs,CompositionMatrix-method
setMethod("pairs", c(x = "CompositionMatrix"), pairs.CompositionMatrix)

#' @export
#' @method pairs GroupedComposition
pairs.GroupedComposition <- function(x, ..., margin = NULL,
                                    palette_color = palette_color_discrete(),
                                    palette_symbol = palette_shape()) {
  ## Aesthetics
  lvl <- group_names(x)
  col <- palette_color(lvl)
  pch <- palette_symbol(lvl)

  isopleuros::ternary_pairs(x, margin = margin, col = col, pch = pch, ...)
  invisible(x)
}

#' @export
#' @rdname pairs
#' @aliases pairs,GroupedComposition-method
setMethod("pairs", c(x = "GroupedComposition"), pairs.GroupedComposition)

# SPLIT
#' @include AllGenerics.R
NULL

#' @export
#' @method split CompositionMatrix
split.CompositionMatrix <- function(x, f, drop = FALSE, ...) {
  lapply(
    X = split(x = seq_len(nrow(x)), f = f, drop = drop, sep = "_", ...),
    FUN = function(ind) x[ind, , drop = FALSE]
  )
}

#' @export
#' @rdname group_split
#' @aliases split,CompositionMatrix-method
setMethod("split", "CompositionMatrix", split.CompositionMatrix)

#' @export
#' @method split LogRatio
split.LogRatio <- function(x, f, drop = FALSE, ...) {
  lapply(
    X = split(x = seq_len(nrow(x)), f = f, drop = drop, sep = "_", ...),
    FUN = function(ind) x[ind, , drop = FALSE]
  )
}

#' @export
#' @rdname group_split
#' @aliases split,LogRatio-method
setMethod("split", "LogRatio", split.LogRatio)

#' @export
#' @rdname group_split
#' @aliases group_split,GroupedComposition-method
setMethod(
  f = "group_split",
  signature = "GroupedComposition",
  definition = function(object, ...) {
    lapply(
      X = group_rows(object),
      FUN = function(ind) ungroup(object[ind, , drop = FALSE])
    )
  }
)

#' @export
#' @rdname group_split
#' @aliases group_split,GroupedLogRatio-method
setMethod(
  f = "group_split",
  signature = "GroupedLogRatio",
  definition = function(object, ...) {
    lapply(
      X = group_rows(object),
      FUN = function(ind) ungroup(object[ind, , drop = FALSE])
    )
  }
)

# VARIANCE
#' @include AllGenerics.R
NULL

# Variance =====================================================================
#' @export
#' @rdname variance
#' @aliases variance,LogRatio-method
setMethod(
  f = "variance",
  signature = c("LogRatio"),
  definition = function(x, row_weights = NULL, column_weights = TRUE) {
    n <- nrow(x)
    m <- ncol(x)

    w_row <- rep(1 / n, n)
    if (length(row_weights) == n) {
      arkhe::assert_length(row_weights, n)
      arkhe::assert_positive(row_weights, strict = FALSE)
      w_row <- row_weights / sum(row_weights) # Sum up to 1
    }

    w_col <- if (isTRUE(column_weights)) weights(x) else rep(1 / m, m)
    if (length(column_weights) == m) {
      arkhe::assert_length(column_weights, m)
      arkhe::assert_positive(column_weights, strict = FALSE)
      w_col <- column_weights / sum(column_weights) # Sum up to 1
    }

    z <- sweep(x, MARGIN = 2, STATS = colSums(x * w_row), FUN = "-")
    z <- colSums(diag(w_row) %*% z^2 %*% diag(w_col))
    names(z) <- colnames(x)
    z
  }
)

# Total variance ===============================================================
#' @export
#' @describeIn variance_total The total variance of compositional data is the
#'  trace of the [centred log-ratio covariance][covariance()] matrix
#'  (i.e. *totvar1* in Aitchison 1997).
#' @aliases variance_total,CompositionMatrix-method
setMethod(
  f = "variance_total",
  signature = c("CompositionMatrix"),
  definition = function(x, sd = FALSE) {
    z <- sum(diag(covariance(x, center = TRUE)))
    if (sd) z <- sqrt((1 / (ncol(x) - 1)) * z)
    z
  }
)

#' @export
#' @describeIn variance_total Computes the total log-ratio variance. This is
#'  identical to the weighted sum-of-squared distances between samples
#'  (i.e. *totvar2* in Aitchison 1997).
#' @aliases variance_total,LogRatio-method
setMethod(
  f = "variance_total",
  signature = c("LogRatio"),
  definition = function(x, row_weights = NULL, column_weights = TRUE) {
    sum(variance(x, row_weights = row_weights, column_weights = column_weights))
  }
)


# STATISTICS: QUANTILE
#' @include AllGenerics.R
NULL

# Quantile =====================================================================
#' @export
#' @method quantile CompositionMatrix
quantile.CompositionMatrix <- function(x, ..., probs = seq(0, 1, 0.25),
                                       na.rm = FALSE, names = TRUE) {
  apply(
    X = x,
    MARGIN = 2,
    FUN = stats::quantile,
    probs = probs,
    na.rm = na.rm,
    names = names,
    ...
  )
}

#' @export
#' @rdname quantile
#' @aliases quantile,CompositionMatrix-method
setMethod("quantile", "CompositionMatrix", quantile.CompositionMatrix)

# AGGREGATE
#' @include AllGenerics.R
NULL

# Aggregate ====================================================================
#' @export
#' @method aggregate CompositionMatrix
aggregate.CompositionMatrix <- function(x, by, FUN, ...,
                                        simplify = TRUE, drop = TRUE) {
  x <- group(x, by = by, drop_levels = drop, verbose = FALSE)
  aggregate(x, FUN, ..., simplify = simplify)
}

#' @export
#' @rdname aggregate
#' @aliases aggregate,CompositionMatrix-method
setMethod("aggregate", "CompositionMatrix", aggregate.CompositionMatrix)

#' @export
#' @method aggregate GroupedComposition
aggregate.GroupedComposition <- function(x, FUN, ..., simplify = TRUE) {
  ## Grouping
  aggr <- lapply(
    X = group_rows(x),
    FUN = function(i, data, fun, ...) fun(data[i, , drop = FALSE], ...),
    data = x,
    fun = FUN,
    ...
  )

  has_dim <- vapply(
    X = aggr,
    FUN = function(x) !is.null(nrow(x)) && nrow(x) > 1,
    FUN.VALUE = logical(1)
  )

  if (any(has_dim) || !simplify) return(aggr)
  aggr <- do.call(rbind, aggr)
  rownames(aggr) <- group_levels(x)
  aggr
}

#' @export
#' @rdname aggregate
#' @aliases aggregate,GroupedComposition-method
setMethod("aggregate", "GroupedComposition", aggregate.GroupedComposition)

# SCALE
#' @include AllGenerics.R
NULL

# Scale ========================================================================
#' @export
#' @method scale CompositionMatrix
scale.CompositionMatrix <- function(x, center = TRUE, scale = TRUE) {
  if (isFALSE(center) & isFALSE(scale)) return(x)

  y <- x
  if (!isFALSE(center)) {
    if (isTRUE(center)) center <- mean(x)
    arkhe::assert_type(center, "numeric")
    arkhe::assert_length(center, NCOL(x))

    y <- perturbation(y, 1 / center)
  }

  if (!isFALSE(scale)) {
    if (isTRUE(scale)) scale <- sqrt(mean(diag(covariance(x, center = TRUE))))
    arkhe::assert_type(scale, "numeric")

    y <- powering(y, 1 / scale)
  }

  y
}

#' @export
#' @rdname scale
#' @aliases scale,CompositionMatrix-method
setMethod("scale", "CompositionMatrix", scale.CompositionMatrix)

# COVARIANCE
#' @include AllGenerics.R
NULL

# Covariance ===================================================================
#' @export
#' @rdname covariance
#' @aliases covariance,CompositionMatrix-method
setMethod(
  f = "covariance",
  signature = c(x = "CompositionMatrix"),
  definition = function(x, center = TRUE, method = "pearson") {
    x <- if (center) transform_clr(x) else transform_alr(x)
    methods::callGeneric(x = x, method = method)
  }
)

#' @export
#' @describeIn covariance Computes the log-ratio covariance matrix
#'  (Aitchison 1986, definition 4.5).
#' @aliases covariance,ALR-method
setMethod(
  f = "covariance",
  signature = c(x = "ALR"),
  definition = function(x, method = "pearson") {
    stats::cov(x, method = method)
  }
)

#' @export
#' @describeIn covariance Computes the centered log-ratio covariance matrix
#'  (Aitchison 1986, definition 4.6).
#' @aliases covariance,ALR-method
setMethod(
  f = "covariance",
  signature = c(x = "CLR"),
  definition = function(x, method = "pearson") {
    stats::cov(x, method = method)
  }
)

# SHOW
#' @include AllGenerics.R
NULL

setMethod(
  f = "show",
  signature = "CompositionMatrix",
  definition = function(object) {
    m <- nrow(object)
    p <- ncol(object)
    mtx <- methods::as(object, "matrix")

    txt_dim <- sprintf("%d x %d", m, p)
    txt_mtx <- sprintf("<%s: %s>", class(object), txt_dim)

    cat(
      txt_mtx,
      utils::capture.output(mtx),
      sep = "\n"
    )
    invisible(object)
  }
)

setMethod(
  f = "show",
  signature = "LogRatio",
  definition = function(object) {
    m <- nrow(object)
    p <- ncol(object)
    mtx <- methods::as(object, "matrix")

    txt_dim <- sprintf("%d x %d", m, p)
    txt_mtx <- sprintf("<%s: %s>", get_transformation(object), txt_dim)

    cat(
      txt_mtx,
      utils::capture.output(mtx),
      sep = "\n"
    )
    invisible(object)
  }
)

1		# COERCION
2		#' @include AllGenerics.R
3		NULL
4
5		# To CompositionMatrix =========================================================
6		#' @export
7		#' @rdname as_composition
8		#' @aliases as_composition,numeric-method
9		setMethod(
10		f = "as_composition",
11		signature = c(from = "numeric"),
12		definition = function(from) {
13	17x	from <- matrix(data = from, nrow = 1, ncol = length(from))
14	17x	methods::callGeneric(from)
15		}
16		)
17
18		#' @export
19		#' @rdname as_composition
20		#' @aliases as_composition,matrix-method
21		setMethod(
22		f = "as_composition",
23		signature = c(from = "matrix"),
24		definition = function(from) {
25		## Validation
26	33x	arkhe::assert_type(from, "numeric")
27
28		## Make row/column names
29	33x	lab <- make_names(x = NULL, n = nrow(from), prefix = "S")
30	33x	rownames(from) <- if (has_rownames(from)) rownames(from) else lab
31	33x	colnames(from) <- make_names(x = colnames(from), n = ncol(from), prefix = "V")
32
33		## Close
34	33x	totals <- rowSums(from, na.rm = TRUE)
35	33x	from <- from / totals
36
37	33x	.CompositionMatrix(from, totals = unname(totals))
38		}
39		)
40
41		#' @export
42		#' @rdname as_composition
43		#' @aliases as_composition,data.frame-method
44		setMethod(
45		f = "as_composition",
46		signature = c(from = "data.frame"),
47		definition = function(from, parts = NULL, groups = NULL, autodetect = TRUE,
48		verbose = getOption("nexus.verbose")) {
49		## Clean row/column names
50	18x	lab <- make_names(x = NULL, n = nrow(from), prefix = "S")
51	18x	rownames(from) <- if (has_rownames(from)) rownames(from) else lab
52	18x	colnames(from) <- make_names(x = colnames(from), n = ncol(from), prefix = "V")
53
54		## Remove non-numeric columns
55	18x	if (is.null(parts)) {
56	17x	if (isTRUE(autodetect)) {
57	17x	parts <- arkhe::detect(from, f = is.numeric, margin = 2)
58	17x	if (isTRUE(verbose)) {
59	!	n <- sum(parts)
60	!	what <- ngettext(n, "Found %g part (%s)", "Found %g parts (%s)")
61	!	cols <- paste0(colnames(from)[parts], collapse = ", ")
62	!	message(sprintf(what, n, cols))
63		}
64		} else {
65	!	arkhe::assert_filled(parts)
66		}
67		} else {
68	1x	if (is.numeric(parts)) parts <- seq_len(ncol(from)) %in% parts
69	!	if (is.character(parts)) parts <- colnames(from) %in% parts
70		}
71	18x	coda <- from[, parts, drop = FALSE]
72	18x	arkhe::assert_filled(coda)
73
74		## Build matrix
75	18x	coda <- data.matrix(coda, rownames.force = NA)
76	18x	totals <- rowSums(coda, na.rm = TRUE)
77	18x	coda <- coda / totals
78
79	18x	z <- .CompositionMatrix(coda, totals = unname(totals))
80	17x	if (is.null(groups)) return(z)
81
82		## Group names
83	1x	grp <- from[, groups, drop = FALSE]
84	1x	group(z, by = grp, verbose = verbose)
85		}
86		)
87
88		# To amounts ===================================================================
89		#' @export
90		#' @rdname as_amounts
91		#' @aliases as_amounts,CompositionMatrix-method
92		setMethod(
93		f = "as_amounts",
94		signature = c(from = "CompositionMatrix"),
95		definition = function(from) {
96	2x	methods::as(from, "matrix") * totals(from)
97		}
98		)
99
100		# To data.frame ================================================================
101		#' @method as.data.frame CompositionMatrix
102		#' @export
103		as.data.frame.CompositionMatrix <- function(x, row.names = rownames(x),
104		optional = FALSE, ...) {
105	!	as.data.frame(methods::as(x, "matrix"), row.names = row.names, optional = optional)
106		}
107
108		#' @export
109		#' @rdname as.data.frame
110		#' @aliases as.data.frame,CompositionMatrix-method
111		setMethod("as.data.frame", "CompositionMatrix", as.data.frame.CompositionMatrix)
112
113		#' @method as.data.frame GroupedComposition
114		#' @export
115		as.data.frame.GroupedComposition <- function(x, row.names = rownames(x),
116		optional = FALSE, ...,
117		group_var = ".group",
118		group_after = 0) {
119	1x	z <- data.frame(
120	1x	methods::as(x, "matrix"),
121	1x	row.names = row.names,
122	1x	check.names = !optional
123		)
124	1x	arkhe::append_column(z, group_names(x), after = group_after, var = group_var)
125		}
126
127		#' @export
128		#' @rdname as.data.frame
129		#' @aliases as.data.frame,GroupedComposition-method
130		setMethod("as.data.frame", "GroupedComposition", as.data.frame.GroupedComposition)
131
132		#' @method as.data.frame LogRatio
133		#' @export
134		as.data.frame.LogRatio <- function(x, row.names = rownames(x),
135		optional = FALSE, ...) {
136	!	z <- as.data.frame(methods::as(x, "matrix"), row.names = row.names, optional = optional)
137
138	!	lab <- labels(x)
139	!	for (j in seq_len(ncol(z))) attr(z[[j]], "label") <- lab[[j]]
140	!	z
141		}
142
143		#' @export
144		#' @rdname as.data.frame
145		#' @aliases as.data.frame,LogRatio-method
146		setMethod("as.data.frame", "LogRatio", as.data.frame.LogRatio)
147
148		#' @method as.data.frame GroupedLogRatio
149		#' @export
150		as.data.frame.GroupedLogRatio <- function(x, row.names = rownames(x),
151		optional = FALSE, ...,
152		group_var = ".group",
153		group_after = 0) {
154	!	z <- data.frame(
155	!	methods::as(x, "matrix"),
156	!	row.names = row.names,
157	!	check.names = !optional
158		)
159
160	!	lab <- labels(x)
161	!	for (j in seq_len(ncol(z))) attr(z[[j]], "label") <- lab[[j]]
162
163	!	arkhe::append_column(z, group_names(x), after = group_after, var = group_var)
164		}
165
166		#' @export
167		#' @rdname as.data.frame
168		#' @aliases as.data.frame,GroupedLogRatio-method
169		setMethod("as.data.frame", "GroupedLogRatio", as.data.frame.GroupedLogRatio)
170
171		#' @method as.data.frame OutlierIndex
172		#' @export
173		as.data.frame.OutlierIndex <- function(x, row.names = rownames(x),
174		optional = FALSE, ...) {
175	!	as.data.frame(x@standard, row.names = row.names, optional = optional)
176		}
177
178		#' @export
179		#' @rdname as.data.frame
180		#' @aliases as.data.frame,OutlierIndex-method
181		setMethod("as.data.frame", "OutlierIndex", as.data.frame.OutlierIndex)

1		# OPERATORS
2		#' @include AllGenerics.R
3		NULL
4
5		# Closure ======================================================================
6		#' @export
7		#' @rdname closure
8		#' @aliases closure,numeric-method
9		setMethod(
10		f = "closure",
11		signature = c(x = "numeric"),
12		definition = function(x, total = 1, na.rm = TRUE) {
13	41x	x * total / sum(x, na.rm = na.rm)
14		}
15		)
16
17		#' @export
18		#' @rdname closure
19		#' @aliases closure,matrix-method
20		setMethod(
21		f = "closure",
22		signature = c(x = "matrix"),
23		definition = function(x, total = 1, na.rm = TRUE) {
24	1x	x * total / rowSums(x, na.rm = na.rm)
25		}
26		)
27
28		# Perturbation =================================================================
29		#' @export
30		#' @rdname arithmetic
31		#' @aliases `%perturbe%`,CompositionMatrix,CompositionMatrix-method
32		setMethod(
33		f = "%perturbe%",
34		signature = c(x = "CompositionMatrix", y = "CompositionMatrix"),
35		definition = function(x, y) {
36	5x	arkhe::assert_dimensions(y, dim(x))
37
38	!	if (all(x <= 0)) x <- 1 / x
39	1x	if (all(y <= 0)) y <- 1 / y
40
41	5x	z <- x * y
42	5x	z <- as_composition(z)
43	5x	rownames(z) <- rownames(x)
44
45	5x	if (is_grouped(x)) {
46	!	x@group_indices <- group_indices(x)
47	!	x@group_levels <- group_levels(x)
48		}
49
50	5x	z
51		}
52		)
53
54		#' @export
55		#' @rdname perturbation
56		#' @aliases perturbation,numeric,numeric-method
57		setMethod(
58		f = "perturbation",
59		signature = c(x = "numeric", y = "numeric"),
60		definition = function(x, y) {
61	1x	arkhe::assert_length(y, length(x))
62	1x	closure(x * y)
63		}
64		)
65
66		#' @export
67		#' @rdname perturbation
68		#' @aliases perturbation,CompositionMatrix,numeric-method
69		setMethod(
70		f = "perturbation",
71		signature = c(x = "CompositionMatrix", y = "numeric"),
72		definition = function(x, y) {
73	2x	y <- matrix(data = y, nrow = nrow(x), ncol = length(y), byrow = TRUE)
74	2x	x %perturbe% as_composition(y)
75		}
76		)
77
78		#' @export
79		#' @rdname perturbation
80		#' @aliases perturbation,CompositionMatrix,CompositionMatrix-method
81		setMethod(
82		f = "perturbation",
83		signature = c(x = "CompositionMatrix", y = "matrix"),
84		definition = function(x, y) {
85	1x	x %perturbe% as_composition(y)
86		}
87		)
88
89		# Powering =====================================================================
90		#' @export
91		#' @rdname arithmetic
92		#' @aliases `%power%`,CompositionMatrix,numeric-method
93		setMethod(
94		f = "%power%",
95		signature = c(x = "CompositionMatrix", y = "numeric"),
96		definition = function(x, y) {
97	3x	arkhe::assert_length(y, 1L)
98
99	3x	z <- x ^ y
100	3x	z <- as_composition(z)
101	3x	rownames(z) <- rownames(x)
102
103	3x	if (is_grouped(x)) {
104	!	x@group_indices <- group_indices(x)
105	!	x@group_levels <- group_levels(x)
106		}
107
108	3x	z
109		}
110		)
111
112		#' @export
113		#' @rdname arithmetic
114		#' @aliases `%power%`,numeric,CompositionMatrix-method
115		setMethod(
116		f = "%power%",
117		signature = c(x = "numeric", y = "CompositionMatrix"),
118		definition = function(x, y) {
119	1x	methods::callGeneric(x = y, y = x)
120		}
121		)
122
123		#' @export
124		#' @rdname powering
125		#' @aliases powering,numeric,numeric-method
126		setMethod(
127		f = "powering",
128		signature = c(x = "numeric", a = "numeric"),
129		definition = function(x, a) {
130	1x	arkhe::assert_length(a, 1L)
131	1x	closure(x ^ a)
132		}
133		)
134
135		#' @export
136		#' @rdname powering
137		#' @aliases powering,numeric,numeric-method
138		setMethod(
139		f = "powering",
140		signature = c(x = "CompositionMatrix", a = "numeric"),
141		definition = function(x, a) {
142	1x	x %power% a
143		}
144		)
145
146		# Scalar product ===============================================================
147		#' @export
148		#' @rdname scalar
149		#' @aliases scalar,numeric,numeric-method
150		setMethod(
151		f = "scalar",
152		signature = c(x = "numeric", y = "numeric"),
153		definition = function(x, y) {
154	302x	n <- length(x)
155	302x	arkhe::assert_length(y, n)
156
157	302x	D <- seq_len(n)
158	302x	z <- 0
159	302x	for (i in D) {
160	1506x	j <- utils::tail(D, -i)
161	1506x	z <- z + sum(log(x[i] / x[j]) * log(y[i] / y[j]))
162		}
163	302x	(1 / n) * z
164		}
165		)
166
167		#' @export
168		#' @rdname scalar
169		#' @aliases scalar,CompositionMatrix,CompositionMatrix-method
170		setMethod(
171		f = "scalar",
172		signature = c(x = "CompositionMatrix", y = "CompositionMatrix"),
173		definition = function(x, y) {
174	1x	arkhe::assert_dimensions(y, dim(x))
175	1x	m <- nrow(x)
176
177	1x	z <- numeric(m)
178	1x	for (i in seq_len(m)) {
179	1x	z[i] <- scalar(x[i, , drop = TRUE], y[i, , drop = TRUE])
180		}
181	1x	z
182		}
183		)
184
185		#' Norm of a Composition
186		#'
187		#' @param x A [`CompositionMatrix-class`] object.
188		#' @return A [`numeric`] vector.
189		#' @keywords internal
190		#' @noRd
191		norm <- function(x) {
192	300x	sqrt(scalar(x, x))
193		}
194
195		#' Aitchison Distance
196		#'
197		#' @param x A [`CompositionMatrix-class`] object.
198		#' @param diag A [`logical`] scalar: should the diagonal of the distance matrix
199		#' be printed?
200		#' @param upper A [`logical`] scalar: should the upper triangle of the distance
201		#' matrix be printed?
202		#' @return A [`dist`] object.
203		#' @keywords internal
204		#' @noRd
205		aitchison <- function(x, diag = FALSE, upper = FALSE) {
206	1x	m <- nrow(x)
207	1x	spl <- rownames(x)
208
209	1x	d <- utils::combn(
210	1x	x = seq_len(m),
211	1x	m = 2,
212	1x	FUN = function(i, coda) {
213	300x	x <- coda[i[1], , drop = TRUE]
214	300x	y <- coda[i[2], , drop = TRUE]
215	300x	norm(x / y)
216		},
217	1x	coda = x
218		)
219
220		## Matrix of results
221	1x	mtx <- matrix(data = 0, nrow = m, ncol = m, dimnames = list(spl, spl))
222	1x	mtx[lower.tri(mtx, diag = FALSE)] <- d
223	1x	mtx <- t(mtx)
224	1x	mtx[lower.tri(mtx, diag = FALSE)] <- d
225
226	1x	stats::as.dist(mtx, diag = diag, upper = upper)
227		}

1		# STATISTICS: MEAN
2		#' @include AllGenerics.R
3		NULL
4
5		# Mean =========================================================================
6		#' @export
7		#' @method mean CompositionMatrix
8		mean.CompositionMatrix <- function(x, ..., ignore_na = FALSE, ignore_zero = TRUE) {
9	35x	m <- apply(
10	35x	X = x,
11	35x	MARGIN = 2,
12	35x	FUN = gmean,
13	35x	ignore_na = ignore_na,
14	35x	ignore_zero = ignore_zero,
15	35x	simplify = TRUE
16		)
17	35x	m <- closure(m)
18	35x	names(m) <- colnames(x)
19	35x	m
20		}
21
22		#' @export
23		#' @rdname mean
24		#' @aliases mean,CompositionMatrix-method
25		setMethod("mean", "CompositionMatrix", mean.CompositionMatrix)
26
27		#' Geometric Mean
28		#'
29		#' @param x A [`numeric`] vector.
30		#' @param trim A length-one [`numeric`] vector specifying the fraction (0 to 0.5)
31		#' of observations to be trimmed from each end of `x` before the mean is
32		#' computed.
33		#' @param ignore_na A [`logical`] scalar: should [missing values][NA] be
34		#' stripped before the computation proceeds?
35		#' @param ignore_zero A [`logical`] scalar: should zeros be stripped before the
36		#' computation proceeds?
37		#' @return A [`numeric`] vector.
38		#' @keywords internal
39		gmean <- function(x, trim = 0, ignore_na = FALSE, ignore_zero = TRUE) {
40	10x	if (ignore_na) x <- x[is.finite(x)]
41	235x	if (ignore_zero) x <- x[x > 0]
42	237x	exp(mean(log(unclass(x)), trim = trim))
43		}

1		# GROUPS
2		#' @include AllGenerics.R
3		NULL
4
5		.group_extract <- function(object, which) {
6		## Validation
7	1x	arkhe::assert_type(which, "character")
8	!	if (!any_assigned(object)) stop("No group is defined.", call. = FALSE)
9
10	1x	ok <- group_names(object) %in% which
11	1x	if (!any(ok)) {
12	!	msg <- ngettext(length(which), "No sample belongs to the group: %s.",
13	!	"No sample belongs to the groups: %s.")
14	!	message(sprintf(msg, paste0(dQuote(which), collapse = ", ")))
15	!	return(object)
16		}
17
18	1x	object[ok, , drop = FALSE]
19		}
20
21		#' @export
22		#' @rdname group_subset
23		#' @aliases group_subset,GroupedComposition-method
24		setMethod(
25		f = "group_subset",
26		signature = c("GroupedComposition"),
27		definition = .group_extract
28		)
29
30		#' @export
31		#' @rdname group_subset
32		#' @aliases group_subset,GroupedLogRatio-method
33		setMethod(
34		f = "group_subset",
35		signature = c("GroupedLogRatio"),
36		definition = .group_extract
37		)
38
39		# Groups =======================================================================
40		#' Compute Groups
41		#'
42		#' @param x A ([`list`] of) [`factors`] for which interaction is to be computed.
43		#' @param drop_levels A [`logical`] scalar: should unused factor levels be
44		#' dropped?
45		#' @param allow_na A [`logical`] scalar: should `NA` be considered an extra
46		#' level?
47		#' @return A [`factor`] vector.
48		#' @keywords internal
49		#' @noRd
50		compute_groups <- function(x, drop_levels = TRUE, allow_na = TRUE) {
51	19x	if (!is.list(x)) x <- list(x)
52	21x	x <- interaction(x, sep = "_")
53	21x	if (drop_levels) x <- droplevels(x)
54	21x	if (allow_na) x <- addNA(x, ifany = TRUE)
55
56	21x	x
57		}
58
59		#' @export
60		#' @rdname group
61		#' @aliases group,CompositionMatrix-method
62		setMethod(
63		f = "group",
64		signature = "CompositionMatrix",
65		definition = function(object, by, verbose = getOption("nexus.verbose"), ...) {
66		## Compute groups
67	21x	by <- compute_groups(by, ...)
68
69		## Validation
70	21x	arkhe::assert_length(by, nrow(object))
71	20x	if (nlevels(by) == 0) {
72	!	stop(tr_("Nothing to group by."), call. = FALSE)
73		}
74	20x	if (isTRUE(verbose)) {
75	1x	if (nlevels(by) == nrow(object)) {
76	!	message(tr_("As many groups as individuals."))
77		}
78
79	1x	n <- nlevels(by)
80	1x	what <- ngettext(n, "Found %g group (%s)", "Found %g groups (%s)")
81	1x	grp <- paste0(levels(by), collapse = ", ")
82	1x	message(sprintf(what, n, grp))
83		}
84
85	20x	.GroupedComposition(
86	20x	object,
87	20x	group_indices = as.integer(by),
88	20x	group_levels = levels(by)
89		)
90		}
91		)
92
93		#' @export
94		#' @rdname group
95		#' @aliases group,GroupedComposition-method
96		setMethod(
97		f = "group",
98		signature = "GroupedComposition",
99		definition = function(object, by, add = FALSE,
100		verbose = getOption("nexus.verbose"), ...) {
101		## Compute groups
102	4x	if (isTRUE(add)) {
103	1x	if (!is.list(by)) by <- list(by)
104	1x	by <- c(list(group_factor(object)), by)
105		}
106	4x	methods::callNextMethod(object, by = by, verbose = verbose, ...)
107		}
108		)
109
110		#' @export
111		#' @rdname group
112		#' @aliases ungroup,GroupedComposition-method
113		setMethod(
114		f = "ungroup",
115		signature = "GroupedComposition",
116		definition = function(object) {
117	21x	methods::as(object, "CompositionMatrix", strict = TRUE)
118		}
119		)
120
121		#' @export
122		#' @rdname group
123		#' @aliases ungroup,GroupedLR-method
124		setMethod(
125		f = "ungroup",
126		signature = "GroupedLR",
127		definition = function(object) {
128	!	methods::as(object, "LR", strict = TRUE)
129		}
130		)
131
132		#' @export
133		#' @rdname group
134		#' @aliases ungroup,GroupedCLR-method
135		setMethod(
136		f = "ungroup",
137		signature = "GroupedCLR",
138		definition = function(object) {
139	!	methods::as(object, "CLR", strict = TRUE)
140		}
141		)
142
143		#' @export
144		#' @rdname group
145		#' @aliases ungroup,GroupedALR-method
146		setMethod(
147		f = "ungroup",
148		signature = "GroupedALR",
149		definition = function(object) {
150	!	methods::as(object, "ALR", strict = TRUE)
151		}
152		)
153
154		#' @export
155		#' @rdname group
156		#' @aliases ungroup,GroupedILR-method
157		setMethod(
158		f = "ungroup",
159		signature = "GroupedILR",
160		definition = function(object) {
161	!	methods::as(object, "ILR", strict = TRUE)
162		}
163		)
164
165		#' @export
166		#' @rdname group
167		#' @aliases ungroup,GroupedPLR-method
168		setMethod(
169		f = "ungroup",
170		signature = "GroupedPLR",
171		definition = function(object) {
172	!	methods::as(object, "PLR", strict = TRUE)
173		}
174		)
175
176		# Metadata =====================================================================
177		#' @export
178		#' @describeIn group_metadata returns a [`character`] vector giving the group
179		#' names.
180		#' @aliases group_levels,ReferenceGroups-method
181		setMethod(
182		f = "group_levels",
183		signature = "ReferenceGroups",
184	225x	definition = function(object) object@group_levels
185		)
186
187		#' @export
188		#' @describeIn group_metadata returns a [`character`] vector giving the name of
189		#' the group that each observation belongs to.
190		#' @aliases group_names,ReferenceGroups-method
191		setMethod(
192		f = "group_names",
193		signature = "ReferenceGroups",
194	112x	definition = function(object) group_levels(object)[group_indices(object)]
195		)
196
197		group_factor <- function(object) {
198	100x	factor(
199	100x	x = group_names(object),
200	100x	levels = group_levels(object),
201	100x	exclude = NULL
202		)
203		}
204
205		#' @export
206		#' @describeIn group_metadata returns an [`integer`] vector giving the group
207		#' that each value belongs to.
208		#' @aliases group_indices,ReferenceGroups-method
209		setMethod(
210		f = "group_indices",
211		signature = "ReferenceGroups",
212	114x	definition = function(object) object@group_indices
213		)
214
215		#' @export
216		#' @describeIn group_metadata returns a `list` of [`integer`] vectors giving the
217		#' observation that each group contains.
218		#' @aliases group_rows,ReferenceGroups-method
219		setMethod(
220		f = "group_rows",
221		signature = "ReferenceGroups",
222		definition = function(object) {
223	18x	i <- group_factor(object)
224	18x	split(seq_along(i), f = i)
225		}
226		)
227
228		#' @export
229		#' @describeIn group_metadata gives the total number of groups.
230		#' @aliases group_length,ReferenceGroups-method
231		setMethod(
232		f = "group_length",
233		signature = "ReferenceGroups",
234	5x	definition = function(object) length(group_levels(object))
235		)
236
237		#' @export
238		#' @describeIn group_metadata gives the size of each group.
239		#' @aliases group_size,ReferenceGroups-method
240		setMethod(
241		f = "group_size",
242		signature = "ReferenceGroups",
243	5x	definition = function(object) lengths(group_rows(object))
244		)
245
246		#' @export
247		#' @describeIn group_metadata returns a [`logical`] vector specifying whether or
248		#' not an observation belongs to a group.
249		#' @aliases is_assigned,ReferenceGroups-method
250		setMethod(
251		f = "is_assigned",
252		signature = "ReferenceGroups",
253	5x	definition = function(object) !is.na(group_names(object))
254		)
255
256		#' @export
257		#' @describeIn group_metadata returns an [`logical`] scalar specifying if any
258		#' observation belongs to a group.
259		#' @aliases any_assigned,ReferenceGroups-method
260		setMethod(
261		f = "any_assigned",
262		signature = "ReferenceGroups",
263	3x	definition = function(object) any(is_assigned(object))
264		)
265
266		#' @export
267		#' @describeIn group_metadata returns an [`logical`] scalar specifying if all
268		#' observations belong to a group.
269		#' @aliases all_assigned,ReferenceGroups-method
270		setMethod(
271		f = "all_assigned",
272		signature = "ReferenceGroups",
273	1x	definition = function(object) all(is_assigned(object))
274		)

1		# HISTOGRAM
2		#' @include AllGenerics.R
3		NULL
4
5		# CompositionMatrix ============================================================
6		#' @export
7		#' @method hist CompositionMatrix
8		hist.CompositionMatrix <- function(x, ..., select = 1,
9		breaks = "Sturges",
10		freq = FALSE, labels = FALSE,
11		main = NULL, sub = NULL,
12		ann = graphics::par("ann"),
13		axes = TRUE, frame.plot = axes) {
14		## Validation
15	1x	if (is.character(select)) select <- match(select, labels(x))
16	2x	arkhe::assert_length(select, 1)
17
18		## Graphical parameters
19	2x	cex.lab <- list(...)$cex.lab %\|\|% graphics::par("cex.lab")
20	2x	col.lab <- list(...)$col.lab %\|\|% graphics::par("col.lab")
21	2x	font.lab <- list(...)$font.lab %\|\|% graphics::par("font.lab")
22	2x	cex.main <- list(...)$cex.main %\|\|% graphics::par("cex.main")
23	2x	col.main <- list(...)$col.main %\|\|% graphics::par("col.main")
24	2x	font.main <- list(...)$font.main %\|\|% graphics::par("font.main")
25
26		## Compute univariate ilr transformation
27	2x	z <- univariate_ilr(x)
28
29		## Select one compositonal part
30	2x	xi <- x[, select, drop = TRUE]
31	2x	zi <- z[, select, drop = TRUE]
32
33		## Compute axis in percent
34	2x	lab_i <- pretty(xi, n = 6)
35	2x	lab_i <- lab_i[lab_i > 0]
36	2x	at_i <- univariate_ilr(lab_i)
37
38		## Plot histogram
39	2x	h <- graphics::hist(x = zi, breaks = breaks, plot = FALSE)
40	2x	xlim <- range(at_i, h$breaks, finite = TRUE)
41	2x	plot(h, freq = freq, xlim = xlim, labels = labels, ...,
42	2x	main = main, sub = sub, xlab = NULL, ylab = NULL, axes = FALSE)
43
44		## Construct axis
45	2x	if (axes) {
46	2x	graphics::axis(side = 1, xpd = NA, las = 1)
47	2x	graphics::axis(side = 3, at = at_i, labels = label_percent(lab_i),
48	2x	xpd = NA, las = 1)
49	2x	graphics::axis(side = 2, xpd = NA, las = 1)
50		}
51
52		## Plot frame
53	2x	if (frame.plot) {
54	2x	graphics::box()
55		}
56
57		## Add annotation
58	2x	if (ann) {
59	2x	xlab <- labels(x)[select]
60	2x	ylab <- tr_("Frequency")
61	2x	graphics::mtext(sprintf("ilr(%s)", xlab), side = 1, line = 3,
62	2x	cex = cex.lab, col = col.lab, font = font.lab)
63	2x	graphics::mtext(sprintf("%s %%", xlab), side = 3, line = 3,
64	2x	cex = cex.lab, col = col.lab, font = font.lab)
65	2x	graphics::mtext(ylab, side = 2, line = 3,
66	2x	cex = cex.lab, col = col.lab, font = font.lab)
67		}
68
69	2x	invisible(x)
70		}
71
72		#' @export
73		#' @rdname hist
74		#' @aliases hist,CompositionMatrix-method
75		setMethod("hist", c(x = "CompositionMatrix"), hist.CompositionMatrix)

1		# ACCESSORS
2		#' @include AllGenerics.R
3		NULL
4
5		# Extract ======================================================================
6		.subscript1 <- function(x, i) {
7	5x	x@.Data[i]
8		}
9
10		.subscript2 <- function(x, i, j, drop) {
11		## Rows
12	203x	if (missing(i)) i <- seq_len(nrow(x))
13	!	if (is.character(i)) i <- match(i, dimnames(x)[1L])
14	867x	totals <- totals(x)[i]
15
16		## Columns
17	658x	if (missing(j)) j <- seq_len(ncol(x))
18
19		## Subset
20	867x	z <- x@.Data[i, j, drop = drop]
21	747x	if (drop) return(z)
22
23		## /!\ Subcomposition /!\
24		# if (ncol(z) < ncol(x)) {
25		# tot <- rowSums(z, na.rm = TRUE)
26		# totals <- totals * tot
27		# z <- z / tot
28		# }
29
30	120x	if (is_grouped(x)) {
31	66x	g <- droplevels(group_factor(x)[i])
32	66x	methods::initialize(
33	66x	x, z,
34	66x	totals = totals,
35	66x	group_indices = as.integer(g),
36	66x	group_levels = levels(g)
37		)
38		} else {
39	54x	methods::initialize(x, z, totals = totals)
40		}
41		}
42
43		wrong_dimensions <- function(i, j) {
44	7x	msg <- sprintf("M[%s%s]: incorrect number of dimensions.", i, j)
45	7x	stop(msg, call. = FALSE)
46		}
47
48		## CompositionMatrix -----------------------------------------------------------
49		#' @export
50		#' @rdname subset
51		#' @aliases [,CompositionMatrix,missing,missing,missing-method
52		setMethod(
53		f = "[",
54		signature = c(x = "CompositionMatrix", i = "missing", j = "missing", drop = "missing"),
55		definition = function(x, i, j, ..., drop) {
56	5x	na <- nargs()
57	2x	if (na == 2L) return(x) # x[]
58	2x	if (na == 3L) return(x) # x[, ]
59	1x	wrong_dimensions(".", ".") # x[, , ], etc.
60		}
61		)
62
63		#' @export
64		#' @rdname subset
65		#' @aliases [,CompositionMatrix,missing,missing,logical-method
66		setMethod(
67		f = "[",
68		signature = c(x = "CompositionMatrix", i = "missing", j = "missing", drop = "logical"),
69		definition = function(x, i, j, ..., drop) {
70	!	na <- nargs()
71	!	if (na < 4L) return(x) # x[drop=], x[, drop=], x[drop=, ]
72	!	if (na == 4L) {
73	!	x <- if (drop) x@.Data else x
74	!	return(x) # x[, , drop=], x[, drop=, ], x[drop=, , ]
75		}
76	!	wrong_dimensions(".", ".") # x[, , , drop=], etc.
77		}
78		)
79
80		#' @export
81		#' @rdname subset
82		#' @aliases [,CompositionMatrix,index,missing,missing-method
83		setMethod(
84		f = "[",
85		signature = c(x = "CompositionMatrix", i = "index", j = "missing", drop = "missing"),
86		definition = function(x, i, j, ..., drop) {
87	4x	na <- nargs()
88	4x	if (na == 2L) {
89	2x	x <- .subscript1(x, i)
90	2x	return(x) # x[i=]
91		}
92	2x	if (na == 3L) {
93		#/!\ DROP /!\
94	1x	x <- .subscript2(x, i, , drop = FALSE)
95	1x	return(x) # x[i=, ], x[, i=]
96		}
97	1x	wrong_dimensions("i", ".") # x[i=, , ], etc.
98		}
99		)
100
101		#' @export
102		#' @rdname subset
103		#' @aliases [,CompositionMatrix,index,missing,logical-method
104		setMethod(
105		f = "[",
106		signature = c(x = "CompositionMatrix", i = "index", j = "missing", drop = "logical"),
107		definition = function(x, i, j, ..., drop) {
108	659x	na <- nargs()
109	659x	if (na == 3L) {
110	1x	x <- .subscript1(x, i)
111	1x	return(x) # x[i=, drop=]
112		}
113	658x	if (na == 4L) {
114	657x	x <- .subscript2(x, i, , drop = drop)
115	657x	return(x) # x[i=, , drop=], x[, i=, drop=]
116		}
117	1x	wrong_dimensions("i", ".") # x[i=, , , drop=], etc.
118		}
119		)
120
121		#' @export
122		#' @rdname subset
123		#' @aliases [,CompositionMatrix,missing,index,missing-method
124		setMethod(
125		f = "[",
126		signature = c(x = "CompositionMatrix", i = "missing", j = "index", drop = "missing"),
127		definition = function(x, i, j, ..., drop) {
128	23x	na <- nargs()
129	23x	if (na == 2L) {
130	1x	x <- .subscript1(x, j) # x[j=]
131	1x	return(x)
132		}
133	22x	if (na == 3L) {
134		# /!\ DROP /!\
135	21x	x <- .subscript2(x, , j, drop = FALSE) # x[j=, ], x[, j=]
136	21x	return(x)
137		}
138	1x	wrong_dimensions(".", "j") # x[, j=, ], etc.
139		}
140		)
141
142		#' @export
143		#' @rdname subset
144		#' @aliases [,CompositionMatrix,missing,index,logical-method
145		setMethod(
146		f = "[",
147		signature = c(x = "CompositionMatrix", i = "missing", j = "index", drop = "logical"),
148		definition = function(x, i, j, ..., drop) {
149	184x	na <- nargs()
150	184x	if (na == 3L) {
151	1x	x <- .subscript1(x, j) # x[j=, drop=]
152	1x	return(x)
153		}
154	183x	if (na == 4L) {
155	182x	x <- .subscript2(x, , j, drop = drop) # x[j=, , drop=], x[, j=, drop=]
156	182x	return(x)
157		}
158	1x	wrong_dimensions(".", "j") # x[, j=, , drop=], etc.
159		}
160		)
161
162		#' @export
163		#' @rdname subset
164		#' @aliases [,CompositionMatrix,index,index,missing-method
165		setMethod(
166		f = "[",
167		signature = c(x = "CompositionMatrix", i = "index", j = "index", drop = "missing"),
168		definition = function(x, i, j, ..., drop) {
169	2x	na <- nargs()
170	2x	if (na == 3L) {
171		# /!\ DROP /!\
172	1x	x <- .subscript2(x, i, j, drop = FALSE) # x[i=, j=], x[j=, i=]
173	1x	return(x)
174		}
175	1x	wrong_dimensions("i", "j") # x[i=, j=, ], etc.
176		}
177		)
178
179		#' @export
180		#' @rdname subset
181		#' @aliases [,CompositionMatrix,index,index,logical-method
182		setMethod(
183		f = "[",
184		signature = c(x = "CompositionMatrix", i = "index", j = "index", drop = "logical"),
185		definition = function(x, i, j, ..., drop) {
186	6x	na <- nargs()
187	6x	if (na == 4L) {
188	5x	x <- .subscript2(x, i, j, drop = drop) # x[i=, j=, drop=], x[j=, i=, drop=]
189	5x	return(x)
190		}
191	1x	wrong_dimensions("i", "j") # x[i=, j=, , drop=], etc.
192		}
193		)
194
195		# Replace ======================================================================
196		## [<- -------------------------------------------------------------------------
197		#' @export
198		#' @rdname subset
199		#' @aliases [<-,CompositionMatrix-method
200		setMethod(
201		f = "[<-",
202		signature = c(x = "CompositionMatrix"),
203		function(x, i, j, ..., value) {
204	5x	z <- methods::callNextMethod()
205	5x	methods::validObject(z)
206	5x	z
207		}
208		)
209
210		## [[<- ------------------------------------------------------------------------
211		#' @export
212		#' @rdname subset
213		#' @aliases [[<-,CompositionMatrix-method
214		setMethod(
215		f = "[[<-",
216		signature = c(x = "CompositionMatrix"),
217		function(x, i, j, ..., value) {
218	2x	z <- methods::callNextMethod()
219	2x	methods::validObject(z)
220	2x	z
221		}
222		)
223
224		# Transpose ====================================================================
225		#' @export
226		#' @rdname t
227		#' @aliases t,CompositionMatrix-method
228		setMethod(
229		f = "t",
230		signature = c(x = "CompositionMatrix"),
231		function(x) {
232	1x	t(x@.Data)
233		}
234		)
235
236		#' @export
237		#' @rdname t
238		#' @aliases t,LogRatio-method
239		setMethod(
240		f = "t",
241		signature = c(x = "LogRatio"),
242		function(x) {
243	1x	t(x@.Data)
244		}
245		)

1		# VARIATION MATRIX
2		#' @include AllGenerics.R
3		NULL
4
5		# Variation matrix =============================================================
6		#' @export
7		#' @rdname variation
8		#' @aliases variation,CompositionMatrix-method
9		setMethod(
10		f = "variation",
11		signature = c(x = "CompositionMatrix"),
12		definition = function(x) {
13	2x	J <- ncol(x)
14	2x	parts <- colnames(x)
15
16	2x	varia <- utils::combn(
17	2x	x = seq_len(J),
18	2x	m = 2,
19	2x	FUN = function(i, coda) {
20	20x	r <- coda[, i[[1]], drop = TRUE] / coda[, i[[2]], drop = TRUE]
21	20x	z <- log(r, base = exp(1))
22	20x	stats::var(z)
23		},
24	2x	coda = x
25		)
26
27	2x	mtx <- matrix(data = 0, nrow = J, ncol = J)
28	2x	mtx[lower.tri(mtx, diag = FALSE)] <- varia
29	2x	mtx <- t(mtx)
30	2x	mtx[lower.tri(mtx, diag = FALSE)] <- varia
31
32	2x	dimnames(mtx) <- list(parts, parts)
33	2x	mtx
34		}
35		)
36
37		#' @export
38		#' @rdname pip
39		#' @aliases pip,CompositionMatrix-method
40		setMethod(
41		f = "pip",
42		signature = c(x = "CompositionMatrix"),
43		definition = function(x) {
44	1x	v <- variation(x)
45	1x	1 / (1 + sqrt(v))
46		}
47		)
48
49		# Variation array ==============================================================
50		# @export
51		# @rdname variation_array
52		# @aliases variation_array,CompositionMatrix-method
53		# setMethod(
54		# f = "variation_array",
55		# signature = c(object = "CompositionMatrix"),
56		# definition = function(object) {
57		# J <- ncol(object)
58		# cbn <- utils::combn(seq_len(J), 2)
59		# varia <- apply(
60		# X = cbn,
61		# MARGIN = 2,
62		# FUN = function(j, x) {
63		# mean(log(x[, j[1]] / x[, j[2]]))
64		# },
65		# x = object
66		# )
67		#
68		# mtx <- variation(object)
69		# mtx[lower.tri(mtx, diag = FALSE)] <- varia
70		# mtx
71		# }
72		# )

1		# DATA TRANSFORMATION: INVERSE
2		#' @include AllGenerics.R
3		NULL
4
5		# Backtransform ================================================================
6		## CLR -------------------------------------------------------------------------
7		#' @export
8		#' @rdname transform_inverse
9		#' @aliases transform_inverse,CLR,missing-method
10		setMethod(
11		f = "transform_inverse",
12		signature = c(object = "CLR", origin = "missing"),
13		definition = function(object) {
14	3x	y <- methods::as(object, "matrix") # Drop slots
15	3x	y <- exp(y)
16	3x	y <- y / rowSums(y)
17
18	3x	dimnames(y) <- list(rownames(object), object@parts)
19	3x	.CompositionMatrix(y, totals = totals(object))
20		}
21		)
22
23		#' @export
24		#' @rdname transform_inverse
25		#' @aliases transform_inverse,GroupedCLR,missing-method
26		setMethod(
27		f = "transform_inverse",
28		signature = c(object = "GroupedCLR", origin = "missing"),
29		definition = function(object) {
30	!	z <- methods::callNextMethod()
31	!	.GroupedComposition(z, group_indices = group_indices(object),
32	!	group_levels = group_levels(object))
33		}
34		)
35
36		## ALR -------------------------------------------------------------------------
37		#' @export
38		#' @rdname transform_inverse
39		#' @aliases transform_inverse,ALR,missing-method
40		setMethod(
41		f = "transform_inverse",
42		signature = c(object = "ALR", origin = "missing"),
43		definition = function(object) {
44	3x	y <- exp(object)
45	3x	y <- y / (1 + rowSums(y))
46	3x	z <- 1 - rowSums(y)
47
48	3x	y <- cbind(y, z)
49	3x	dimnames(y) <- list(rownames(object), object@parts)
50	3x	y <- y[, object@order]
51
52	3x	.CompositionMatrix(y, totals = totals(object))
53		}
54		)
55
56		#' @export
57		#' @rdname transform_inverse
58		#' @aliases transform_inverse,GroupedALR,missing-method
59		setMethod(
60		f = "transform_inverse",
61		signature = c(object = "GroupedALR", origin = "missing"),
62		definition = function(object) {
63	!	z <- methods::callNextMethod()
64	!	.GroupedComposition(z, group_indices = group_indices(object),
65	!	group_levels = group_levels(object))
66		}
67		)
68
69		## ILR -------------------------------------------------------------------------
70		#' @export
71		#' @rdname transform_inverse
72		#' @aliases transform_inverse,ILR,missing-method
73		setMethod(
74		f = "transform_inverse",
75		signature = c(object = "ILR", origin = "missing"),
76		definition = function(object) {
77	2x	y <- tcrossprod(object@.Data, object@base)
78	2x	y <- exp(y)
79	2x	y <- y / rowSums(y)
80
81	2x	dimnames(y) <- list(rownames(object), object@parts)
82	2x	y <- y[, object@order]
83
84	2x	.CompositionMatrix(y, totals = totals(object))
85		}
86		)
87
88		#' @export
89		#' @rdname transform_inverse
90		#' @aliases transform_inverse,GroupedILR,missing-method
91		setMethod(
92		f = "transform_inverse",
93		signature = c(object = "GroupedILR", origin = "missing"),
94		definition = function(object) {
95	!	z <- methods::callNextMethod()
96	!	.GroupedComposition(z, group_indices = group_indices(object),
97	!	group_levels = group_levels(object))
98		}
99		)
100
101		#' @export
102		#' @rdname transform_inverse
103		#' @aliases transform_inverse,matrix,ILR-method
104		setMethod(
105		f = "transform_inverse",
106		signature = c(object = "matrix", origin = "ILR"),
107		definition = function(object, origin) {
108	1x	y <- tcrossprod(object, origin@base)
109	1x	y <- exp(y)
110	1x	y <- y / rowSums(y)
111
112	1x	dimnames(y) <- list(rownames(object), origin@parts)
113	1x	y <- y[, origin@order]
114
115	1x	y
116		}
117		)

1		# MUTATORS
2		#' @include AllGenerics.R
3		NULL
4
5		# Not exported
6		get_transformation <- function(x) {
7	!	if (methods::is(x, "LR")) return("Pairwise Log-Ratio")
8	!	if (methods::is(x, "CLR")) return("Centered Log-Ratio")
9	!	if (methods::is(x, "ALR")) return("Additive Log-Ratio")
10	2x	if (methods::is(x, "ILR")) return("Isometric Log-Ratio")
11	!	if (methods::is(x, "PLR")) return("Pivot Log-Ratio")
12		}
13
14		# Predicates ===================================================================
15		#' @export
16		#' @rdname CompositionMatrix-class
17		is_composition <- function(object) {
18	4x	methods::is(object, "CompositionMatrix")
19		}
20
21		#' @export
22		#' @rdname LogRatio-class
23		is_logratio <- function(object) {
24	!	methods::is(object, "LogRatio")
25		}
26
27		#' @export
28		#' @rdname ReferenceGroups-class
29		is_grouped <- function(object) {
30	137x	methods::is(object, "ReferenceGroups")
31		}
32
33		# Getter =======================================================================
34		#' @export
35		#' @method labels CompositionMatrix
36		labels.CompositionMatrix <- function(object, ...) {
37	3x	colnames(object)
38		}
39
40		#' @export
41		#' @rdname mutators
42		#' @aliases labels,CompositionMatrix-method
43		setMethod("labels", "CompositionMatrix", labels.CompositionMatrix)
44
45		#' @export
46		#' @method labels LogRatio
47		labels.LogRatio <- function(object, ...) {
48	2x	object@ratio
49		}
50
51		#' @export
52		#' @rdname mutators
53		#' @aliases labels,LogRatio-method
54		setMethod("labels", "LogRatio", labels.LogRatio)
55
56		# Weights ======================================================================
57		#' @export
58		#' @method weights ALR
59		weights.ALR <- function(object, ...) {
60	!	w <- object@weights
61	!	w[-1] * w[1]
62		}
63
64		#' @export
65		#' @rdname mutators
66		#' @aliases weights,ALR-method
67		setMethod("weights", "ALR", weights.ALR)
68
69		#' @export
70		#' @method weights LR
71		weights.LR <- function(object, ...) {
72	2x	w <- object@weights
73	2x	w <- utils::combn(
74	2x	x = w,
75	2x	m = 2,
76	2x	FUN = function(x) Reduce(`*`, x),
77	2x	simplify = FALSE
78		)
79	2x	unlist(w)
80		}
81
82		#' @export
83		#' @rdname mutators
84		#' @aliases weights,LR-method
85		setMethod("weights", "LR", weights.LR)
86
87		#' @export
88		#' @method weights LogRatio
89		weights.LogRatio <- function(object, ...) {
90	5x	object@weights
91		}
92
93		#' @export
94		#' @rdname mutators
95		#' @aliases weights,LogRatio-method
96		setMethod("weights", "LogRatio", weights.LogRatio)
97
98		# Totals =======================================================================
99		#' @export
100		#' @rdname totals
101		#' @aliases totals,CompositionMatrix-method
102	937x	setMethod("totals", "CompositionMatrix", function(object) object@totals)
103
104		#' @export
105		#' @rdname totals
106		#' @aliases totals,LogRatio-method
107	15x	setMethod("totals", "LogRatio", function(object) object@totals)
108
109		#' @export
110		#' @rdname totals
111		setMethod(
112		f = "totals<-",
113		signature = "CompositionMatrix",
114		definition = function(object, value) {
115	3x	object@totals <- if (is.null(value)) rowSums(object) else as.numeric(value)
116	3x	methods::validObject(object)
117	1x	object
118		}
119		)

1		# DATA SUMMARY: DESCRIBE
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname describe
7		#' @aliases describe,CompositionMatrix-method
8		setMethod(
9		f = "describe",
10		signature = c(x = "CompositionMatrix"),
11		definition = function(x) {
12		## Dimensions
13	!	m <- nrow(x)
14	!	p <- ncol(x)
15	!	rows <- sprintf(ngettext(m, "%d composition", "%d compositions"), m)
16	!	cols <- sprintf(ngettext(p, "with %d part", "with %d parts"), p)
17	!	msg_tbl <- sprintf("%s %s:", rows, cols)
18
19		## Message
20	!	cat(msg_tbl)
21	!	.describe(x)
22
23	!	invisible(x)
24		}
25		)
26
27		#' @export
28		#' @rdname describe
29		#' @aliases describe,GroupedComposition-method
30		setMethod(
31		f = "describe",
32		signature = c(x = "GroupedComposition"),
33		definition = function(x) {
34		## Dimensions
35	!	m <- nrow(x)
36	!	p <- ncol(x)
37
38	!	rows <- sprintf(ngettext(m, "%d composition", "%d compositions"), m)
39	!	cols <- sprintf(ngettext(p, "with %d part", "with %d parts"), p)
40	!	msg_tbl <- sprintf("%s %s:", rows, cols)
41
42		## Groups
43	!	i <- group_length(x)
44	!	ls_grp <- paste0(dQuote(group_levels(x)), collapse = ", ")
45	!	msg_grp <- sprintf(ngettext(i, "%d group", "%d groups"), i)
46	!	msg_grp <- sprintf("%s: %s.", msg_grp, ls_grp)
47
48	!	j <- sum(!is_assigned(x))
49	!	msg_ung <- sprintf(ngettext(j, "%d unassigned sample.", "%d unassigned samples."), j)
50
51	!	cat(msg_tbl, msg_grp, msg_ung, sep = "\n* ")
52	!	.describe(x)
53
54	!	invisible(x)
55		}
56		)
57
58		.describe <- function(x) {
59	!	m <- nrow(x)
60	!	p <- ncol(x)
61
62		## Missing values
63	!	n_NA <- sum(count(x, f = is.na))
64	!	m_NA <- sum(detect(x, f = is.na, margin = 1))
65	!	p_NA <- sum(detect(x, f = is.na, margin = 2))
66	!	pc <- label_percent(c(m_NA / m, p_NA / p), digits = 1, trim = TRUE)
67
68	!	msg_NA <- sprintf(ngettext(n_NA, "\n%d missing value:", "\n%d missing values:"), n_NA)
69
70	!	rows_NA <- ngettext(m_NA, "%d observation (%s) contains missing values.",
71	!	"%d observations (%s) contain missing values.")
72	!	msg_row_NA <- sprintf(rows_NA, m_NA, pc[[1]])
73
74	!	cols_NA <- ngettext(p_NA, "%d variable (%s) contains missing values.",
75	!	"%d variables (%s) contain missing values.")
76	!	msg_col_NA <- sprintf(cols_NA, p_NA, pc[[2]])
77
78		## Constant columns
79	!	p_var <- sum(detect(x, f = function(x) is_unique(x), margin = 2))
80	!	cols_var <- ngettext(p_var, "%d part with no variance.",
81	!	"%d parts with no variance.")
82	!	msg_col_var <- sprintf(cols_var, p_var)
83
84		## Sparsity
85	!	spa <- sparsity(x, count = FALSE)
86	!	msg_spa <- sprintf(tr_("%s of values are zero."), label_percent(spa, digits = 1))
87
88	!	cat(msg_NA, msg_row_NA, msg_col_NA, sep = "\n* ")
89	!	cat(tr_("\nData checking:"), msg_spa, msg_col_var, sep = "\n* ")
90		}

1		# DISTANCES
2		#' @include AllGenerics.R
3		NULL
4
5		# Distances ====================================================================
6		#' @export
7		#' @method dist CompositionMatrix
8		dist.CompositionMatrix <- function(x, method = "euclidean",
9		diag = FALSE, upper = FALSE, p = 2) {
10	2x	x <- transform_clr(x)
11	2x	stats::dist(x, method = method, diag = diag, upper = upper, p = p)
12		}
13
14		#' @export
15		#' @rdname dist
16		#' @aliases dist,CompositionMatrix-method
17		setMethod("dist", "CompositionMatrix", dist.CompositionMatrix)
18
19		# Mahalanobis ==================================================================
20		#' @export
21		#' @method mahalanobis CompositionMatrix
22		mahalanobis.CompositionMatrix <- function(x, center, cov, ..., robust = TRUE,
23		method = c("mve", "mcd")) {
24		## Transformation
25	1x	x <- transform_ilr(x)
26	1x	mahalanobis(x, center, cov, ..., robust = robust, method = method)
27		}
28
29		#' @export
30		#' @rdname mahalanobis
31		#' @aliases mahalanobis,CompositionMatrix-method
32		setMethod("mahalanobis", "CompositionMatrix", mahalanobis.CompositionMatrix)
33
34		#' @export
35		#' @method mahalanobis ILR
36		mahalanobis.ILR <- function(x, center, cov, ..., robust = TRUE,
37		method = c("mve", "mcd")) {
38
39	1x	if (missingORnull(center) \| missingORnull(cov)) {
40	1x	if (!robust) method <- "classical" # Standard estimators
41	!	else method <- match.arg(method, several.ok = FALSE) # Robust estimators
42	1x	v <- MASS::cov.rob(x, method = method, ...)
43		}
44
45	1x	est <- list(center = NULL, cov = NULL)
46	1x	est$center <- if (missingORnull(center)) v$center else center
47	1x	est$cov <- if (missingORnull(cov)) v$cov else cov
48
49	1x	message(v$sing)
50
51	1x	stats::mahalanobis(x, center = est$center, cov = est$cov)
52		}
53
54		#' @export
55		#' @rdname mahalanobis
56		#' @aliases mahalanobis,ILR-method
57		setMethod("mahalanobis", "ILR", mahalanobis.ILR)

1		# BARPLOT
2		#' @include AllGenerics.R
3		NULL
4
5		# CompositionMatrix ============================================================
6		#' @export
7		#' @method barplot CompositionMatrix
8		barplot.CompositionMatrix <- function(height, ...,
9		order_columns = FALSE, order_rows = NULL,
10		decreasing = TRUE, names = TRUE,
11		space = 0.2, offset = 0.025,
12		palette_color = palette_color_discrete(),
13		border = NA, axes = TRUE, legend = TRUE) {
14		## Validation
15	4x	if (ncol(height) < 2) {
16	!	stop(tr_("At least two compositional parts are needed."), call. = FALSE)
17		}
18
19		## Prepare data
20	4x	xy <- prepare_barplot(height, order_columns = order_columns,
21	4x	order_rows = order_rows, decreasing = decreasing,
22	4x	offset = offset)
23	4x	parts <- factor(xy$data$column, levels = colnames(height))
24	4x	col <- palette_color(parts)
25	4x	n <- nrow(height)
26
27		## Graphical parameters
28	4x	cex.axis <- list(...)$cex.axis %\|\|% graphics::par("cex.axis")
29	4x	col.axis <- list(...)$col.axis %\|\|% graphics::par("col.axis")
30	4x	font.axis <- list(...)$font.axis %\|\|% graphics::par("font.axis")
31
32		## Save and restore
33	4x	mfrow <- graphics::par("mfrow")
34	4x	mar <- graphics::par("mar")
35	4x	nlines <- height2line("M", cex = cex.axis)
36	4x	mar[1] <- 3
37	4x	mar[2] <- if (names) width2line(rownames(height), cex = cex.axis) + 0.1 else nlines
38	4x	mar[3] <- nlines
39	4x	mar[4] <- nlines
40
41	4x	old_par <- graphics::par(mfrow = mfrow, mar = mar)
42	4x	on.exit(graphics::par(old_par))
43
44		## Open new window
45	4x	grDevices::dev.hold()
46	4x	on.exit(grDevices::dev.flush(), add = TRUE)
47	4x	graphics::plot.new()
48
49		## Set plotting coordinates
50	4x	space <- 1 / n * 0.5 - space / n * 0.5
51	4x	xlim <- range(0, 1)
52	4x	ylim <- range(xy$data$y)
53	4x	ylim <- ylim + c(0, 2 * (offset + space) * legend) + c(-1, 1) * space
54	4x	graphics::plot.window(xlim = xlim, ylim = ylim, xaxs = "i", yaxs = "i")
55
56		## Plot
57	4x	graphics::rect(
58	4x	xleft = xy$data$xmin,
59	4x	ybottom = xy$data$y - space,
60	4x	xright = xy$data$xmax,
61	4x	ytop = xy$data$y + space,
62	4x	col = col,
63	4x	border = border
64		)
65
66		## Construct axis
67	4x	if (axes) {
68	4x	at <- graphics::axTicks(side = 1)
69	4x	graphics::axis(side = 1, at = at, labels = label_percent(at),
70	4x	xpd = NA, las = 1, cex.axis = cex.axis, col.axis = col.axis,
71	4x	font.axis = font.axis)
72	4x	graphics::mtext(text = names(xy$groups), side = 4, line = 0, at = xy$groups,
73	4x	cex = cex.axis, col = col.axis, font = font.axis)
74	4x	if (names) {
75	4x	graphics::axis(side = 2, at = unique(xy$data$y), labels = unique(xy$data$row),
76	4x	las = 2, lty = 0, cex.axis = cex.axis, col.axis = col.axis,
77	4x	font.axis = font.axis)
78		}
79		}
80
81		## Add legend
82	4x	if (legend) {
83	4x	graphics::rect(
84	4x	xleft = cumsum(xy$mean) - xy$mean,
85	4x	ybottom = max(xy$data$y) + offset + space,
86	4x	xright = cumsum(xy$mean),
87	4x	ytop = max(xy$data$y) + 2 * (offset + space),
88	4x	col = unique(col),
89	4x	border = border
90		)
91	4x	lab <- label_chemical(names(xy$mean))
92	4x	graphics::mtext(text = parse(text = lab), side = 3, line = 0,
93	4x	at = cumsum(xy$mean) - xy$mean / 2,
94	4x	cex = cex.axis, col = unique(col), font = font.axis)
95		}
96
97	4x	invisible(height)
98		}
99
100		#' @export
101		#' @rdname barplot
102		#' @aliases barplot,CompositionMatrix-method
103		setMethod("barplot", c(height = "CompositionMatrix"), barplot.CompositionMatrix)
104
105		prepare_barplot <- function(x, order_rows = NULL, order_columns = FALSE,
106		decreasing = TRUE, offset = 0.025) {
107		## Relative frequencies
108	4x	n <- nrow(x)
109	4x	x <- x / rowSums(x)
110
111		## Validation
112	4x	stopifnot(is_composition(x))
113
114		## Row order
115	4x	if (!is.null(order_rows)) {
116	1x	j <- x[, order_rows, drop = TRUE]
117	1x	i <- order(j, decreasing = decreasing)
118	1x	x <- x[i, , drop = FALSE]
119		}
120
121		## Columns order
122	4x	center <- mean(x)
123	4x	if (order_columns) {
124	2x	col_order <- order(center, decreasing = FALSE)
125	2x	center <- center[col_order]
126	2x	x <- x[, col_order, drop = FALSE]
127		}
128
129		## Grouping
130	4x	if (!is_grouped(x)) {
131	3x	x <- group(x, by = rep(NA, n), verbose = FALSE)
132		}
133	4x	spl <- group_split(x)
134	4x	z <- do.call(rbind, spl)
135
136		## Build a long table
137	4x	row <- row(z, as.factor = TRUE)
138	4x	col <- col(z, as.factor = TRUE)
139	4x	data <- data.frame(
140	4x	row = as.vector(row),
141	4x	column = as.vector(col),
142	4x	value = as.vector(z)
143		)
144
145	4x	xmax <- t(apply(X = z, MARGIN = 1, FUN = cumsum))
146	4x	xmin <- xmax - z
147	4x	data$xmin <- as.vector(xmin)
148	4x	data$xmax <- as.vector(xmax)
149	4x	data$y <- as.vector(n + 1 - as.numeric(row)) / n # Reverse levels order
150
151		## Offset
152	4x	n_grp <- group_length(x)
153	4x	n_spl <- group_size(x)
154	4x	offset <- rev(seq_len(n_grp)) * offset - offset
155	4x	data$y <- data$y + rep(offset, n_spl)[as.numeric(row)]
156
157	4x	list(
158	4x	data = data,
159	4x	mean = center,
160	4x	groups = 1 - cumsum(n_spl) / n + n_spl / n * 0.5 + offset
161		)
162		}

1		# OUTLIERS
2		#' @include AllGenerics.R
3		NULL
4
5		# Find =========================================================================
6		#' @export
7		#' @rdname detect_outlier
8		#' @aliases detect_outlier,CompositionMatrix,missing-method
9		setMethod(
10		f = "detect_outlier",
11		signature = c(object = "CompositionMatrix", reference = "missing"),
12		definition = function(object, ..., robust = TRUE, method = c("mve", "mcd"),
13		quantile = 0.975) {
14	1x	methods::callGeneric(object, reference = object, robust = robust,
15	1x	method = method, quantile = quantile)
16		}
17		)
18
19		#' @export
20		#' @rdname detect_outlier
21		#' @aliases detect_outlier,CompositionMatrix,CompositionMatrix-method
22		setMethod(
23		f = "detect_outlier",
24		signature = c(object = "CompositionMatrix", reference = "CompositionMatrix"),
25		definition = function(object, reference, ..., quantile = 0.975, robust = TRUE,
26		method = c("mve", "mcd")) {
27		## Validation
28	1x	if (!identical(colnames(object), colnames(reference))) {
29	!	stop(tr_("Column names do not match!"), call. = FALSE)
30		}
31
32		## Transformation
33	1x	obj <- transform_ilr(object)
34	1x	ref <- transform_ilr(reference)
35
36		## Clean
37	1x	n <- nrow(ref)
38	1x	p <- ncol(ref)
39	1x	if (n < (p + 1)) {
40	!	msg <- tr_("Sample size is too small (%d).")
41	!	stop(sprintf(msg, n), call. = FALSE)
42		}
43	1x	if (n < (2 * p)) {
44	!	msg <- tr_("Possibly too small sample size (%d).")
45	!	warning(sprintf(msg, n), call. = FALSE)
46		}
47
48		## Compute center and spread + Mahalanobis distance
49		## Standard estimators
50	1x	estc <- list(center = colMeans(ref, na.rm = TRUE), cov = cov(ref))
51	1x	dc <- stats::mahalanobis(obj, center = estc$center, cov = estc$cov)
52
53		## Robust estimators
54	1x	dr <- rep(NA_real_, nrow(obj))
55	1x	if (robust) {
56	!	method <- match.arg(method, several.ok = FALSE)
57	!	estr <- MASS::cov.rob(ref, method = method, ...)
58	!	dr <- stats::mahalanobis(obj, center = estr$center, cov = estr$cov)
59		}
60
61		## Threshold
62	1x	limit <- stats::qchisq(p = quantile, df = p)
63
64	1x	.OutlierIndex(
65	1x	samples = rownames(obj),
66	1x	standard = sqrt(dc),
67	1x	robust = sqrt(dr),
68	1x	limit = sqrt(limit),
69	1x	dof = p
70		)
71		}
72		)
73
74		#' @export
75		#' @rdname detect_outlier
76		#' @aliases is_outlier,OutlierIndex-method
77		setMethod(
78		f = "is_outlier",
79		signature = c("OutlierIndex"),
80		definition = function(object, robust = TRUE) {
81	1x	d <- if (robust) object@robust else object@standard
82	1x	out <- d > object@limit
83	1x	names(out) <- object@samples
84	1x	out
85		}
86		)
87
88		# Plot =========================================================================
89		#' @export
90		#' @method plot OutlierIndex
91		plot.OutlierIndex <- function(x, ...,
92		type = c("dotchart", "distance"),
93		robust = TRUE,
94		symbols = c(16, 1, 3),
95		xlim = NULL, ylim = NULL,
96		xlab = NULL, ylab = NULL,
97		main = NULL, sub = NULL,
98		ann = graphics::par("ann"),
99		axes = TRUE, frame.plot = axes,
100		panel.first = NULL, panel.last = NULL,
101		legend = list(x = "topleft")) {
102		## Get data
103	1x	dc <- x@standard
104	1x	dr <- x@robust
105	1x	dof <- x@dof
106	1x	limit <- x@limit
107	1x	n <- length(dc)
108
109		## Validation
110	1x	if (all(is.na(dr))) {
111	1x	robust <- FALSE
112	1x	type <- "dotchart"
113		}
114	1x	type <- match.arg(type, several.ok = FALSE)
115
116		## Graphical parameters
117	1x	shape <- rep(symbols[[1L]], n)
118	!	if (robust \|\| type == "distance") shape[dr > limit] <- symbols[[3L]]
119	1x	if (!robust \|\| type == "distance") shape[dc > limit] <- symbols[[2L]]
120
121	1x	cy <- if (robust) dr else dc
122	1x	dlab <- ifelse(robust, tr_("Robust Mahalanobis distance"),
123	1x	tr_("Standard Mahalanobis distance"))
124	1x	ylab <- ylab %\|\|% dlab
125
126	1x	if (type == "dotchart") {
127	1x	asp <- NA
128	1x	cx <- seq_along(dc)
129	1x	xlab <- xlab %\|\|% tr_("Index")
130	1x	panel <- function() {
131	1x	graphics::points(x = cx, y = cy, pch = shape, ...)
132	1x	graphics::abline(h = limit, lty = 1)
133		}
134		}
135	1x	if (type == "distance") {
136	!	asp <- 1
137	!	cx <- dc
138	!	cy <- dr
139	!	xlab <- xlab %\|\|% tr_("Standard Mahalanobis distance")
140	!	ylab <- ylab %\|\|% tr_("Robust Mahalanobis distance")
141	!	panel <- function() {
142	!	graphics::points(x = cx, y = cy, pch = shape, ...)
143	!	graphics::abline(h = limit, lty = 1)
144	!	graphics::abline(v = limit, lty = 1)
145	!	graphics::abline(a = 0, b = 1, lty = 2, col = "darkgrey")
146		}
147		}
148
149		## Open new window
150	1x	grDevices::dev.hold()
151	1x	on.exit(grDevices::dev.flush(), add = TRUE)
152	1x	graphics::plot.new()
153
154		## Set plotting coordinates
155	1x	xlim <- xlim %\|\|% range(cx, finite = TRUE)
156	1x	ylim <- ylim %\|\|% range(cy, finite = TRUE)
157	1x	graphics::plot.window(xlim = xlim, ylim = ylim, asp = asp)
158
159		## Evaluate pre-plot expressions
160	1x	panel.first
161
162		## Plot
163	1x	panel()
164
165		## Evaluate post-plot and pre-axis expressions
166	1x	panel.last
167
168		## Construct Axis
169	1x	if (axes) {
170	1x	graphics::axis(side = 1, las = 1)
171	1x	graphics::axis(side = 2, las = 1)
172		}
173
174		## Plot frame
175	1x	if (frame.plot) {
176	1x	graphics::box()
177		}
178
179		## Add annotation
180	1x	if (ann) {
181	1x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab)
182		}
183
184		## Add legend
185	1x	if (is.list(legend)) {
186	1x	if (type == "distance") {
187	!	lab <- c(tr_("No outlier"), tr_("Robust only"), tr_("Both"))
188	!	pch <- symbols
189		} else {
190	1x	lab <- c(tr_("No outlier"), tr_("Outlier"))
191	1x	pch <- symbols[-2 - !robust]
192		}
193	1x	args <- list(x = "topleft", legend = lab, pch = pch, bty = "n", xpd = NA)
194	1x	args <- utils::modifyList(args, legend)
195	1x	do.call(graphics::legend, args = args)
196		}
197
198	1x	invisible(x)
199		}
200
201		#' @export
202		#' @rdname plot_outlier
203		#' @aliases plot,OutlierIndex,missing-method
204		setMethod("plot", c(x = "OutlierIndex", y = "missing"), plot.OutlierIndex)

1		# BIND
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname bind
7		#' @aliases rbind2,CompositionMatrix,CompositionMatrix-method
8		setMethod(
9		f = "rbind2",
10		signature = c(x = "CompositionMatrix", y = "CompositionMatrix"),
11		definition = function(x, y) {
12	4x	mtx_x <- methods::as(x, "matrix")
13	4x	mtx_y <- methods::as(y, "matrix")
14
15	4x	spl <- c(rownames(x), rownames(y))
16	4x	if (any(duplicated(spl))) {
17	!	warning(tr_("Duplicated rownames!"), call. = FALSE)
18	!	spl <- make.unique(spl, sep = "_")
19		}
20
21	4x	z <- rbind(mtx_x, mtx_y)
22	4x	rownames(z) <- spl
23
24	4x	.CompositionMatrix(z, totals = c(totals(x), totals(y)))
25		}
26		)
27
28		#' @export
29		#' @rdname bind
30		#' @aliases rbind2,GroupedComposition,GroupedComposition-method
31		setMethod(
32		f = "rbind2",
33		signature = c(x = "GroupedComposition", y = "GroupedComposition"),
34		definition = function(x, y) {
35	!	stop(tr_("Combining grouped matrices by rows is not currently supported."),
36	!	call. = FALSE)
37		}
38		)

1		# REPLACE
2		#' @include AllGenerics.R
3		NULL
4
5		# Zeros ========================================================================
6		#' @export
7		#' @rdname replace_zero
8		#' @aliases replace_zero,CompositionMatrix-method
9		setMethod(
10		f = "replace_zero",
11		signature = c(x = "CompositionMatrix"),
12		definition = function(x, value, delta = 2/3) {
13		## Validation
14	1x	D <- ncol(x)
15	!	if (is.null(value)) return(x)
16	!	if (length(value) == 1) value <- rep(value, D)
17	1x	if (length(value) > 1) arkhe::assert_length(value, D)
18	!	if (length(delta) > 1) arkhe::assert_length(delta, D)
19
20	1x	sigma <- value * delta
21	1x	r <- apply(X = x, MARGIN = 1, FUN = zero_multiplicative, sigma = sigma)
22
23	1x	methods::initialize(x, t(r))
24		}
25		)
26
27		# zero_additive <- function(x, sigma) {
28		# D <- length(x)
29		#
30		# is_zero <- x == 0
31		# Z <- sum(is_zero)
32		#
33		# x[is_zero] <- (sigma * (Z + 1) * (D - Z)) / D^2
34		# x[!is_zero] <- x[!is_zero] - (sigma * (Z + 1) * Z) / D^2
35		#
36		# x
37		# }
38		zero_multiplicative <- function(x, sigma) {
39	9x	is_zero <- x == 0 & !is.na(x)
40	9x	x[is_zero] <- sigma[is_zero]
41	9x	x[!is_zero] <- x[!is_zero] * (1 - sum(sigma[is_zero]) / 1)
42	9x	x
43		}
44
45		# Missing values ===============================================================
46		#' @export
47		#' @rdname replace_NA
48		#' @aliases replace_NA,CompositionMatrix-method
49		setMethod(
50		f = "replace_NA",
51		signature = c(x = "CompositionMatrix"),
52		definition = function(x, value) {
53		## Validation
54	2x	D <- ncol(x)
55	!	if (is.null(value)) return(x)
56	2x	if (length(value) == 1) value <- rep(value, D)
57	2x	if (length(value) > 1) arkhe::assert_length(value, D)
58
59	2x	r <- apply(X = x, MARGIN = 1, FUN = missing_multiplicative, sigma = value)
60
61	2x	methods::initialize(x, t(r))
62		}
63		)
64
65		missing_multiplicative <- function(x, sigma) {
66	18x	is_missing <- is.na(x)
67	18x	x[is_missing] <- sigma[is_missing]
68	18x	x[!is_missing] <- x[!is_missing] * (1 - sum(sigma[is_missing])) / sum(x[!is_missing])
69	18x	x
70		}

1		# MIXED-MODE ANALYSIS
2		#' @include AllGenerics.R
3		NULL
4
5		# First approach ===============================================================
6		#' @export
7		#' @describeIn mix First approach of mixed-mode analysis.
8		#' @aliases mix,matrix,matrix-method
9		setMethod(
10		f = "mix",
11		signature = c(x = "matrix", y = "matrix"),
12		definition = function(x, y, lambda = 1, ...) {
13		## Validation
14	!	arkhe::assert_package("cluster")
15	!	stopifnot(nrow(y) == nrow(x))
16
17	!	X <- vector(mode = "list", length = lambda + 1)
18	!	X[[1]] <- y
19	!	lambda <- as.integer(lambda[[1L]])
20	!	for (i in seq_len(lambda)) X[[i + 1]] <- x
21	!	X <- do.call(cbind, X)
22
23	!	d <- cluster::daisy(X, metric = "gower", ...)
24	!	as.dist(d)
25		}
26		)
27
28		# Second approach ==============================================================
29		#' @export
30		#' @describeIn mix Second approach of mixed-mode analysis.
31		#' @aliases mix,dist,dist-method
32		setMethod(
33		f = "mix",
34		signature = c(x = "dist", y = "dist"),
35		definition = function(x, y, mu = 0.5) {
36		## Validation
37	!	stopifnot(mu >= 0 & mu <= 1)
38
39	!	x <- as.matrix(x)
40	!	y <- as.matrix(y)
41
42	!	d <- mu * x + (1 - mu) * y
43	!	as.dist(d)
44		}
45		)

1		# DATA TRANSFORMATION: ADDITIVE LOG RATIO
2		#' @include AllGenerics.R
3		NULL
4
5		# ALR ==========================================================================
6		alr_base <- function(D) {
7	13x	V <- diag(1, nrow = D, ncol = D - 1)
8	13x	V[D, ] <- -1
9
10	13x	V
11		}
12
13		#' @export
14		#' @rdname transform_alr
15		#' @aliases transform_alr,CompositionMatrix-method
16		setMethod(
17		f = "transform_alr",
18		signature = c(object = "CompositionMatrix"),
19		definition = function(object, j = ncol(object), weights = FALSE) {
20	11x	D <- ncol(object)
21	11x	parts <- colnames(object)
22
23		## Reorder (move denominator)
24	11x	j <- if (is.character(j)) which(parts == j) else as.integer(j)
25	11x	ordering <- c(which(j != seq_len(D)), j)
26	11x	parts <- parts[ordering]
27	11x	z <- object[, ordering, drop = FALSE]
28
29		## Compute ratios
30	11x	base <- alr_base(D)
31	11x	alr <- log(z, base = exp(1)) %*% base
32	11x	rownames(alr) <- rownames(object)
33	11x	colnames(alr) <- paste(parts[-D], parts[D], sep = "_")
34
35		## Compute weights
36	11x	weights <- make_weights(object, weights = weights)
37
38	11x	.ALR(
39	11x	alr,
40	11x	parts = parts,
41	11x	ratio = paste(parts[-D], parts[D], sep = "/"),
42	11x	order = order(ordering),
43	11x	base = base,
44	11x	weights = weights,
45	11x	totals = totals(object)
46		)
47		}
48		)
49
50		#' @export
51		#' @rdname transform_alr
52		#' @aliases transform_alr,GroupedComposition-method
53		setMethod(
54		f = "transform_alr",
55		signature = c(object = "GroupedComposition"),
56		definition = function(object, j = ncol(object), weights = FALSE) {
57	!	z <- methods::callNextMethod()
58	!	.GroupedALR(z, group_indices = group_indices(object),
59	!	group_levels = group_levels(object))
60		}
61		)
62
63		#' @export
64		#' @rdname transform_alr
65		#' @aliases transform_alr,CLR-method
66		setMethod(
67		f = "transform_alr",
68		signature = c(object = "CLR"),
69		definition = function(object, j = ncol(object)) {
70	2x	D <- ncol(object)
71	2x	parts <- object@parts
72
73		## Reorder (move denominator)
74	2x	j <- if (is.character(j)) which(parts == j) else as.integer(j)
75	2x	ordering <- c(which(j != seq_len(D)), j)
76	2x	parts <- parts[ordering]
77	2x	z <- object[, ordering, drop = FALSE]
78
79		## Compute ratios
80	2x	base <- alr_base(D)
81	2x	alr <- z %*% base
82	2x	rownames(alr) <- rownames(object)
83	2x	colnames(alr) <- paste(parts[-D], parts[D], sep = "_")
84
85	2x	.ALR(
86	2x	alr,
87	2x	parts = parts,
88	2x	ratio = paste(parts[-D], parts[D], sep = "/"),
89	2x	order = order(ordering),
90	2x	base = base,
91	2x	weights = object@weights,
92	2x	totals = totals(object)
93		)
94		}
95		)
96
97		#' @export
98		#' @rdname transform_alr
99		#' @aliases transform_alr,GroupedCLR-method
100		setMethod(
101		f = "transform_alr",
102		signature = c(object = "GroupedCLR"),
103		definition = function(object, j = ncol(object), weights = FALSE) {
104	!	z <- methods::callNextMethod()
105	!	.GroupedALR(z, group_indices = group_indices(object),
106	!	group_levels = group_levels(object))
107		}
108		)

1		# DATA TRANSFORMATION: CENTERED LOG RATIO
2		#' @include AllGenerics.R
3		NULL
4
5		# CLR ==========================================================================
6		clr_base <- function(D, weights = rep(1 / D, D)) {
7	23x	diag(D) - matrix(data = weights, nrow = D, ncol = D)
8		}
9
10		#' @export
11		#' @rdname transform_clr
12		#' @aliases transform_clr,CompositionMatrix-method
13		setMethod(
14		f = "transform_clr",
15		signature = c(object = "CompositionMatrix"),
16		definition = function(object, weights = FALSE) {
17	20x	D <- ncol(object)
18	20x	parts <- colnames(object)
19
20	20x	weights <- make_weights(object, weights = weights)
21	20x	base <- clr_base(D, weights = weights)
22	20x	clr <- log(object, base = exp(1)) %*% base
23	20x	dimnames(clr) <- dimnames(object)
24
25	20x	.CLR(
26	20x	clr,
27	20x	parts = parts,
28	20x	ratio = parts,
29	20x	order = seq_len(D),
30	20x	base = base,
31	20x	weights = weights,
32	20x	totals = totals(object)
33		)
34		}
35		)
36
37		#' @export
38		#' @rdname transform_clr
39		#' @aliases transform_clr,GroupedComposition-method
40		setMethod(
41		f = "transform_clr",
42		signature = c(object = "GroupedComposition"),
43		definition = function(object, weights = FALSE) {
44	!	z <- methods::callNextMethod()
45	!	.GroupedCLR(z, group_indices = group_indices(object),
46	!	group_levels = group_levels(object))
47		}
48		)
49
50		#' @export
51		#' @rdname transform_clr
52		#' @aliases transform_clr,ALR-method
53		setMethod(
54		f = "transform_clr",
55		signature = c(object = "ALR"),
56		definition = function(object) {
57	3x	D <- ncol(object) + 1
58	3x	w <- object@weights
59
60	3x	base <- clr_base(D, weights = w)
61	3x	clr <- object %*% base[-D, ]
62	3x	dimnames(clr) <- list(rownames(object), object@parts)
63
64	3x	.CLR(
65	3x	clr,
66	3x	parts = object@parts,
67	3x	ratio = object@parts,
68	3x	order = seq_len(D),
69	3x	base = base,
70	3x	weights = w,
71	3x	totals = totals(object)
72		)
73		}
74		)
75
76		#' @export
77		#' @rdname transform_clr
78		#' @aliases transform_clr,GroupedALR-method
79		setMethod(
80		f = "transform_clr",
81		signature = c(object = "GroupedALR"),
82		definition = function(object) {
83	!	z <- methods::callNextMethod()
84	!	.GroupedCLR(z, group_indices = group_indices(object),
85	!	group_levels = group_levels(object))
86		}
87		)

1		# DATA TRANSFORMATION: LOG RATIO
2		#' @include AllGenerics.R
3		NULL
4
5		# LR ===========================================================================
6		#' @export
7		#' @rdname transform_lr
8		#' @aliases transform_lr,CompositionMatrix-method
9		setMethod(
10		f = "transform_lr",
11		signature = c(object = "CompositionMatrix"),
12		definition = function(object, weights = FALSE) {
13	6x	J <- ncol(object)
14	6x	parts <- colnames(object)
15
16		## Compute weights
17	6x	weights <- make_weights(object, weights = weights)
18
19		## Computes ratios
20	6x	jj <- utils::combn(seq_len(J), 2, simplify = FALSE)
21	6x	lr <- matrix(data = 0, nrow = nrow(object), ncol = length(jj))
22	6x	for (i in seq_along(jj)) {
23	46x	a <- jj[[i]][[1]]
24	46x	b <- jj[[i]][[2]]
25	46x	r <- object[, a, drop = TRUE] / object[, b, drop = TRUE]
26	46x	lr[, i] <- log(r, base = exp(1))
27		}
28
29		## Make names
30	6x	ratio <- unlist(utils::combn(parts, 2, FUN = paste, collapse = "/", simplify = FALSE))
31	6x	rownames(lr) <- rownames(object)
32	6x	colnames(lr) <- ratio
33
34	6x	.LR(
35	6x	lr,
36	6x	parts = parts,
37	6x	ratio = ratio,
38	6x	order = seq_len(J),
39	6x	weights = weights,
40	6x	totals = totals(object)
41		)
42		}
43		)
44
45		#' @export
46		#' @rdname transform_lr
47		#' @aliases transform_lr,GroupedComposition-method
48		setMethod(
49		f = "transform_lr",
50		signature = c(object = "GroupedComposition"),
51		definition = function(object, weights = FALSE) {
52	!	z <- methods::callNextMethod()
53	!	.GroupedLR(z, group_indices = group_indices(object),
54	!	group_levels = group_levels(object))
55		}
56		)

1		# COMPOSITIONAL MEAN
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname condense
7		#' @aliases condense,CompositionMatrix-method
8		setMethod(
9		f = "condense",
10		signature = "CompositionMatrix",
11		definition = function(x, by, ignore_na = FALSE, ignore_zero = TRUE,
12		verbose = FALSE, ...) {
13	3x	x <- group(x, by = by)
14	3x	y <- methods::callGeneric(x = x, ignore_na = ignore_na,
15	3x	ignore_zero = ignore_zero, verbose = verbose, ...)
16	3x	ungroup(y)
17		}
18		)
19
20		#' @export
21		#' @rdname condense
22		#' @aliases condense,GroupedComposition-method
23		setMethod(
24		f = "condense",
25		signature = "GroupedComposition",
26		definition = function(x, by = NULL, ignore_na = FALSE, ignore_zero = TRUE,
27		verbose = FALSE, ...) {
28		## Grouping
29	5x	grp <- group_factor(x)
30	1x	if (!is.null(by)) x <- group(x, by = by, verbose = verbose)
31
32		## Compute mean
33	5x	z <- aggregate(
34	5x	x = x,
35	5x	FUN = mean,
36	5x	ignore_na = ignore_na,
37	5x	ignore_zero = ignore_zero,
38	5x	simplify = TRUE
39		)
40	5x	tot <- tapply(X = totals(x), INDEX = group_factor(x), FUN = mean)
41
42	5x	z <- .CompositionMatrix(z, totals = as.numeric(tot))
43	5x	group(z, by = flatten_chr(x = grp, by = group_factor(x)), verbose = verbose)
44		}
45		)
46
47		flatten_chr <- function(x, by) {
48	5x	x <- as.character(x)
49	5x	z <- tapply(X = x, INDEX = by, FUN = unique, simplify = FALSE)
50	5x	z <- vapply(X = z, FUN = paste0, FUN.VALUE = character(1), collapse = ":")
51	5x	z
52		}

1		# HELPERS
2
3		## https://michaelchirico.github.io/potools/articles/developers.html
4		tr_ <- function(...) {
5	7x	enc2utf8(gettext(paste0(...), domain = "R-nexus"))
6		}
7
8		missingORnull <- function(x) {
9	4x	missing(x) \|\| is.null(x)
10		}
11
12		has_rownames <- function(x) {
13	51x	.row_names_info(x, type = 1L) > 0L &&
14	51x	!is.na(.row_names_info(x, type = 0L)[[1L]])
15		}
16
17		make_names <- function(x, n = length(x), prefix = "X") {
18	102x	x <- if (n > 0) x %\|\|% paste0(prefix, seq_len(n)) else character(0)
19	102x	x <- make.unique(x, sep = "_")
20	102x	x
21		}
22
23		#' Plotting Dimensions of Character Strings
24		#'
25		#' Convert string length in inch to number of (margin) lines.
26		#' @param x A [`character`] vector of string whose length is to be calculated.
27		#' @param ... Further parameter to be passed to [graphics::strwidth()]`, such as
28		#' `cex`.
29		#' @return
30		#' A [`numeric`] vector (maximum string width in units of margin lines).
31		#' @note For internal use only.
32		#' @family graphic tools
33		#' @keywords internal
34		#' @noRd
35		width2line <- function(x, ...) {
36	4x	(max(graphics::strwidth(x, units = "inch", ...)) /
37	4x	graphics::par("cin")[2] + graphics::par("mgp")[2]) * graphics::par("cex")
38		}
39		height2line <- function(x, ...) {
40	4x	(max(graphics::strheight(x, units = "inch", ...)) /
41	4x	graphics::par("cin")[2] + graphics::par("mgp")[2]) * graphics::par("cex")
42		}
43
44		#' Label Percentages
45		#'
46		#' @param x A [`numeric`] vector.
47		#' @param digits An [`integer`] indicating the number of decimal places.
48		#' If `NULL` (the default), breaks will have the minimum number of digits
49		#' needed to show the difference between adjacent values.
50		#' @param trim A [`logical`] scalar. If `FALSE` (the default), values are
51		#' right-justified to a common width (see [base::format()]).
52		#' @return A [`character`] vector.
53		#' @keywords internal
54		#' @noRd
55		label_percent <- function(x, digits = NULL, trim = FALSE) {
56	6x	i <- !is.na(x)
57	6x	y <- x[i]
58	6x	y <- abs(y) * 100
59	6x	y <- format(y, trim = trim, digits = digits)
60	6x	y <- paste0(y, "%")
61	6x	x[i] <- y
62	6x	x
63		}
64
65		#' Label Chemical Formula
66		#'
67		#' @param x A [`character`] vector.
68		#' @return A [`character`] vector.
69		#' @keywords internal
70		#' @noRd
71		label_chemical <- function(x, digits = NULL, trim = FALSE) {
72	!	if (!all(is_chemical(x))) return(x)
73	5x	x <- gsub(pattern = "([[:digit:]]+)", replacement = "[\\1]", x = x, fixed = FALSE)
74	5x	x <- gsub(pattern = "([[:alpha:]]+)", replacement = "\"\\1\"", x = x, fixed = FALSE)
75	5x	gsub(pattern = "]\"", replacement = "]*\"", x = x, fixed = FALSE)
76		}
77
78		#' Column Weights
79		#'
80		#' Computes column weights.
81		#' @param x A `numeric` [`matrix`].
82		#' @param weights A [`logical`] scalar: should varying weights (column means)
83		#' be computed? If `FALSE` (the default), equally-weighted parts are used.
84		#' Alternatively, a positive [`numeric`] vector of weights can be specified
85		#' (will be rescaled to sum to \eqn{1}).
86		#' @return A [`numeric`] vector.
87		#' @keywords internal
88		#' @noRd
89		make_weights <- function(x, weights = FALSE) {
90	37x	D <- ncol(x)
91
92	37x	w <- if (isTRUE(weights)) colMeans(x) else rep(1 / D, D)
93	37x	if (is.numeric(weights)) {
94	3x	arkhe::assert_length(weights, D)
95	3x	arkhe::assert_positive(weights, strict = TRUE)
96	3x	w <- weights / sum(weights) # Sum up to 1
97		}
98
99	37x	unname(w)
100		}

1		# PCA
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @describeIn pca PCA of centered log-ratio, i.e. log-ratio analysis (LRA).
7		#' @aliases pca,CompositionMatrix-method
8		setMethod(
9		f = "pca",
10		signature = c("CompositionMatrix"),
11		definition = function(object, center = TRUE, scale = FALSE, rank = NULL,
12		sup_row = NULL, sup_col = NULL,
13		weight_row = NULL, weight_col = NULL) {
14	1x	message(tr_("PCA of centered log-ratio (CLR)."))
15	1x	x <- transform_clr(object)
16	1x	methods::callGeneric(object = x, center = center, scale = scale,
17	1x	rank = rank, sup_row = sup_row, sup_col = sup_col,
18	1x	weight_row = weight_row, weight_col = weight_col)
19		}
20		)
21
22		#' @export
23		#' @rdname pca
24		#' @aliases pca,LogRatio-method
25		setMethod(
26		f = "pca",
27		signature = c("LogRatio"),
28		definition = function(object, center = TRUE, scale = FALSE, rank = NULL,
29		sup_row = NULL, sup_col = NULL,
30		weight_row = NULL, weight_col = NULL) {
31	1x	methods::callNextMethod(object, center = center, scale = scale,
32	1x	rank = rank, sup_row = sup_row,
33	1x	sup_col = sup_col, weight_row = weight_row,
34	1x	weight_col = weight_col)
35		}
36		)

1		# MARGIN
2		#' @include AllGenerics.R
3		NULL
4
5		# Margin =======================================================================
6		#' @export
7		#' @rdname margin
8		#' @aliases margin,CompositionMatrix-method
9		setMethod(
10		f = "margin",
11		signature = c("CompositionMatrix"),
12		definition = function(x, parts = c(1, 2), name = "*") {
13		## Validation
14	1x	p <- NCOL(x)
15	1x	parts <- unique(parts)
16	1x	if (is.character(parts)) parts <- match(parts, colnames(x))
17	!	if (p <= length(parts)) return(x)
18	1x	if (p == length(parts) + 1) {
19	!	d <- seq_len(p)
20	!	d <- c(d[parts], d[-parts])
21	!	return(x[, d, drop = FALSE])
22		}
23
24	1x	rest <- x[, -parts, drop = FALSE]
25	1x	star <- apply(X = rest, MARGIN = 1, FUN = gmean)
26	1x	mar <- cbind(x[, parts, drop = FALSE], star)
27	1x	colnames(mar) <- c(colnames(x)[parts], name[[1L]])
28
29	1x	clo <- closure(mar)
30	1x	methods::initialize(x, clo)
31		}
32		)

1		# PLOT
2		#' @include AllGenerics.R
3		NULL
4
5		# LogRatio =====================================================================
6		#' @export
7		#' @method plot LogRatio
8		plot.LogRatio <- function(x, ..., factor = 1, amount = NULL,
9		palette_color = palette_color_discrete(),
10		palette_symbol = palette_shape(),
11		xlab = NULL, ylab = NULL,
12		main = NULL, sub = NULL, ann = graphics::par("ann"),
13		axes = TRUE, frame.plot = axes,
14		legend = list(x = "topright")) {
15		## Get data
16	2x	xy <- data.frame(
17	2x	x = jitter(as.vector(col(x)), factor = factor, amount = amount),
18	2x	y = as.vector(x)
19		)
20
21		## Graphical parameters
22	2x	col <- list(...)$col %\|\|% graphics::par("col")
23	2x	bg <- list(...)$bg %\|\|% graphics::par("bg")
24	2x	pch <- list(...)$pch %\|\|% graphics::par("pch")
25	2x	cex.axis <- list(...)$cex.axis %\|\|% graphics::par("cex.axis")
26	2x	col.axis <- list(...)$col.axis %\|\|% graphics::par("col.axis")
27	2x	font.axis <- list(...)$font.axis %\|\|% graphics::par("font.axis")
28
29		## Grouping
30	2x	if (is_grouped(x)) {
31	1x	lvl <- group_names(x)
32	1x	col <- palette_color(lvl)
33	1x	bg <- grDevices::adjustcolor(col, alpha.f = 0.5)
34	1x	pch <- palette_symbol(lvl)
35		}
36
37		## Open new window
38	2x	grDevices::dev.hold()
39	2x	on.exit(grDevices::dev.flush(), add = TRUE)
40	2x	graphics::plot.new()
41
42		## Set plotting coordinates
43	2x	xlim <- range(xy$x)
44	2x	ylim <- range(xy$y)
45	2x	graphics::plot.window(xlim = xlim, ylim = ylim)
46
47		## Plot
48	2x	graphics::points(x = xy$x, y = xy$y, pch = pch, col = col, bg = bg)
49
50		## Construct axis
51	2x	if (axes) {
52	2x	graphics::axis(side = 1, at = seq_len(ncol(x)),
53	2x	labels = labels(x), las = 1,
54	2x	cex.axis = cex.axis, col.axis = col.axis,
55	2x	font.axis = font.axis)
56	2x	graphics::axis(side = 2, las = 1)
57		}
58
59		## Plot frame
60	2x	if (frame.plot) {
61	2x	graphics::box()
62		}
63
64		## Add annotation
65	2x	if (ann) {
66	2x	ylab <- ylab %\|\|% get_transformation(x)
67	2x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
68		}
69
70		## Add legend
71	2x	if (is.list(legend) && is_grouped(x)) {
72		## Compute legend position
73	1x	args <- list(x = "topright", legend = unique(lvl), pch = unique(pch),
74	1x	col = unique(col), bg = unique(bg), bty = "n")
75	1x	args <- utils::modifyList(args, legend)
76
77		## Plot legend
78	1x	do.call(graphics::legend, args = args)
79		}
80
81	2x	invisible(x)
82		}
83
84		#' @export
85		#' @rdname plot
86		#' @aliases plot,LogRatio,missing-method
87		setMethod("plot", c(x = "LogRatio", y = "missing"), plot.LogRatio)

1		# GRAPH
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname as_graph
7		#' @aliases as_graph,LR-method
8		setMethod(
9		f = "as_graph",
10		signature = c(object = "LR"),
11		definition = function(object) {
12		## Validation
13	!	arkhe::assert_package("igraph")
14
15	!	ratio <- labels(object)
16	!	edges <- do.call(rbind, strsplit(ratio, "/"))
17	!	edges <- edges[, c(2, 1)]
18	!	igraph::graph_from_edgelist(edges, directed = FALSE)
19		}
20		)
21
22		#' @export
23		#' @rdname as_graph
24		#' @aliases as_graph,ALR-method
25		setMethod(
26		f = "as_graph",
27		signature = c(object = "ALR"),
28		definition = function(object) {
29		## Validation
30	!	arkhe::assert_package("igraph")
31
32	!	ratio <- labels(object)
33	!	edges <- do.call(rbind, strsplit(ratio, "/"))
34	!	edges <- edges[, c(2, 1)]
35	!	igraph::graph_from_edgelist(edges, directed = TRUE)
36		}
37		)
38
39		#' @export
40		#' @rdname as_graph
41		#' @aliases as_graph,ILR-method
42		setMethod(
43		f = "as_graph",
44		signature = c(object = "ILR"),
45		definition = function(object) {
46		## Validation
47	!	arkhe::assert_package("igraph")
48
49	!	ratio <- labels(object)
50	!	ratio <- gsub(pattern = "[\\(\\)]", replacement = "", x = ratio)
51	!	edges <- lapply(
52	!	X = strsplit(ratio, "/"),
53	!	FUN = function(x) {
54	!	a <- unlist(strsplit(x[[1]], ","))
55	!	b <- unlist(strsplit(x[[2]], ","))
56	!	expand.grid(b, a)
57		}
58		)
59	!	edges <- do.call(rbind, edges)
60	!	igraph::graph_from_data_frame(edges, directed = TRUE)
61		}
62		)

1		# DATA TRANSFORMATION: ISOMETRIC LOG RATIO
2		#' @include AllGenerics.R
3		NULL
4
5		# ILR ==========================================================================
6		ilr_base <- function(D, method = "basic") {
7		## Validation
8	12x	method <- match.arg(method, several.ok = FALSE)
9
10	12x	seq_parts <- seq_len(D - 1)
11
12		## Original ILR transformation defined by Egozcue et al. 2003
13	12x	if (method == "basic") {
14		## Helmert matrix (rotation matrix)
15	12x	H <- stats::contr.helmert(D) # D x D-1
16	12x	H <- t(H) / sqrt((seq_parts + 1) * seq_parts) # D-1 x D
17
18		## Center
19	12x	M <- diag(x = 1, nrow = D) - matrix(data = 1 / D, nrow = D, ncol = D)
20	12x	V <- tcrossprod(M, H)
21		}
22
23	12x	V
24		}
25
26		#' @export
27		#' @rdname transform_ilr
28		#' @aliases transform_ilr,CompositionMatrix,missing-method
29		setMethod(
30		f = "transform_ilr",
31		signature = c(object = "CompositionMatrix"),
32		definition = function(object) {
33	9x	weights <- rep(1 / ncol(object), ncol(object))
34	9x	base <- ilr_base(D = ncol(object), method = "basic")
35	9x	.transform_ilr(object, base, weights)
36		}
37		)
38
39		#' @export
40		#' @rdname transform_ilr
41		#' @aliases transform_ilr,GroupedComposition,missing-method
42		setMethod(
43		f = "transform_ilr",
44		signature = c(object = "GroupedComposition"),
45		definition = function(object) {
46	1x	z <- methods::callNextMethod()
47	1x	.GroupedILR(z, group_indices = group_indices(object),
48	1x	group_levels = group_levels(object))
49		}
50		)
51
52		#' @export
53		#' @rdname transform_ilr
54		#' @aliases transform_ilr,CLR,missing-method
55		setMethod(
56		f = "transform_ilr",
57		signature = c(object = "CLR"),
58		definition = function(object) {
59	2x	weights <- object@weights
60	2x	base <- ilr_base(D = ncol(object), method = "basic")
61	2x	object@.Data <- exp(object@.Data)
62	2x	.transform_ilr(object, base, weights)
63		}
64		)
65
66		#' @export
67		#' @rdname transform_ilr
68		#' @aliases transform_ilr,GroupedCLR,missing-method
69		setMethod(
70		f = "transform_ilr",
71		signature = c(object = "GroupedCLR"),
72		definition = function(object) {
73	!	z <- methods::callNextMethod()
74	!	.GroupedILR(z, group_indices = group_indices(object),
75	!	group_levels = group_levels(object))
76		}
77		)
78
79		#' @export
80		#' @rdname transform_ilr
81		#' @aliases transform_ilr,ALR,missing-method
82		setMethod(
83		f = "transform_ilr",
84		signature = c(object = "ALR"),
85		definition = function(object) {
86	1x	object <- transform_clr(object)
87	1x	methods::callGeneric(object)
88		}
89		)
90
91		#' @export
92		#' @rdname transform_ilr
93		#' @aliases transform_ilr,GroupedALR,missing-method
94		setMethod(
95		f = "transform_ilr",
96		signature = c(object = "GroupedALR"),
97		definition = function(object) {
98	!	z <- methods::callNextMethod()
99	!	.GroupedILR(z, group_indices = group_indices(object),
100	!	group_levels = group_levels(object))
101		}
102		)
103
104		.transform_ilr <- function(object, base, weights) {
105	11x	D <- ncol(object)
106	11x	seq_parts <- seq_len(D - 1)
107	11x	parts <- colnames(object)
108
109		## Rotated and centered values
110	11x	y <- log(object, base = exp(1))
111	11x	ilr <- y %*% base
112
113	11x	ratio <- vapply(
114	11x	X = seq_parts,
115	11x	FUN = function(i, k) {
116	36x	sprintf("(%s)/%s", paste0(k[seq_len(i)], collapse = ","), k[i + 1])
117		},
118	11x	FUN.VALUE = character(1),
119	11x	k = parts
120		)
121	11x	colnames(ilr) <- paste0("Z", seq_parts)
122	11x	rownames(ilr) <- rownames(object)
123
124	11x	.ILR(
125	11x	ilr,
126	11x	parts = parts,
127	11x	ratio = ratio,
128	11x	order = seq_len(D),
129	11x	base = base,
130	11x	weights = weights,
131	11x	totals = totals(object)
132		)
133		}
134
135		# Univariate ILR ===============================================================
136		#' @export
137		#' @rdname univariate_ilr
138		#' @aliases univariate_ilr,numeric-method
139		setMethod(
140		f = "univariate_ilr",
141		signature = c(object = "numeric"),
142		definition = function(object) {
143	12x	sqrt(1 / 2) * log(object / (1 - object))
144		}
145		)
146
147		#' @export
148		#' @rdname univariate_ilr
149		#' @aliases univariate_ilr,matrix-method
150		setMethod(
151		f = "univariate_ilr",
152		signature = c(object = "matrix"),
153		definition = function(object) {
154	2x	apply(X = object, MARGIN = 2, FUN = univariate_ilr)
155		}
156		)

1		# DATA TRANSFORMATION: PIVOT LOG RATIO
2		#' @include AllGenerics.R
3		NULL
4
5		# Pivot ========================================================================
6		#' @export
7		#' @rdname transform_plr
8		#' @aliases transform_plr,CompositionMatrix-method
9		setMethod(
10		f = "transform_plr",
11		signature = c(object = "CompositionMatrix"),
12		definition = function(object, pivot = 1) {
13	5x	J <- ncol(object)
14	5x	parts <- colnames(object)
15
16		## Reorder
17	5x	pivot <- if (is.character(pivot)) which(parts == pivot) else as.integer(pivot)
18	5x	ordering <- c(pivot, which(pivot != seq_len(J)))
19	5x	parts <- parts[ordering]
20	5x	obj <- object[, ordering, drop = FALSE]
21
22	5x	x <- seq_len(J - 1)
23	5x	balances <- diag(sqrt((J - x) / (J - x + 1)))
24	5x	z <- 1 / matrix(data = seq_len(J) - J, nrow = J, ncol = J)
25	5x	z[lower.tri(z)] <- 0
26	5x	diag(z) <- 1
27	5x	z <- z[-nrow(z), ]
28
29	5x	H <- t(balances %*% z)
30	5x	plr <- log(obj, base = exp(1)) %*% H
31
32	5x	ratio <- vapply(
33	5x	X = seq_len(J - 1),
34	5x	FUN = function(i, parts) {
35	16x	j <- length(parts)
36	16x	sprintf("%s/(%s)", parts[1], paste0(parts[(i+1):j], collapse = ","))
37		},
38	5x	FUN.VALUE = character(1),
39	5x	parts = parts
40		)
41	5x	colnames(plr) <- paste0("Z", seq_len(J - 1))
42	5x	rownames(plr) <- rownames(object)
43
44	5x	.PLR(
45	5x	plr,
46	5x	parts = parts,
47	5x	ratio = ratio,
48	5x	order = order(ordering),
49	5x	base = H,
50	5x	weights = rep(1 / J, J),
51	5x	totals = totals(object)
52		)
53		}
54		)
55
56		#' @export
57		#' @rdname transform_plr
58		#' @aliases transform_plr,GroupedComposition-method
59		setMethod(
60		f = "transform_plr",
61		signature = c(object = "GroupedComposition"),
62		definition = function(object, pivot = 1) {
63	!	z <- methods::callNextMethod()
64	!	.GroupedPLR(z, group_indices = group_indices(object),
65	!	group_levels = group_levels(object))
66		}
67		)

1		# CHEMISTRY
2		#' @include AllGenerics.R
3		NULL
4
5		is_chemical <- function(object) {
6	6x	pattern <- "^([A-Z]{1}[a-z]?[1-9]*)+$"
7	6x	z <- grepl(pattern, x = object)
8	6x	names(z) <- object
9	6x	z
10		}
11
12		#' @export
13		#' @rdname chemistry
14		#' @aliases is_oxide,character-method
15		setMethod(
16		f = "is_oxide",
17		signature = c(object = "character"),
18		definition = function(object) {
19	3x	pattern <- "^[A-Z]{1}[a-z]?[1-9][O]{1}[1-9]$"
20	3x	z <- grepl(pattern, x = object)
21	3x	names(z) <- object
22	3x	z
23		}
24		)
25
26		#' @export
27		#' @rdname chemistry
28		#' @aliases is_oxide,CompositionMatrix-method
29		setMethod(
30		f = "is_oxide",
31		signature = c(object = "CompositionMatrix"),
32		definition = function(object) {
33	!	methods::callGeneric(colnames(object))
34		}
35		)
36
37		#' @export
38		#' @rdname chemistry
39		#' @aliases is_element_major,CompositionMatrix-method
40		setMethod(
41		f = "is_element_major",
42		signature = c(object = "CompositionMatrix"),
43		definition = function(object, min = 1 / 100, max = Inf) {
44	1x	.element_threshold(object, min = min, max = max)
45		}
46		)
47
48		#' @export
49		#' @rdname chemistry
50		#' @aliases is_element_minor,CompositionMatrix-method
51		setMethod(
52		f = "is_element_minor",
53		signature = c(object = "CompositionMatrix"),
54		definition = function(object, min = 0.1 / 100, max = 1 / 100) {
55	1x	.element_threshold(object, min = min, max = max)
56		}
57		)
58
59		#' @export
60		#' @rdname chemistry
61		#' @aliases is_element_trace,CompositionMatrix-method
62		setMethod(
63		f = "is_element_trace",
64		signature = c(object = "CompositionMatrix"),
65		definition = function(object, min = -Inf, max = 0.1 / 100) {
66	!	.element_threshold(object, min = min, max = max)
67		}
68		)
69
70		.element_threshold <- function(x, min = -Inf, max = Inf) {
71	2x	x <- mean(x)
72	2x	x >= min & x < max
73		}

1		# PAIRS
2		#' @include AllGenerics.R
3		NULL
4
5		# CompositionMatrix ============================================================
6		#' @export
7		#' @method pairs CompositionMatrix
8		pairs.CompositionMatrix <- function(x, margin = NULL, ...) {
9	1x	isopleuros::ternary_pairs(x, margin = margin, ...)
10	1x	invisible(x)
11		}
12
13		#' @export
14		#' @rdname pairs
15		#' @aliases pairs,CompositionMatrix-method
16		setMethod("pairs", c(x = "CompositionMatrix"), pairs.CompositionMatrix)
17
18		#' @export
19		#' @method pairs GroupedComposition
20		pairs.GroupedComposition <- function(x, ..., margin = NULL,
21		palette_color = palette_color_discrete(),
22		palette_symbol = palette_shape()) {
23		## Aesthetics
24	1x	lvl <- group_names(x)
25	1x	col <- palette_color(lvl)
26	1x	pch <- palette_symbol(lvl)
27
28	1x	isopleuros::ternary_pairs(x, margin = margin, col = col, pch = pch, ...)
29	1x	invisible(x)
30		}
31
32		#' @export
33		#' @rdname pairs
34		#' @aliases pairs,GroupedComposition-method
35		setMethod("pairs", c(x = "GroupedComposition"), pairs.GroupedComposition)

1		# SPLIT
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @method split CompositionMatrix
7		split.CompositionMatrix <- function(x, f, drop = FALSE, ...) {
8	!	lapply(
9	!	X = split(x = seq_len(nrow(x)), f = f, drop = drop, sep = "_", ...),
10	!	FUN = function(ind) x[ind, , drop = FALSE]
11		)
12		}
13
14		#' @export
15		#' @rdname group_split
16		#' @aliases split,CompositionMatrix-method
17		setMethod("split", "CompositionMatrix", split.CompositionMatrix)
18
19		#' @export
20		#' @method split LogRatio
21		split.LogRatio <- function(x, f, drop = FALSE, ...) {
22	!	lapply(
23	!	X = split(x = seq_len(nrow(x)), f = f, drop = drop, sep = "_", ...),
24	!	FUN = function(ind) x[ind, , drop = FALSE]
25		)
26		}
27
28		#' @export
29		#' @rdname group_split
30		#' @aliases split,LogRatio-method
31		setMethod("split", "LogRatio", split.LogRatio)
32
33		#' @export
34		#' @rdname group_split
35		#' @aliases group_split,GroupedComposition-method
36		setMethod(
37		f = "group_split",
38		signature = "GroupedComposition",
39		definition = function(object, ...) {
40	6x	lapply(
41	6x	X = group_rows(object),
42	6x	FUN = function(ind) ungroup(object[ind, , drop = FALSE])
43		)
44		}
45		)
46
47		#' @export
48		#' @rdname group_split
49		#' @aliases group_split,GroupedLogRatio-method
50		setMethod(
51		f = "group_split",
52		signature = "GroupedLogRatio",
53		definition = function(object, ...) {
54	!	lapply(
55	!	X = group_rows(object),
56	!	FUN = function(ind) ungroup(object[ind, , drop = FALSE])
57		)
58		}
59		)