nexus coverage - 84.66%

Files
Source

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

# 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 labels
#' @aliases labels,CompositionMatrix-method
setMethod("labels", "CompositionMatrix", labels.CompositionMatrix)

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

#' @export
#' @rdname labels
#' @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 weights
#' @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 weights
#' @aliases weights,LR-method
setMethod("weights", "LR", weights.LR)

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

#' @export
#' @rdname weights
#' @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
  }
)

# 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 <- closure(x * y)
    z <- methods::initialize(x, z)
    rownames(z) <- rownames(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 <- closure(x ^ y)
    z <- methods::initialize(x, z)
    rownames(z) <- rownames(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)
}

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

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

# BOXPLOT
#' @include AllGenerics.R
NULL

# LogRatio =====================================================================
#' @export
#' @method boxplot LogRatio
boxplot.LogRatio <- function(x, ..., range = 1.5, width = NULL, varwidth = FALSE,
                             notch = FALSE, outline = TRUE,
                             plot = TRUE, horizontal = FALSE,
                             xlab = NULL, ylab = NULL, main = NULL, sub = NULL,
                             ann = graphics::par("ann")) {
  z <- as.data.frame(x)
  if (horizontal) {
    xlab <- xlab %||% get_transformation(x)
  } else {
    ylab <- ylab %||% get_transformation(x)
  }
  box <- graphics::boxplot(z, ..., range = range, width = width,
                           varwidth = varwidth, notch = notch, outline = outline,
                           names = labels(x), plot = plot, horizontal = horizontal,
                           xlab = xlab, ylab = ylab, main = main, sub = sub,
                           ann = ann, log = "", las = 1)
  if (!plot) return(invisible(box))

  invisible(x)
}

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

# GroupedLogRatio ==============================================================
#' @export
#' @method boxplot GroupedLogRatio
boxplot.GroupedLogRatio <- function(x, ..., range = 1.5, width = NULL,
                                    varwidth = FALSE, notch = FALSE,
                                    outline = TRUE, plot = TRUE,
                                    horizontal = FALSE, color = NULL,
                                    xlab = NULL, ylab = NULL,
                                    main = NULL, sub = NULL,
                                    ann = graphics::par("ann"),
                                    legend = list(x = "topright")) {
  ## Graphical parameters
  lvl <- group_levels(x)
  col <- khroma::palette_color_discrete(color)(lvl)
  bg <- grDevices::adjustcolor(col, alpha.f = 0.5)
  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")

  ## Prepare data
  df <- data.frame(
    x = as.numeric(x),
    y = interaction(
      factor(rep(labels(x), each = nrow(x)), levels = labels(x)),
      factor(rep(group_names(x), times = ncol(x)), levels = group_levels(x)),
      lex.order = TRUE
    )
  )

  ## Plot
  box <- graphics::boxplot(
    formula = x ~ y, data = df, ...,
    range = range,
    width = width,
    varwidth = varwidth,
    notch = notch,
    outline = outline,
    plot = plot,
    border = col,
    col = bg,
    horizontal = horizontal,
    log = "",
    ann = FALSE,
    xaxt = ifelse(!horizontal, "n", "s"),
    yaxt = ifelse(horizontal, "n", "s"),
    las = 1
  )
  if (!plot) return(invisible(box))

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

  ## Construct axis
  n_group <- length(group_levels(x))
  graphics::axis(
    side = ifelse(horizontal, 2, 1),
    at = seq(from = (n_group + 1) / 2, to = ncol(x) * n_group, by = n_group),
    labels = labels(x),
    cex.axis = cex.axis,
    col.axis = col.axis,
    font.axis = font.axis,
    las = 1
  )

  ## Add legend
  if (is.list(legend) && length(legend) > 0) {
    args <- list(x = "topright", legend = lvl, fill = col, bty = "n")
    args <- utils::modifyList(args, legend)
    do.call(graphics::legend, args = args)
  }

  invisible(x)
}

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

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

# 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 SUMMARY: DESCRIBE
#' @include AllGenerics.R
NULL

#' @export
#' @rdname describe
#' @aliases describe,CompositionMatrix-method
setMethod(
  f = "describe",
  signature = c(x = "CompositionMatrix"),
  definition = function(x) {
    ## Variables
    msg_parts <- describe_coda(x)

    ## Missing values
    msg_miss <- describe_missing(x)

    ## Check
    msg_val <- describe_check(x)

    cat(msg_parts, msg_miss, msg_val, sep = "\n")

    invisible(x)
  }
)

#' @export
#' @rdname describe
#' @aliases describe,GroupedComposition-method
setMethod(
  f = "describe",
  signature = c(x = "GroupedComposition"),
  definition = function(x) {
    ## Variables
    msg_parts <- describe_coda(x)

    ## Groups
    msg_group <- describe_groups(x)

    ## Missing values
    msg_miss <- describe_missing(x)

    ## Check
    msg_val <- describe_check(x)

    cat(paste0(msg_parts, msg_group), msg_miss, msg_val, sep = "\n")

    invisible(x)
  }
)

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

  rows <- sprintf(ngettext(m, "%d composition", "%d compositions"), m)
  title <- sprintf("%s:", rows)

  cols <- paste0(dQuote(labels(x)), collapse = ", ")
  msg <- sprintf(ngettext(p, "%d part", "%d parts"), p)
  msg <- sprintf("\n* %s: %s.", msg, cols)

  paste0(title, msg, collapse = "")
}
describe_groups <- function(x) {

  i <- group_n(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)

  paste0(sprintf("\n* %s.", c(msg_grp, msg_ung)), collapse = "")
}
describe_missing <- function(x) {
  m <- nrow(x)
  p <- ncol(x)

  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)

  title <- sprintf(ngettext(n_NA, "%d missing value:", "%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]])

  msg <- paste0(sprintf("\n* %s.", c(msg_row_NA, msg_col_NA)), collapse = "")
  paste0("\n", title, msg, collapse = "")
}
describe_check <- function(x) {
  title <- tr_("Data checking:")

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

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

  msg <- paste0(sprintf("\n* %s.", c(msg_spa, msg_col_var)), collapse = "")
  paste0("\n", title, msg, collapse = "")
}

# PLOT
#' @include AllGenerics.R
NULL

# LogRatio =====================================================================
#' @export
#' @method plot LogRatio
plot.LogRatio <- function(x, ..., jitter_factor = 1, jitter_amount = NULL,
                          xlab = NULL, ylab = NULL,
                          main = NULL, sub = NULL, ann = graphics::par("ann"),
                          axes = TRUE, frame.plot = axes) {
  ## Get data
  xy <- data.frame(
    x = jitter(as.vector(col(x)), factor = jitter_factor, amount = jitter_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")

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

  invisible(x)
}

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

#' @export
#' @method plot GroupedLogRatio
plot.GroupedLogRatio <- function(x, ..., jitter_factor = 1, jitter_amount = NULL,
                                 color = NULL, symbol = NULL,
                                 xlab = NULL, ylab = NULL,
                                 main = NULL, sub = NULL,
                                 ann = graphics::par("ann"),
                                 axes = TRUE, frame.plot = axes,
                                 legend = list(x = "topright")) {

  ## Graphical parameters
  lvl <- group_names(x)
  col <- khroma::palette_color_discrete(color)(lvl)
  bg <- grDevices::adjustcolor(col, alpha.f = 0.5)
  pch <- khroma::palette_shape(symbol)(lvl)

  ## Plot
  plot(
    ungroup(x),
    col = col, pch = pch, bg = bg,
    jitter_factor = jitter_factor,
    jitter_amount = jitter_amount,
    xlab = xlab, ylab = ylab,
    main = main, sub = sub, ann = ann,
    axes = axes, frame.plot = frame.plot
  )

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

  invisible(x)
}

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

# 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),
                group_ordered = is_ordered(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),
                group_ordered = is_ordered(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),
                group_ordered = is_ordered(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)
  }
)

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


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

# 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,
                                      color = NULL, 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 <- khroma::palette_color_discrete(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
  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(mar = mar)
  on.exit(graphics::par(old_par), add = TRUE)

  ## 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,
                            verbose = getOption("nexus.verbose")) {
  ## Remove missing values
  x <- arkhe::remove_NA(x, margin = 2, verbose = verbose)

  ## 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_n(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
  )
}

# 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, exclude = NULL)[i])
    methods::initialize(
      x, z,
      totals = totals,
      group_indices = as.integer(g),
      group_levels = levels(g),
      group_ordered = is.ordered(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)
  }
)

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

# 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) [`factor`]s 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 <- if (length(x) > 1) interaction(x, sep = "_") else x[[1L]]
  }
  x <- as.factor(x)
  if (drop_levels) x <- droplevels(x)
  if (allow_na) x <- addNA(x, ifany = TRUE)

  x
}

#' Validate Groups
#'
#' @param object A [`matrix`]-like object.
#' @param by A ([`list`] of) [`factor`]s for which interaction is to be computed.
#' @param verbose A [`logical`] scalar: should \R report extra information
#'  on progress?
#' @return Invisibly returns `by`.
#' @keywords internal
#' @noRd
validate_groups <- function(object, by, verbose = getOption("nexus.verbose")) {
  arkhe::assert_type(by, "integer")
  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))
  }

  invisible(by)
}

#' @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
    validate_groups(object, by, verbose = verbose)

    .GroupedComposition(
      object,
      group_indices = as.integer(by),
      group_levels = levels(by),
      group_ordered = is.ordered(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, exclude = NULL)), 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
#' @rdname group_names
#' @aliases group_levels,ReferenceGroups-method
setMethod(
  f = "group_levels",
  signature = "ReferenceGroups",
  definition = function(object) object@group_levels
)

#' @export
#' @rdname group_names
#' @aliases group_names,ReferenceGroups-method
setMethod(
  f = "group_names",
  signature = "ReferenceGroups",
  definition = function(object) group_levels(object)[group_indices(object)]
)

is_ordered <- function(object) {
  isTRUE(object@group_ordered)
}

#' @export
#' @rdname group_names
#' @aliases group_factor,ReferenceGroups-method
setMethod(
  f = "group_factor",
  signature = "ReferenceGroups",
  definition = function(object, exclude = NA) {
    factor(
      x = group_names(object),
      levels = group_levels(object),
      exclude = exclude,
      ordered = is_ordered(object)
    )
  }
)

#' @export
#' @rdname group_names
#' @aliases group_indices,ReferenceGroups-method
setMethod(
  f = "group_indices",
  signature = "ReferenceGroups",
  definition = function(object) object@group_indices
)

#' @export
#' @rdname group_names
#' @aliases group_rows,ReferenceGroups-method
setMethod(
  f = "group_rows",
  signature = "ReferenceGroups",
  definition = function(object) {
    i <- group_factor(object, exclude = NULL)
    split(seq_along(i), f = i)
  }
)

#' @export
#' @rdname group_names
#' @aliases group_n,ReferenceGroups-method
setMethod(
  f = "group_n",
  signature = "ReferenceGroups",
  definition = function(object) length(group_levels(object))
)

#' @export
#' @rdname group_names
#' @aliases group_size,ReferenceGroups-method
setMethod(
  f = "group_size",
  signature = "ReferenceGroups",
  definition = function(object) lengths(group_rows(object))
)

# Predicates ===================================================================
#' @export
#' @rdname is_assigned
#' @aliases is_assigned,ReferenceGroups-method
setMethod(
  f = "is_assigned",
  signature = "ReferenceGroups",
  definition = function(object) !is.na(group_names(object))
)

#' @export
#' @rdname is_assigned
#' @aliases any_assigned,ReferenceGroups-method
setMethod(
  f = "any_assigned",
  signature = "ReferenceGroups",
  definition = function(object) any(is_assigned(object))
)

#' @export
#' @rdname is_assigned
#' @aliases all_assigned,ReferenceGroups-method
setMethod(
  f = "all_assigned",
  signature = "ReferenceGroups",
  definition = function(object) all(is_assigned(object))
)

# 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]
    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, ..., group_var = ".group", group_after = 0) {
  z <- as.data.frame(methods::as(x, "matrix"), row.names = rownames(x), ...)

  if (!is_grouped(x)) return(z)
  arkhe::append_column(z, group_names(x), after = group_after, var = group_var)
}

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

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

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

  if (!is_grouped(x)) return(z)
  arkhe::append_column(z, group_names(x), after = group_after, var = group_var)
}

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

#' @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)

# 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_scalar(scale, "numeric")

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

  y
}

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

# 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),
                group_ordered = is_ordered(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),
                group_ordered = is_ordered(object))
  }
)

# 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),
                        group_ordered = is_ordered(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),
                        group_ordered = is_ordered(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),
                        group_ordered = is_ordered(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
  }
)

# 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),
                group_ordered = is_ordered(object))
  }
)

# 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 --------------------------------------------------------------------
#' Weights
#'
#' Extract weights from an object.
#' @param object An object for which the extraction of weights is meaningful.
# @param value A possible value for the element(s) of `object`.
#' @param ... Currently not used.
# @example inst/examples/ex-mutators.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name weights
#' @rdname weights
#' @aliases get set
NULL

#' Labels
#'
#' Find a suitable set of labels from an object.
#' @param object An object from which to get the labels.
# @param value A possible value for the element(s) of `object`.
#' @param ... Currently not used.
# @example inst/examples/ex-mutators.R
#' @author N. Frerebeau
#' @docType methods
#' @family mutators
#' @name labels
#' @rdname labels
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 `object`.
#' @return
#'  * `totals() <- value` returns an object of the same sort as `object` with
#'    the new row sums assigned.
#'  * `totals()` returns the row sums of `object`.
#' @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 ----------------------------------------------------------------------
#' 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 (only used if `x` is [grouped][group()]).
#' @param group_after An [`integer`] specifying a subscript, after which the new
#'  `group_var` column is to be appended (only used if `x` is [grouped][group()]).
#' @param ... Further parameters to be passed to
#'  [`as.data.frame()`][base::as.data.frame].
#' @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 =======================================================================
#' Group by One or More Variables
#'
#' 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")
)

#' Grouped Data
#'
#' Retrieve the (reference) groups to which the observations belong.
#' @param object A [grouped][group()] \R object.
#' @return
#'  * `is_assigned()` returns a [`logical`] vector specifying whether or
#'    not an observation belongs to a group.
#'  * `any_assigned()` returns an [`logical`] scalar specifying if any
#'    observation belongs to a group.
#'  * `all_assigned()` returns an [`logical`] scalar specifying if all
#'    observations belong to a group.
#' @example inst/examples/ex-group.R
#' @author N. Frerebeau
#' @docType methods
#' @family grouping methods
#' @aliases is_assigned-method
setGeneric(
  name = "is_assigned",
  def = function(object) standardGeneric("is_assigned")
)

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

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

#' Grouping Metadata
#'
#' Retrieve the (reference) groups to which the observations belong.
#' @param object A [grouped][group()] \R object.
#' @param exclude A [`character`] vector of values to be excluded when forming
#'  the set of levels.
#' @param ... Currently not used.
#' @return
#'  * `group_levels()` returns a [`character`] vector giving the group
#'    names.
#'  * `group_size()` returns an [`integer`] vector giving the size of each
#'    group.
#'  * `group_n()` gives the total number of groups.
#'  * `group_names()` returns a [`character`] vector giving the name of
#'    the group that each observation belongs to.
#'  * `group_factor()` returns a [`factor`] vector giving the name of
#'    the group that each observation belongs to.
#'  * `group_indices()` returns an [`integer`] vector giving the group
#'    that each value belongs to.
#'  * `group_rows()` returns a `list` of [`integer`] vectors giving the
#'    observation that each group contains.
#' @example inst/examples/ex-group.R
#' @author N. Frerebeau
#' @docType methods
#' @family grouping methods
#' @aliases group_names-method
setGeneric(
  name = "group_names",
  def = function(object) standardGeneric("group_names")
)

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

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

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

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

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

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

#' 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 color A (named) vector of colors (will be mapped to the group names
#'  of `object`) or a [`function`] that when called with a single argument (an
#'  integer specifying the number of colors) returns a 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 color A (named) vector of colors (will be mapped to the group names
#'  of `object`) or a [`function`] that when called with a single argument (an
#'  integer specifying the number of colors) returns a vector of colors.
#' @param symbol A (named) vector of colors (will be mapped to the group names
#'  of `object`).
#' @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

## Scatter plot ----------------------------------------------------------------
#' Plot Log-Ratios
#'
#' Displays a scatter plot.
#' @param x A [`LogRatio-class`] object.
#' @param jitter_factor,jitter_amount A length-one [`numeric`] vector specifying
#'  the amount of jitter (see [jitter()]).
#' @param color A (named) vector of colors (will be mapped to the group names
#'  of `object`) or a [`function`] that when called with a single argument (an
#'  integer specifying the number of colors) returns a vector of colors.
#' @param symbol A (named) vector of colors (will be mapped to the group names
#'  of `object`).
#' @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

## Boxplot ---------------------------------------------------------------------
#' Boxplot of Log-Ratios
#'
#' Displays box-and-whisker plots of the given (grouped) values.
#' @param x A [`LogRatio-class`] object.
#' @param range A length-one [`numeric`] vector specifying how far the plot
#'  whiskers extend out from the box (see [graphics::boxplot()]).
#' @param width A [`numeric`] vector giving the relative widths of the boxes
#'  making up the plot.
#' @param varwidth A [`logical`] scalar: should the boxes be drawn with widths
#'  proportional to the square-roots of the number of observations in the
#'  groups?
#' @param notch A [`logical`] scalar: should a notch be drawn in each side of
#'  the boxes?
#' @param outline A [`logical`] scalar: should the outliers be drawn?
#' @param plot A [`logical`] scalar: should a boxplot be produced? If `FALSE`,
#'  the summaries which the boxplots are based on are invisibly returned.
#' @param horizontal A [`logical`] scalar: should the boxplots be horizontal?
#' @param color A (named) vector of colors (will be mapped to the group names
#'  of `object`) or a [`function`] that when called with a single argument (an
#'  integer specifying the number of colors) returns a vector of colors.
#' @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 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
#'  `boxplot()` is called for its side-effects: is results in a graphic being
#'  displayed (invisibly return `x`).
#' @seealso [graphics::boxplot()]
#' @example inst/examples/ex-boxplot.R
#' @author N. Frerebeau
#' @docType methods
#' @family plot methods
#' @name boxplot
#' @rdname boxplot
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"
)

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

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

# 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(as.matrix(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)
  }
)

# 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, verbose = verbose)
    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, exclude = NULL)
    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, exclude = NULL),
      FUN = mean
    )

    z <- .CompositionMatrix(z, totals = as.numeric(tot))
    group(z, by = flatten_chr(x = grp, by = group_factor(x, exclude = NULL)),
          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
}

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

# 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),
                group_ordered = is_ordered(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),
                group_ordered = is_ordered(object))
  }
)

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

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

# 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),
               group_ordered = is_ordered(object))
  }
)

# 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) {
    arkhe::assert_colnames(y, colnames(x))

    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) {
    z <- methods::callNextMethod(x, y)
    group(z, by = c(group_names(x), group_names(y)))
  }
)

# 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,
                                     color = NULL, symbol = NULL) {
  ## Aesthetics
  lvl <- group_names(x)
  col <- khroma::palette_color_discrete(color)(lvl)
  pch <- khroma::palette_shape(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])
    )
  }
)

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

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

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

1		# AGGREGATE
2		#' @include AllGenerics.R
3		NULL
4
5		# Aggregate ====================================================================
6		#' @export
7		#' @method aggregate CompositionMatrix
8		aggregate.CompositionMatrix <- function(x, by, FUN, ...,
9		simplify = TRUE, drop = TRUE) {
10	1x	x <- group(x, by = by, drop_levels = drop, verbose = FALSE)
11	1x	aggregate(x, FUN, ..., simplify = simplify)
12		}
13
14		#' @export
15		#' @rdname aggregate
16		#' @aliases aggregate,CompositionMatrix-method
17		setMethod("aggregate", "CompositionMatrix", aggregate.CompositionMatrix)
18
19		#' @export
20		#' @method aggregate GroupedComposition
21		aggregate.GroupedComposition <- function(x, FUN, ..., simplify = TRUE) {
22		## Grouping
23	6x	aggr <- lapply(
24	6x	X = group_rows(x),
25	6x	FUN = function(i, data, fun, ...) fun(data[i, , drop = FALSE], ...),
26	6x	data = x,
27	6x	fun = FUN,
28		...
29		)
30
31	6x	has_dim <- vapply(
32	6x	X = aggr,
33	6x	FUN = function(x) !is.null(nrow(x)) && nrow(x) > 1,
34	6x	FUN.VALUE = logical(1)
35		)
36
37	!	if (any(has_dim) \|\| !simplify) return(aggr)
38	6x	aggr <- do.call(rbind, aggr)
39	6x	rownames(aggr) <- group_levels(x)
40	6x	aggr
41		}
42
43		#' @export
44		#' @rdname aggregate
45		#' @aliases aggregate,GroupedComposition-method
46		setMethod("aggregate", "GroupedComposition", aggregate.GroupedComposition)

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

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		# 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		# BOXPLOT
2		#' @include AllGenerics.R
3		NULL
4
5		# LogRatio =====================================================================
6		#' @export
7		#' @method boxplot LogRatio
8		boxplot.LogRatio <- function(x, ..., range = 1.5, width = NULL, varwidth = FALSE,
9		notch = FALSE, outline = TRUE,
10		plot = TRUE, horizontal = FALSE,
11		xlab = NULL, ylab = NULL, main = NULL, sub = NULL,
12		ann = graphics::par("ann")) {
13	1x	z <- as.data.frame(x)
14	1x	if (horizontal) {
15	!	xlab <- xlab %\|\|% get_transformation(x)
16		} else {
17	1x	ylab <- ylab %\|\|% get_transformation(x)
18		}
19	1x	box <- graphics::boxplot(z, ..., range = range, width = width,
20	1x	varwidth = varwidth, notch = notch, outline = outline,
21	1x	names = labels(x), plot = plot, horizontal = horizontal,
22	1x	xlab = xlab, ylab = ylab, main = main, sub = sub,
23	1x	ann = ann, log = "", las = 1)
24	!	if (!plot) return(invisible(box))
25
26	1x	invisible(x)
27		}
28
29		#' @export
30		#' @rdname boxplot
31		#' @aliases boxplot,LogRatio-method
32		setMethod("boxplot", c(x = "LogRatio"), boxplot.LogRatio)
33
34		# GroupedLogRatio ==============================================================
35		#' @export
36		#' @method boxplot GroupedLogRatio
37		boxplot.GroupedLogRatio <- function(x, ..., range = 1.5, width = NULL,
38		varwidth = FALSE, notch = FALSE,
39		outline = TRUE, plot = TRUE,
40		horizontal = FALSE, color = NULL,
41		xlab = NULL, ylab = NULL,
42		main = NULL, sub = NULL,
43		ann = graphics::par("ann"),
44		legend = list(x = "topright")) {
45		## Graphical parameters
46	1x	lvl <- group_levels(x)
47	1x	col <- khroma::palette_color_discrete(color)(lvl)
48	1x	bg <- grDevices::adjustcolor(col, alpha.f = 0.5)
49	1x	cex.axis <- list(...)$cex.axis %\|\|% graphics::par("cex.axis")
50	1x	col.axis <- list(...)$col.axis %\|\|% graphics::par("col.axis")
51	1x	font.axis <- list(...)$font.axis %\|\|% graphics::par("font.axis")
52
53		## Prepare data
54	1x	df <- data.frame(
55	1x	x = as.numeric(x),
56	1x	y = interaction(
57	1x	factor(rep(labels(x), each = nrow(x)), levels = labels(x)),
58	1x	factor(rep(group_names(x), times = ncol(x)), levels = group_levels(x)),
59	1x	lex.order = TRUE
60		)
61		)
62
63		## Plot
64	1x	box <- graphics::boxplot(
65	1x	formula = x ~ y, data = df, ...,
66	1x	range = range,
67	1x	width = width,
68	1x	varwidth = varwidth,
69	1x	notch = notch,
70	1x	outline = outline,
71	1x	plot = plot,
72	1x	border = col,
73	1x	col = bg,
74	1x	horizontal = horizontal,
75	1x	log = "",
76	1x	ann = FALSE,
77	1x	xaxt = ifelse(!horizontal, "n", "s"),
78	1x	yaxt = ifelse(horizontal, "n", "s"),
79	1x	las = 1
80		)
81	!	if (!plot) return(invisible(box))
82
83		## Add annotation
84	1x	if (ann) {
85	1x	if (horizontal) {
86	!	xlab <- xlab %\|\|% get_transformation(x)
87		} else {
88	1x	ylab <- ylab %\|\|% get_transformation(x)
89		}
90	1x	graphics::title(main = main, sub = sub, xlab = xlab, ylab = ylab)
91		}
92
93		## Construct axis
94	1x	n_group <- length(group_levels(x))
95	1x	graphics::axis(
96	1x	side = ifelse(horizontal, 2, 1),
97	1x	at = seq(from = (n_group + 1) / 2, to = ncol(x) * n_group, by = n_group),
98	1x	labels = labels(x),
99	1x	cex.axis = cex.axis,
100	1x	col.axis = col.axis,
101	1x	font.axis = font.axis,
102	1x	las = 1
103		)
104
105		## Add legend
106	1x	if (is.list(legend) && length(legend) > 0) {
107	1x	args <- list(x = "topright", legend = lvl, fill = col, bty = "n")
108	1x	args <- utils::modifyList(args, legend)
109	1x	do.call(graphics::legend, args = args)
110		}
111
112	1x	invisible(x)
113		}
114
115		#' @export
116		#' @rdname boxplot
117		#' @aliases boxplot,GroupedLogRatio-method
118		setMethod("boxplot", c(x = "GroupedLogRatio"), boxplot.GroupedLogRatio)

1		# SHOW
2		#' @include AllGenerics.R
3		NULL
4
5		setMethod(
6		f = "show",
7		signature = "CompositionMatrix",
8		definition = function(object) {
9	4x	m <- nrow(object)
10	4x	p <- ncol(object)
11	4x	mtx <- methods::as(object, "matrix")
12
13	4x	txt_dim <- sprintf("%d x %d", m, p)
14	4x	txt_mtx <- sprintf("<%s: %s>", class(object), txt_dim)
15
16	4x	cat(
17	4x	txt_mtx,
18	4x	utils::capture.output(mtx),
19	4x	sep = "\n"
20		)
21	4x	invisible(object)
22		}
23		)
24
25		setMethod(
26		f = "show",
27		signature = "LogRatio",
28		definition = function(object) {
29	!	m <- nrow(object)
30	!	p <- ncol(object)
31	!	mtx <- methods::as(object, "matrix")
32
33	!	txt_dim <- sprintf("%d x %d", m, p)
34	!	txt_mtx <- sprintf("<%s: %s>", get_transformation(object), txt_dim)
35
36	!	cat(
37	!	txt_mtx,
38	!	utils::capture.output(mtx),
39	!	sep = "\n"
40		)
41	!	invisible(object)
42		}
43		)

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 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		## Variables
13	1x	msg_parts <- describe_coda(x)
14
15		## Missing values
16	1x	msg_miss <- describe_missing(x)
17
18		## Check
19	1x	msg_val <- describe_check(x)
20
21	1x	cat(msg_parts, msg_miss, msg_val, sep = "\n")
22
23	1x	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		## Variables
35	1x	msg_parts <- describe_coda(x)
36
37		## Groups
38	1x	msg_group <- describe_groups(x)
39
40		## Missing values
41	1x	msg_miss <- describe_missing(x)
42
43		## Check
44	1x	msg_val <- describe_check(x)
45
46	1x	cat(paste0(msg_parts, msg_group), msg_miss, msg_val, sep = "\n")
47
48	1x	invisible(x)
49		}
50		)
51
52		describe_coda <- function(x) {
53	2x	m <- nrow(x)
54	2x	p <- ncol(x)
55
56	2x	rows <- sprintf(ngettext(m, "%d composition", "%d compositions"), m)
57	2x	title <- sprintf("%s:", rows)
58
59	2x	cols <- paste0(dQuote(labels(x)), collapse = ", ")
60	2x	msg <- sprintf(ngettext(p, "%d part", "%d parts"), p)
61	2x	msg <- sprintf("\n* %s: %s.", msg, cols)
62
63	2x	paste0(title, msg, collapse = "")
64		}
65		describe_groups <- function(x) {
66
67	1x	i <- group_n(x)
68	1x	ls_grp <- paste0(dQuote(group_levels(x)), collapse = ", ")
69	1x	msg_grp <- sprintf(ngettext(i, "%d group", "%d groups"), i)
70	1x	msg_grp <- sprintf("%s: %s", msg_grp, ls_grp)
71
72	1x	j <- sum(!is_assigned(x))
73	1x	msg_ung <- sprintf(ngettext(j, "%d unassigned sample", "%d unassigned samples"), j)
74
75	1x	paste0(sprintf("\n* %s.", c(msg_grp, msg_ung)), collapse = "")
76		}
77		describe_missing <- function(x) {
78	2x	m <- nrow(x)
79	2x	p <- ncol(x)
80
81	2x	n_NA <- sum(count(x, f = is.na))
82	2x	m_NA <- sum(detect(x, f = is.na, margin = 1))
83	2x	p_NA <- sum(detect(x, f = is.na, margin = 2))
84	2x	pc <- label_percent(c(m_NA / m, p_NA / p), digits = 1, trim = TRUE)
85
86	2x	title <- sprintf(ngettext(n_NA, "%d missing value:", "%d missing values:"), n_NA)
87
88	2x	rows_NA <- ngettext(m_NA, "%d observation (%s) contains missing values",
89	2x	"%d observations (%s) contain missing values")
90	2x	msg_row_NA <- sprintf(rows_NA, m_NA, pc[[1]])
91
92	2x	cols_NA <- ngettext(p_NA, "%d variable (%s) contains missing values",
93	2x	"%d variables (%s) contain missing values")
94	2x	msg_col_NA <- sprintf(cols_NA, p_NA, pc[[2]])
95
96	2x	msg <- paste0(sprintf("\n* %s.", c(msg_row_NA, msg_col_NA)), collapse = "")
97	2x	paste0("\n", title, msg, collapse = "")
98		}
99		describe_check <- function(x) {
100	2x	title <- tr_("Data checking:")
101
102		## Constant columns
103	2x	p_var <- sum(detect(x, f = function(x) is_unique(x), margin = 2))
104	2x	cols_var <- ngettext(p_var, "%d variable with no variance",
105	2x	"%d variables with no variance")
106	2x	msg_col_var <- sprintf(cols_var, p_var)
107
108		## Sparsity
109	2x	spa <- sparsity(x, count = FALSE)
110	2x	msg_spa <- sprintf(tr_("%s of numeric values are zero"), label_percent(spa, digits = 1))
111
112	2x	msg <- paste0(sprintf("\n* %s.", c(msg_spa, msg_col_var)), collapse = "")
113	2x	paste0("\n", title, msg, collapse = "")
114		}

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

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	1x	group_ordered = is_ordered(object))
50		}
51		)
52
53		#' @export
54		#' @rdname transform_ilr
55		#' @aliases transform_ilr,CLR,missing-method
56		setMethod(
57		f = "transform_ilr",
58		signature = c(object = "CLR"),
59		definition = function(object) {
60	2x	weights <- object@weights
61	2x	base <- ilr_base(D = ncol(object), method = "basic")
62	2x	object@.Data <- exp(object@.Data)
63	2x	.transform_ilr(object, base, weights)
64		}
65		)
66
67		#' @export
68		#' @rdname transform_ilr
69		#' @aliases transform_ilr,GroupedCLR,missing-method
70		setMethod(
71		f = "transform_ilr",
72		signature = c(object = "GroupedCLR"),
73		definition = function(object) {
74	!	z <- methods::callNextMethod()
75	!	.GroupedILR(z, group_indices = group_indices(object),
76	!	group_levels = group_levels(object),
77	!	group_ordered = is_ordered(object))
78		}
79		)
80
81		#' @export
82		#' @rdname transform_ilr
83		#' @aliases transform_ilr,ALR,missing-method
84		setMethod(
85		f = "transform_ilr",
86		signature = c(object = "ALR"),
87		definition = function(object) {
88	1x	object <- transform_clr(object)
89	1x	methods::callGeneric(object)
90		}
91		)
92
93		#' @export
94		#' @rdname transform_ilr
95		#' @aliases transform_ilr,GroupedALR,missing-method
96		setMethod(
97		f = "transform_ilr",
98		signature = c(object = "GroupedALR"),
99		definition = function(object) {
100	!	z <- methods::callNextMethod()
101	!	.GroupedILR(z, group_indices = group_indices(object),
102	!	group_levels = group_levels(object),
103	!	group_ordered = is_ordered(object))
104		}
105		)
106
107		.transform_ilr <- function(object, base, weights) {
108	11x	D <- ncol(object)
109	11x	seq_parts <- seq_len(D - 1)
110	11x	parts <- colnames(object)
111
112		## Rotated and centered values
113	11x	y <- log(object, base = exp(1))
114	11x	ilr <- y %*% base
115
116	11x	ratio <- vapply(
117	11x	X = seq_parts,
118	11x	FUN = function(i, k) {
119	36x	sprintf("(%s)/%s", paste0(k[seq_len(i)], collapse = ","), k[i + 1])
120		},
121	11x	FUN.VALUE = character(1),
122	11x	k = parts
123		)
124	11x	colnames(ilr) <- paste0("Z", seq_parts)
125	11x	rownames(ilr) <- rownames(object)
126
127	11x	.ILR(
128	11x	ilr,
129	11x	parts = parts,
130	11x	ratio = ratio,
131	11x	order = seq_len(D),
132	11x	base = base,
133	11x	weights = weights,
134	11x	totals = totals(object)
135		)
136		}
137
138		# Univariate ILR ===============================================================
139		#' @export
140		#' @rdname univariate_ilr
141		#' @aliases univariate_ilr,numeric-method
142		setMethod(
143		f = "univariate_ilr",
144		signature = c(object = "numeric"),
145		definition = function(object) {
146	12x	sqrt(1 / 2) * log(object / (1 - object))
147		}
148		)
149
150		#' @export
151		#' @rdname univariate_ilr
152		#' @aliases univariate_ilr,matrix-method
153		setMethod(
154		f = "univariate_ilr",
155		signature = c(object = "matrix"),
156		definition = function(object) {
157	2x	apply(X = object, MARGIN = 2, FUN = univariate_ilr)
158		}
159		)

1		# VARIANCE
2		#' @include AllGenerics.R
3		NULL
4
5		# Variance =====================================================================
6		#' @export
7		#' @rdname variance
8		#' @aliases variance,LogRatio-method
9		setMethod(
10		f = "variance",
11		signature = c("LogRatio"),
12		definition = function(x, row_weights = NULL, column_weights = TRUE) {
13	7x	n <- nrow(x)
14	7x	m <- ncol(x)
15
16	7x	w_row <- rep(1 / n, n)
17	7x	if (length(row_weights) == n) {
18	!	arkhe::assert_length(row_weights, n)
19	!	arkhe::assert_positive(row_weights, strict = FALSE)
20	!	w_row <- row_weights / sum(row_weights) # Sum up to 1
21		}
22
23	7x	w_col <- if (isTRUE(column_weights)) weights(x) else rep(1 / m, m)
24	7x	if (length(column_weights) == m) {
25	!	arkhe::assert_length(column_weights, m)
26	!	arkhe::assert_positive(column_weights, strict = FALSE)
27	!	w_col <- column_weights / sum(column_weights) # Sum up to 1
28		}
29
30	7x	z <- sweep(x, MARGIN = 2, STATS = colSums(x * w_row), FUN = "-")
31	7x	z <- colSums(diag(w_row) %% z^2 %% diag(w_col))
32	7x	names(z) <- colnames(x)
33	7x	z
34		}
35		)
36
37		# Total variance ===============================================================
38		#' @export
39		#' @describeIn variance_total The total variance of compositional data is the
40		#' trace of the [centred log-ratio covariance][covariance()] matrix
41		#' (i.e. totvar1 in Aitchison 1997).
42		#' @aliases variance_total,CompositionMatrix-method
43		setMethod(
44		f = "variance_total",
45		signature = c("CompositionMatrix"),
46		definition = function(x, sd = FALSE) {
47	3x	z <- sum(diag(covariance(x, center = TRUE)))
48	1x	if (sd) z <- sqrt((1 / (ncol(x) - 1)) * z)
49	3x	z
50		}
51		)
52
53		#' @export
54		#' @describeIn variance_total Computes the total log-ratio variance. This is
55		#' identical to the weighted sum-of-squared distances between samples
56		#' (i.e. totvar2 in Aitchison 1997).
57		#' @aliases variance_total,LogRatio-method
58		setMethod(
59		f = "variance_total",
60		signature = c("LogRatio"),
61		definition = function(x, row_weights = NULL, column_weights = TRUE) {
62	6x	sum(variance(x, row_weights = row_weights, column_weights = column_weights))
63		}
64		)
65

1		# HELPERS
2
3		## https://michaelchirico.github.io/potools/articles/developers.html
4		tr_ <- function(...) {
5	11x	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	46x	.row_names_info(x, type = 1L) > 0L &&
14	46x	!is.na(.row_names_info(x, type = 0L)[[1L]])
15		}
16
17		make_names <- function(x, n = length(x), prefix = "X") {
18	92x	x <- if (n > 0) x %\|\|% paste0(prefix, seq_len(n)) else character(0)
19	92x	x <- make.unique(x, sep = "_")
20	92x	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	10x	i <- !is.na(x)
57	10x	y <- x[i]
58	10x	y <- abs(y) * 100
59	10x	y <- format(y, trim = trim, digits = digits)
60	10x	y <- paste0(y, "%")
61	10x	x[i] <- y
62	10x	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		# 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		color = NULL, border = NA,
13		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 <- khroma::palette_color_discrete(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	mar <- graphics::par("mar")
34	4x	nlines <- height2line("M", cex = cex.axis)
35	4x	mar[1] <- 3
36	4x	mar[2] <- if (names) width2line(rownames(height), cex = cex.axis) + 0.1 else nlines
37	4x	mar[3] <- nlines
38	4x	mar[4] <- nlines
39
40	4x	old_par <- graphics::par(mar = mar)
41	4x	on.exit(graphics::par(old_par), add = TRUE)
42
43		## Open new window
44	4x	grDevices::dev.hold()
45	4x	on.exit(grDevices::dev.flush(), add = TRUE)
46	4x	graphics::plot.new()
47
48		## Set plotting coordinates
49	4x	space <- 1 / n * 0.5 - space / n * 0.5
50	4x	xlim <- range(0, 1)
51	4x	ylim <- range(xy$data$y)
52	4x	ylim <- ylim + c(0, 2 * (offset + space) * legend) + c(-1, 1) * space
53	4x	graphics::plot.window(xlim = xlim, ylim = ylim, xaxs = "i", yaxs = "i")
54
55		## Plot
56	4x	graphics::rect(
57	4x	xleft = xy$data$xmin,
58	4x	ybottom = xy$data$y - space,
59	4x	xright = xy$data$xmax,
60	4x	ytop = xy$data$y + space,
61	4x	col = col,
62	4x	border = border
63		)
64
65		## Construct axis
66	4x	if (axes) {
67	4x	at <- graphics::axTicks(side = 1)
68	4x	graphics::axis(side = 1, at = at, labels = label_percent(at),
69	4x	xpd = NA, las = 1, cex.axis = cex.axis, col.axis = col.axis,
70	4x	font.axis = font.axis)
71	4x	graphics::mtext(text = names(xy$groups), side = 4, line = 0, at = xy$groups,
72	4x	cex = cex.axis, col = col.axis, font = font.axis)
73	4x	if (names) {
74	4x	graphics::axis(side = 2, at = unique(xy$data$y), labels = unique(xy$data$row),
75	4x	las = 2, lty = 0, cex.axis = cex.axis, col.axis = col.axis,
76	4x	font.axis = font.axis)
77		}
78		}
79
80		## Add legend
81	4x	if (legend) {
82	4x	graphics::rect(
83	4x	xleft = cumsum(xy$mean) - xy$mean,
84	4x	ybottom = max(xy$data$y) + offset + space,
85	4x	xright = cumsum(xy$mean),
86	4x	ytop = max(xy$data$y) + 2 * (offset + space),
87	4x	col = unique(col),
88	4x	border = border
89		)
90	4x	lab <- label_chemical(names(xy$mean))
91	4x	graphics::mtext(text = parse(text = lab), side = 3, line = 0,
92	4x	at = cumsum(xy$mean) - xy$mean / 2,
93	4x	cex = cex.axis, col = unique(col), font = font.axis)
94		}
95
96	4x	invisible(height)
97		}
98
99		#' @export
100		#' @rdname barplot
101		#' @aliases barplot,CompositionMatrix-method
102		setMethod("barplot", c(height = "CompositionMatrix"), barplot.CompositionMatrix)
103
104		prepare_barplot <- function(x, order_rows = NULL, order_columns = FALSE,
105		decreasing = TRUE, offset = 0.025,
106		verbose = getOption("nexus.verbose")) {
107		## Remove missing values
108	4x	x <- arkhe::remove_NA(x, margin = 2, verbose = verbose)
109
110		## Relative frequencies
111	4x	n <- nrow(x)
112	4x	x <- x / rowSums(x)
113
114		## Validation
115	4x	stopifnot(is_composition(x))
116
117		## Row order
118	4x	if (!is.null(order_rows)) {
119	1x	j <- x[, order_rows, drop = TRUE]
120	1x	i <- order(j, decreasing = decreasing)
121	1x	x <- x[i, , drop = FALSE]
122		}
123
124		## Columns order
125	4x	center <- mean(x)
126	4x	if (order_columns) {
127	2x	col_order <- order(center, decreasing = FALSE)
128	2x	center <- center[col_order]
129	2x	x <- x[, col_order, drop = FALSE]
130		}
131
132		## Grouping
133	4x	if (!is_grouped(x)) {
134	3x	x <- group(x, by = rep(NA, n), verbose = FALSE)
135		}
136	4x	spl <- group_split(x)
137	4x	z <- do.call(rbind, spl)
138
139		## Build a long table
140	4x	row <- row(z, as.factor = TRUE)
141	4x	col <- col(z, as.factor = TRUE)
142	4x	data <- data.frame(
143	4x	row = as.vector(row),
144	4x	column = as.vector(col),
145	4x	value = as.vector(z)
146		)
147
148	4x	xmax <- t(apply(X = z, MARGIN = 1, FUN = cumsum))
149	4x	xmin <- xmax - z
150	4x	data$xmin <- as.vector(xmin)
151	4x	data$xmax <- as.vector(xmax)
152	4x	data$y <- as.vector(n + 1 - as.numeric(row)) / n # Reverse levels order
153
154		## Offset
155	4x	n_grp <- group_n(x)
156	4x	n_spl <- group_size(x)
157	4x	offset <- rev(seq_len(n_grp)) * offset - offset
158	4x	data$y <- data$y + rep(offset, n_spl)[as.numeric(row)]
159
160	4x	list(
161	4x	data = data,
162	4x	mean = center,
163	4x	groups = 1 - cumsum(n_spl) / n + n_spl / n * 0.5 + offset
164		)
165		}

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

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		# 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) [`factor`]s 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	22x	if (is.list(x)) {
52	2x	x <- if (length(x) > 1) interaction(x, sep = "_") else x[[1L]]
53		}
54	22x	x <- as.factor(x)
55	22x	if (drop_levels) x <- droplevels(x)
56	22x	if (allow_na) x <- addNA(x, ifany = TRUE)
57
58	22x	x
59		}
60
61		#' Validate Groups
62		#'
63		#' @param object A [`matrix`]-like object.
64		#' @param by A ([`list`] of) [`factor`]s for which interaction is to be computed.
65		#' @param verbose A [`logical`] scalar: should \R report extra information
66		#' on progress?
67		#' @return Invisibly returns `by`.
68		#' @keywords internal
69		#' @noRd
70		validate_groups <- function(object, by, verbose = getOption("nexus.verbose")) {
71	22x	arkhe::assert_type(by, "integer")
72	22x	arkhe::assert_length(by, nrow(object))
73
74	21x	if (nlevels(by) == 0) {
75	!	stop(tr_("Nothing to group by."), call. = FALSE)
76		}
77	21x	if (isTRUE(verbose)) {
78	1x	if (nlevels(by) == nrow(object)) {
79	!	message(tr_("As many groups as individuals."))
80		}
81
82	1x	n <- nlevels(by)
83	1x	what <- ngettext(n, "Found %g group (%s)", "Found %g groups (%s)")
84	1x	grp <- paste0(levels(by), collapse = ", ")
85	1x	message(sprintf(what, n, grp))
86		}
87
88	21x	invisible(by)
89		}
90
91		#' @export
92		#' @rdname group
93		#' @aliases group,CompositionMatrix-method
94		setMethod(
95		f = "group",
96		signature = "CompositionMatrix",
97		definition = function(object, by, verbose = getOption("nexus.verbose"), ...) {
98		## Compute groups
99	22x	by <- compute_groups(by, ...)
100
101		## Validation
102	22x	validate_groups(object, by, verbose = verbose)
103
104	21x	.GroupedComposition(
105	21x	object,
106	21x	group_indices = as.integer(by),
107	21x	group_levels = levels(by),
108	21x	group_ordered = is.ordered(by)
109		)
110		}
111		)
112
113		#' @export
114		#' @rdname group
115		#' @aliases group,GroupedComposition-method
116		setMethod(
117		f = "group",
118		signature = "GroupedComposition",
119		definition = function(object, by, add = FALSE,
120		verbose = getOption("nexus.verbose"), ...) {
121		## Compute groups
122	4x	if (isTRUE(add)) {
123	1x	if (!is.list(by)) by <- list(by)
124	1x	by <- c(list(group_factor(object, exclude = NULL)), by)
125		}
126	4x	methods::callNextMethod(object, by = by, verbose = verbose, ...)
127		}
128		)
129
130		#' @export
131		#' @rdname group
132		#' @aliases ungroup,GroupedComposition-method
133		setMethod(
134		f = "ungroup",
135		signature = "GroupedComposition",
136		definition = function(object) {
137	21x	methods::as(object, "CompositionMatrix", strict = TRUE)
138		}
139		)
140
141		#' @export
142		#' @rdname group
143		#' @aliases ungroup,GroupedLR-method
144		setMethod(
145		f = "ungroup",
146		signature = "GroupedLR",
147		definition = function(object) {
148	!	methods::as(object, "LR", strict = TRUE)
149		}
150		)
151
152		#' @export
153		#' @rdname group
154		#' @aliases ungroup,GroupedCLR-method
155		setMethod(
156		f = "ungroup",
157		signature = "GroupedCLR",
158		definition = function(object) {
159	!	methods::as(object, "CLR", strict = TRUE)
160		}
161		)
162
163		#' @export
164		#' @rdname group
165		#' @aliases ungroup,GroupedALR-method
166		setMethod(
167		f = "ungroup",
168		signature = "GroupedALR",
169		definition = function(object) {
170	!	methods::as(object, "ALR", strict = TRUE)
171		}
172		)
173
174		#' @export
175		#' @rdname group
176		#' @aliases ungroup,GroupedILR-method
177		setMethod(
178		f = "ungroup",
179		signature = "GroupedILR",
180		definition = function(object) {
181	1x	methods::as(object, "ILR", strict = TRUE)
182		}
183		)
184
185		#' @export
186		#' @rdname group
187		#' @aliases ungroup,GroupedPLR-method
188		setMethod(
189		f = "ungroup",
190		signature = "GroupedPLR",
191		definition = function(object) {
192	!	methods::as(object, "PLR", strict = TRUE)
193		}
194		)
195
196		# Metadata =====================================================================
197		#' @export
198		#' @rdname group_names
199		#' @aliases group_levels,ReferenceGroups-method
200		setMethod(
201		f = "group_levels",
202		signature = "ReferenceGroups",
203	234x	definition = function(object) object@group_levels
204		)
205
206		#' @export
207		#' @rdname group_names
208		#' @aliases group_names,ReferenceGroups-method
209		setMethod(
210		f = "group_names",
211		signature = "ReferenceGroups",
212	115x	definition = function(object) group_levels(object)[group_indices(object)]
213		)
214
215		is_ordered <- function(object) {
216	102x	isTRUE(object@group_ordered)
217		}
218
219		#' @export
220		#' @rdname group_names
221		#' @aliases group_factor,ReferenceGroups-method
222		setMethod(
223		f = "group_factor",
224		signature = "ReferenceGroups",
225		definition = function(object, exclude = NA) {
226	101x	factor(
227	101x	x = group_names(object),
228	101x	levels = group_levels(object),
229	101x	exclude = exclude,
230	101x	ordered = is_ordered(object)
231		)
232		}
233		)
234
235		#' @export
236		#' @rdname group_names
237		#' @aliases group_indices,ReferenceGroups-method
238		setMethod(
239		f = "group_indices",
240		signature = "ReferenceGroups",
241	117x	definition = function(object) object@group_indices
242		)
243
244		#' @export
245		#' @rdname group_names
246		#' @aliases group_rows,ReferenceGroups-method
247		setMethod(
248		f = "group_rows",
249		signature = "ReferenceGroups",
250		definition = function(object) {
251	18x	i <- group_factor(object, exclude = NULL)
252	18x	split(seq_along(i), f = i)
253		}
254		)
255
256		#' @export
257		#' @rdname group_names
258		#' @aliases group_n,ReferenceGroups-method
259		setMethod(
260		f = "group_n",
261		signature = "ReferenceGroups",
262	6x	definition = function(object) length(group_levels(object))
263		)
264
265		#' @export
266		#' @rdname group_names
267		#' @aliases group_size,ReferenceGroups-method
268		setMethod(
269		f = "group_size",
270		signature = "ReferenceGroups",
271	5x	definition = function(object) lengths(group_rows(object))
272		)
273
274		# Predicates ===================================================================
275		#' @export
276		#' @rdname is_assigned
277		#' @aliases is_assigned,ReferenceGroups-method
278		setMethod(
279		f = "is_assigned",
280		signature = "ReferenceGroups",
281	6x	definition = function(object) !is.na(group_names(object))
282		)
283
284		#' @export
285		#' @rdname is_assigned
286		#' @aliases any_assigned,ReferenceGroups-method
287		setMethod(
288		f = "any_assigned",
289		signature = "ReferenceGroups",
290	3x	definition = function(object) any(is_assigned(object))
291		)
292
293		#' @export
294		#' @rdname is_assigned
295		#' @aliases all_assigned,ReferenceGroups-method
296		setMethod(
297		f = "all_assigned",
298		signature = "ReferenceGroups",
299	1x	definition = function(object) all(is_assigned(object))
300		)

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	18x	from <- matrix(data = from, nrow = 1, ncol = length(from))
14	18x	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	27x	arkhe::assert_type(from, "numeric")
27
28		## Make row/column names
29	27x	lab <- make_names(x = NULL, n = nrow(from), prefix = "S")
30	27x	rownames(from) <- if (has_rownames(from)) rownames(from) else lab
31	27x	colnames(from) <- make_names(x = colnames(from), n = ncol(from), prefix = "V")
32
33		## Close
34	27x	totals <- rowSums(from, na.rm = TRUE)
35	27x	from <- from / totals
36
37	27x	.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	19x	lab <- make_names(x = NULL, n = nrow(from), prefix = "S")
51	19x	rownames(from) <- if (has_rownames(from)) rownames(from) else lab
52	19x	colnames(from) <- make_names(x = colnames(from), n = ncol(from), prefix = "V")
53
54		## Remove non-numeric columns
55	19x	if (is.null(parts)) {
56	18x	if (isTRUE(autodetect)) {
57	18x	parts <- arkhe::detect(from, f = is.numeric, margin = 2)
58	18x	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	19x	coda <- from[, parts, drop = FALSE]
72	19x	arkhe::assert_filled(coda)
73
74		## Build matrix
75	19x	coda <- data.matrix(coda, rownames.force = NA)
76	19x	totals <- rowSums(coda, na.rm = TRUE)
77	19x	coda <- coda / totals
78
79	19x	z <- .CompositionMatrix(coda, totals = unname(totals))
80	18x	if (is.null(groups)) return(z)
81
82		## Group names
83	1x	grp <- from[groups]
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, ..., group_var = ".group", group_after = 0) {
104	1x	z <- as.data.frame(methods::as(x, "matrix"), row.names = rownames(x), ...)
105
106	!	if (!is_grouped(x)) return(z)
107	1x	arkhe::append_column(z, group_names(x), after = group_after, var = group_var)
108		}
109
110		#' @export
111		#' @rdname as.data.frame
112		#' @aliases as.data.frame,CompositionMatrix-method
113		setMethod("as.data.frame", "CompositionMatrix", as.data.frame.CompositionMatrix)
114
115		#' @method as.data.frame LogRatio
116		#' @export
117		as.data.frame.LogRatio <- function(x, ..., group_var = ".group", group_after = 0) {
118	1x	z <- as.data.frame(methods::as(x, "matrix"), row.names = rownames(x), ...)
119
120		## Add columns labels
121	1x	lab <- labels(x)
122	1x	for (j in seq_len(ncol(z))) attr(z[[j]], "label") <- lab[[j]]
123
124	1x	if (!is_grouped(x)) return(z)
125	!	arkhe::append_column(z, group_names(x), after = group_after, var = group_var)
126		}
127
128		#' @export
129		#' @rdname as.data.frame
130		#' @aliases as.data.frame,LogRatio-method
131		setMethod("as.data.frame", "LogRatio", as.data.frame.LogRatio)
132
133		#' @method as.data.frame OutlierIndex
134		#' @export
135		as.data.frame.OutlierIndex <- function(x, row.names = rownames(x),
136		optional = FALSE, ...) {
137	!	as.data.frame(x@standard, row.names = row.names, optional = optional)
138		}
139
140		#' @export
141		#' @rdname as.data.frame
142		#' @aliases as.data.frame,OutlierIndex-method
143		setMethod("as.data.frame", "OutlierIndex", as.data.frame.OutlierIndex)

1		# SCALE
2		#' @include AllGenerics.R
3		NULL
4
5		# Scale ========================================================================
6		#' @export
7		#' @method scale CompositionMatrix
8		scale.CompositionMatrix <- function(x, center = TRUE, scale = TRUE) {
9	!	if (isFALSE(center) & isFALSE(scale)) return(x)
10
11	1x	y <- x
12	1x	if (!isFALSE(center)) {
13	1x	if (isTRUE(center)) center <- mean(x)
14	1x	arkhe::assert_type(center, "numeric")
15	1x	arkhe::assert_length(center, NCOL(x))
16
17	1x	y <- perturbation(y, 1 / center)
18		}
19
20	1x	if (!isFALSE(scale)) {
21	1x	if (isTRUE(scale)) scale <- sqrt(mean(diag(covariance(x, center = TRUE))))
22	1x	arkhe::assert_scalar(scale, "numeric")
23
24	1x	y <- powering(y, 1 / scale)
25		}
26
27	1x	y
28		}
29
30		#' @export
31		#' @rdname scale
32		#' @aliases scale,CompositionMatrix-method
33		setMethod("scale", "CompositionMatrix", scale.CompositionMatrix)

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	!	group_ordered = is_ordered(object))
48		}
49		)
50
51		#' @export
52		#' @rdname transform_clr
53		#' @aliases transform_clr,ALR-method
54		setMethod(
55		f = "transform_clr",
56		signature = c(object = "ALR"),
57		definition = function(object) {
58	3x	D <- ncol(object) + 1
59	3x	w <- object@weights
60
61	3x	base <- clr_base(D, weights = w)
62	3x	clr <- object %*% base[-D, ]
63	3x	dimnames(clr) <- list(rownames(object), object@parts)
64
65	3x	.CLR(
66	3x	clr,
67	3x	parts = object@parts,
68	3x	ratio = object@parts,
69	3x	order = seq_len(D),
70	3x	base = base,
71	3x	weights = w,
72	3x	totals = totals(object)
73		)
74		}
75		)
76
77		#' @export
78		#' @rdname transform_clr
79		#' @aliases transform_clr,GroupedALR-method
80		setMethod(
81		f = "transform_clr",
82		signature = c(object = "GroupedALR"),
83		definition = function(object) {
84	!	z <- methods::callNextMethod()
85	!	.GroupedCLR(z, group_indices = group_indices(object),
86	!	group_levels = group_levels(object),
87	!	group_ordered = is_ordered(object))
88		}
89		)

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	!	group_ordered = is_ordered(object))
34		}
35		)
36
37		## ALR -------------------------------------------------------------------------
38		#' @export
39		#' @rdname transform_inverse
40		#' @aliases transform_inverse,ALR,missing-method
41		setMethod(
42		f = "transform_inverse",
43		signature = c(object = "ALR", origin = "missing"),
44		definition = function(object) {
45	3x	y <- exp(object)
46	3x	y <- y / (1 + rowSums(y))
47	3x	z <- 1 - rowSums(y)
48
49	3x	y <- cbind(y, z)
50	3x	dimnames(y) <- list(rownames(object), object@parts)
51	3x	y <- y[, object@order]
52
53	3x	.CompositionMatrix(y, totals = totals(object))
54		}
55		)
56
57		#' @export
58		#' @rdname transform_inverse
59		#' @aliases transform_inverse,GroupedALR,missing-method
60		setMethod(
61		f = "transform_inverse",
62		signature = c(object = "GroupedALR", origin = "missing"),
63		definition = function(object) {
64	!	z <- methods::callNextMethod()
65	!	.GroupedComposition(z, group_indices = group_indices(object),
66	!	group_levels = group_levels(object),
67	!	group_ordered = is_ordered(object))
68		}
69		)
70
71		## ILR -------------------------------------------------------------------------
72		#' @export
73		#' @rdname transform_inverse
74		#' @aliases transform_inverse,ILR,missing-method
75		setMethod(
76		f = "transform_inverse",
77		signature = c(object = "ILR", origin = "missing"),
78		definition = function(object) {
79	2x	y <- tcrossprod(object@.Data, object@base)
80	2x	y <- exp(y)
81	2x	y <- y / rowSums(y)
82
83	2x	dimnames(y) <- list(rownames(object), object@parts)
84	2x	y <- y[, object@order]
85
86	2x	.CompositionMatrix(y, totals = totals(object))
87		}
88		)
89
90		#' @export
91		#' @rdname transform_inverse
92		#' @aliases transform_inverse,GroupedILR,missing-method
93		setMethod(
94		f = "transform_inverse",
95		signature = c(object = "GroupedILR", origin = "missing"),
96		definition = function(object) {
97	!	z <- methods::callNextMethod()
98	!	.GroupedComposition(z, group_indices = group_indices(object),
99	!	group_levels = group_levels(object),
100	!	group_ordered = is_ordered(object))
101		}
102		)
103
104		#' @export
105		#' @rdname transform_inverse
106		#' @aliases transform_inverse,matrix,ILR-method
107		setMethod(
108		f = "transform_inverse",
109		signature = c(object = "matrix", origin = "ILR"),
110		definition = function(object, origin) {
111	1x	y <- tcrossprod(object, origin@base)
112	1x	y <- exp(y)
113	1x	y <- y / rowSums(y)
114
115	1x	dimnames(y) <- list(rownames(object), origin@parts)
116	1x	y <- y[, origin@order]
117
118	1x	y
119		}
120		)

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	!	group_ordered = is_ordered(object))
67		}
68		)

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		# 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		# 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(as.matrix(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		# 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, verbose = verbose)
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, exclude = NULL)
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(
41	5x	X = totals(x),
42	5x	INDEX = group_factor(x, exclude = NULL),
43	5x	FUN = mean
44		)
45
46	5x	z <- .CompositionMatrix(z, totals = as.numeric(tot))
47	5x	group(z, by = flatten_chr(x = grp, by = group_factor(x, exclude = NULL)),
48	5x	verbose = verbose)
49		}
50		)
51
52		flatten_chr <- function(x, by) {
53	5x	x <- as.character(x)
54	5x	z <- tapply(X = x, INDEX = by, FUN = unique, simplify = FALSE)
55	5x	z <- vapply(X = z, FUN = paste0, FUN.VALUE = character(1), collapse = ":")
56	5x	z
57		}

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		# 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	!	group_ordered = is_ordered(object))
61		}
62		)
63
64		#' @export
65		#' @rdname transform_alr
66		#' @aliases transform_alr,CLR-method
67		setMethod(
68		f = "transform_alr",
69		signature = c(object = "CLR"),
70		definition = function(object, j = ncol(object)) {
71	2x	D <- ncol(object)
72	2x	parts <- object@parts
73
74		## Reorder (move denominator)
75	2x	j <- if (is.character(j)) which(parts == j) else as.integer(j)
76	2x	ordering <- c(which(j != seq_len(D)), j)
77	2x	parts <- parts[ordering]
78	2x	z <- object[, ordering, drop = FALSE]
79
80		## Compute ratios
81	2x	base <- alr_base(D)
82	2x	alr <- z %*% base
83	2x	rownames(alr) <- rownames(object)
84	2x	colnames(alr) <- paste(parts[-D], parts[D], sep = "_")
85
86	2x	.ALR(
87	2x	alr,
88	2x	parts = parts,
89	2x	ratio = paste(parts[-D], parts[D], sep = "/"),
90	2x	order = order(ordering),
91	2x	base = base,
92	2x	weights = object@weights,
93	2x	totals = totals(object)
94		)
95		}
96		)
97
98		#' @export
99		#' @rdname transform_alr
100		#' @aliases transform_alr,GroupedCLR-method
101		setMethod(
102		f = "transform_alr",
103		signature = c(object = "GroupedCLR"),
104		definition = function(object, j = ncol(object), weights = FALSE) {
105	!	z <- methods::callNextMethod()
106	!	.GroupedALR(z, group_indices = group_indices(object),
107	!	group_levels = group_levels(object),
108	!	group_ordered = is_ordered(object))
109		}
110		)

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