isopleuros coverage - 86.28%

Files
Source

# TERNARY TITLE
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_title
ternary_title <- function(main = NULL, sub = NULL, xlab = NULL, ylab = NULL,
                          zlab = NULL, line = NA, outer = FALSE, ...) {
  ## 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")

  ## Axes labels
  xlab <- grDevices::as.graphicsAnnot(xlab)
  ylab <- grDevices::as.graphicsAnnot(ylab)
  zlab <- grDevices::as.graphicsAnnot(zlab)

  if (!is.null(xlab)) {
    graphics::text(x = 0, y = 0, label = xlab, pos = 1,
                   col = col.lab, cex = cex.lab, font = font.lab)
  }
  if (!is.null(ylab)) {
    graphics::text(x = 1, y = 0, label = ylab, pos = 1,
                   col = col.lab, cex = cex.lab, font = font.lab)
  }
  if (!is.null(zlab)) {
    graphics::text(x = 0.5, y = the$top, label = zlab, pos = 3,
                   col = col.lab, cex = cex.lab, font = font.lab)
  }

  ## Title
  graphics::title(main = main, sub = sub, xlab = NULL, ylab = NULL,
                  line = line, outer = outer)

  invisible(NULL)
}

# CONTOUR
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_contour
#' @aliases ternary_contour,numeric,numeric,numeric-method
setMethod(
  f = "ternary_contour",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, value, n = 50, nlevels = 10,
                        levels = pretty(range(value, na.rm = TRUE), nlevels),
                        ilr = TRUE, method = "linear", extrapolate = FALSE,
                        palette = function(i) grDevices::hcl.colors(i, "YlOrRd", rev = TRUE),
                        ...) {
    ## Calculate contour lines
    xy <- coordinates_contour(x = x, y = y, z = z, value = value, n = n,
                              nlevels = nlevels, levels = levels,
                              ilr = ilr, method = method,
                              extrapolate = extrapolate)

    ## Get contour levels
    lvl <- vapply(X = xy, FUN = getElement, FUN.VALUE = numeric(1),
                  name = "level")

    ## Colors
    ## (number of levels may differ from nlevels due to interp())
    col <- palette(length(unique(lvl)))
    names(col) <- unique(lvl)
    col <- col[as.character(lvl)]

    ## Plot
    for (i in seq_along(xy)) {
      ## Get contour
      level <- xy[[i]]

      ## Inverse ILR transform
      tern <- cbind(level$x, level$y)
      tern <- if (ilr) ilr_inv(tern) else coordinates_cartesian(tern)

      ## Plot ternary lines
      ternary_lines(tern, col = col[[i]], ...)
    }

    invisible(list(levels = lvl, colors = col))
  }
)

#' @export
#' @rdname ternary_contour
#' @aliases ternary_contour,ANY,missing,missing-method
setMethod(
  f = "ternary_contour",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, value, n = 50, nlevels = 10,
                        levels = pretty(range(value, na.rm = TRUE), nlevels),
                        ilr = TRUE, method = "linear", extrapolate = FALSE,
                        palette = function(i) grDevices::hcl.colors(i, "YlOrRd", rev = TRUE),
                        ...) {
    xyz <- grDevices::xyz.coords(x)
    coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z, value = value,
                                   n = n, nlevels = nlevels, levels = levels,
                                   ilr = ilr, method = method,
                                   extrapolate = extrapolate,
                                   palette = palette, ...)
    invisible(coords)
  }
)

#' Calculate Contour Lines
#'
#' Computes contours coordinates.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param value A [`numeric`] [`matrix`] containing the values to be plotted.
#' @param n A length-one [`numeric`] specifying the number of grid points.
#' @param nlevels A length-one [`numeric`] vector specifying the number of
#'  contour levels desired. Only used if `levels` is `NULL`.
#' @param levels A [`numeric`] vector of levels at which to draw contour lines.
#' @param ilr A [`logical`] scalar: should interpolation be computed in ILR
#'  space? If `FALSE`, interpolation is computed in Cartesian space.
#' @param method A [`character`] string: specifying the method for interpolation
#'  (see [interp::interp()]).
#' @param extrapolate A [`logical`] scalar: should extrapolation be used outside
#'  of the convex hull determined by the data points (see [interp::interp()])?
#' @param ... Further parameters to be passed to [interp::interp()].
#' @return
#'  A [`list`] of contours, each itself a list with elements
#'  (see [grDevices::contourLines()]):
#'  \tabular{ll}{
#'   `level` \tab The contour level. \cr
#'   `x`     \tab The (ILR) x-coordinates of the contour. \cr
#'   `y`     \tab The (ILR) y-coordinates of the contour. \cr
#'  }
#' @keywords internal
#' @noRd
coordinates_contour <- function(x, y, z, value, n = 50, nlevels = 10,
                                levels = pretty(range(value, na.rm = TRUE), nlevels),
                                ilr = TRUE, method = "linear", extrapolate = FALSE,
                                ...) {
  ## Validation
  assert_package("interp")
  assert_length(value, length(x))

  ## ILR vs Cartesian
  coda <- cbind(x, y, z)
  ratio <- if (ilr) ilr(coda) else do.call(cbind, coordinates_ternary(coda))

  ## Remove NA/Inf (if any)
  ok <- apply(X = ratio, MARGIN = 1, FUN = function(x) all(is.finite(x)))
  ratio <- ratio[ok, , drop = FALSE]
  value <- value[ok]

  ## Interpolate
  xlim <- expand_range(ratio[, 1], mult = 0.2)
  ylim <- expand_range(ratio[, 2], mult = 0.2)
  interp <- interp::interp(
    x = ratio[, 1],
    y = ratio[, 2],
    z = value,
    xo = seq(xlim[1L], xlim[2L], length.out = n),
    yo = seq(ylim[1L], ylim[2L], length.out = n),
    method = method,
    extrap = extrapolate,
    ...
  )

  ## Compute contours
  grDevices::contourLines(
    x = interp$x,
    y = interp$y,
    z = interp$z,
    nlevels = nlevels,
    levels = levels
  )
}

# TERNARY COORDINATES
#' @include AllGenerics.R
NULL

# Ternary to cartesian =========================================================
#' @export
#' @rdname coordinates_ternary
#' @aliases coordinates_ternary,numeric,numeric,numeric-method
setMethod(
  f = "coordinates_ternary",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, center = FALSE, scale = FALSE,
                        missing = getOption("isopleuros.missing")) {
    ## Validation
    n <- length(x)
    assert_length(y, n)
    assert_length(z, n)
    assert_center(center)
    assert_scale(scale)

    ## Missing values
    if (missing) {
      x[is.na(x)] <- 0
      y[is.na(y)] <- 0
      z[is.na(z)] <- 0
    }

    total <- x + y + z
    na <- is.na(x) | is.na(y) | is.na(z)
    zero <- total == 0

    x <- x[!na & !zero]
    y <- y[!na & !zero]
    z <- z[!na & !zero]
    total <- total[!na & !zero]

    ## Validation
    if (any(x < 0 | y < 0 | z < 0)) {
      stop(tr_("Positive values are expected."), call. = FALSE)
    }

    coord <- matrix(data = c(x, y, z), ncol = 3) / total
    coord <- scale(coord, center = center, scale = scale)

    list(
      x = coord$y + coord$z / 2,
      y = coord$z * sqrt(3) / 2,
      center = coord$center,
      scale = coord$scale
    )
  }
)

#' @export
#' @rdname coordinates_ternary
#' @aliases coordinates_ternary,ANY,missing,missing-method
setMethod(
  f = "coordinates_ternary",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, xlab = NULL, ylab = NULL, zlab = NULL,
                        center = FALSE, scale = FALSE,
                        missing = getOption("isopleuros.missing")) {
    xyz <- grDevices::xyz.coords(x, xlab = xlab, ylab = ylab, zlab = zlab)
    methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
                         center = center, scale = scale, missing = missing)
  }
)

# Cartesian to ternary =========================================================
#' @export
#' @rdname coordinates_cartesian
#' @aliases coordinates_cartesian,numeric,numeric-method
setMethod(
  f = "coordinates_cartesian",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    ## Validation
    n <- length(x)
    assert_length(y, n)

    k <- y * 2 / sqrt(3)
    j <- x - k / 2
    i <- 1 - (j + k)

    list(
      x = i,
      y = j,
      z = k
    )
  }
)

#' @export
#' @rdname coordinates_cartesian
#' @aliases coordinates_cartesian,ANY,missing-method
setMethod(
  f = "coordinates_cartesian",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, xlab = NULL, ylab = NULL) {
    xy <- grDevices::xy.coords(x, xlab = xlab, ylab = ylab)
    methods::callGeneric(x = xy$x, y = xy$y)
  }
)

# Scale ========================================================================
#' Center and Scale
#'
#' @param x,y,z The x, y and z coordinates of a set of points. Both y and z can
#'  be left at `NULL`. In this case, an attempt is made to interpret x in a way
#'  suitable for plotting.
#' @param center A [`logical`] scalar or a [`numeric`] vector giving the center.
#' @param scale A [`logical`] scalar or a length-one [`numeric`] vector giving a
#'  scaling factor.
#' @return
#'  A [`list`] with the components:
#'  \tabular{ll}{
#'   `x` \tab A [`numeric`] vector of x values. \cr
#'   `y` \tab A [`numeric`] vector of y values. \cr
#'   `z` \tab A [`numeric`] vector of z values. \cr
#'   `center` \tab A [`numeric`] vector giving the center. \cr
#'   `scale` \tab A [`numeric`] vector giving the scale factor. \cr
#'  }
#' @keywords internal
#' @noRd
scale <- function(x, y = NULL, z = NULL, center = TRUE, scale = TRUE) {
  xyz <- grDevices::xyz.coords(x = x, y = y, z = z)
  xyz <- matrix(data = c(xyz$x, xyz$y, xyz$z), ncol = 3)

  y <- xyz
  if (!isFALSE(center) && !is.null(center)) {
    if (isTRUE(center)) {
      center <- apply(X = xyz, MARGIN = 2, FUN = gmean, simplify = TRUE)
      center <- center / sum(center)
    }
    assert_length(center, NCOL(xyz))

    y <- t(t(y) / center)
    y <- y / rowSums(y)
  } else {
    center <- rep(1, NCOL(xyz))
  }

  if (!isFALSE(scale) && !is.null(scale)) {
    if (isTRUE(scale)) {
      scale <- sqrt(mean(diag(stats::cov(clr(xyz)))))
    }
    assert_length(scale, 1)

    y <- y^(1 / scale)
    y <- y / rowSums(y)
  } else {
    scale <- 1
  }

  list(
    x = y[, 1],
    y = y[, 2],
    z = y[, 3],
    center = center,
    scale = scale
  )
}

#' 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 na.rm A [`logical`] scalar: should `NA` values be stripped before the
#'  computation proceeds?
#' @return A [`numeric`] vector.
#' @keywords internal
#' @noRd
gmean <- function(x, trim = 0, na.rm = FALSE) {
  if (na.rm) x <- x[is.finite(x)]
  x <- x[x > 0]
  exp(mean(log(x), trim = trim))
}

# Centered Log-Ratios ==========================================================
#' Centered Log-Ratios (CLR)
#'
#' Computes CLR transformation.
#' @param x A [`numeric`] `matrix`.
#' @keywords internal
#' @noRd
clr <- function(x) {
  J <- ncol(x)
  clr <- log(x, base = exp(1)) %*% diag(J)

  clr
}

#' Inverse Centered Log-Ratios Transformation
#'
#' Computes inverse CLR transformation.
#' @param x A [`numeric`] `matrix` of log ratios.
#' @keywords internal
#' @noRd
clr_inv <- function(x) {
  y <- exp(x)
  y <- y / rowSums(y)

  y
}

# Isometric Log-Ratios =========================================================
#' Isometric Log-Ratios (ILR)
#'
#' Computes ILR transformation.
#' @param x A [`numeric`] `matrix`.
#' @keywords internal
#' @noRd
ilr <- function(x) {
  D <- ncol(x)
  H <- ilr_base(D)

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

  ilr
}

#' Canonical Basis for Isometric Log-Ratio transformation
#'
#' Computes the canonical basis in the CLR plane used for ILR transformation.
#' @param a A [`numeric`] value giving the number of parts of the simplex.
#' @keywords internal
#' @noRd
ilr_base <- function(n) {
  seq_parts <- seq_len(n - 1)

  ## Helmert matrix (rotation matrix)
  H <- stats::contr.helmert(n)                  # n x n-1
  H <- t(H) / sqrt((seq_parts + 1) * seq_parts) # n-1 x n

  ## Center
  m <- diag(x = 1, nrow = n) - matrix(data = 1 / n, nrow = n, ncol = n)
  H <- tcrossprod(m, H)

  H
}

#' Inverse Isometric Log-Ratio Transformation
#'
#' Computes inverse ILR transformation.
#' @param x A [`numeric`] `matrix` of log ratios.
#' @keywords internal
#' @noRd
ilr_inv <- function(x) {
  D <- ncol(x) + 1
  H <- ilr_base(D)

  y <- tcrossprod(x, H)
  y <- exp(y)
  y <- y / rowSums(y)

  y
}

# DENSITY
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_density
#' @aliases ternary_density,numeric,numeric,numeric-method
setMethod(
  f = "ternary_density",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, h = NULL, n = 25, nlevels = 10, levels = NULL,
                        palette = function(i) grDevices::hcl.colors(i, "YlOrRd", rev = TRUE),
                        ...) {
    ## Calculate density contour lines
    xy <- coordinates_kde(x = x, y = y, z = z, h = h, n = n,
                          nlevels = nlevels, levels = levels)

    ## Get contour levels
    lvl <- vapply(X = xy, FUN = getElement, FUN.VALUE = numeric(1),
                  name = "level")

    ## Colors
    ## (number of levels may differ from nlevels due to pretty())
    col <- palette(length(unique(lvl)))
    names(col) <- unique(lvl)
    col <- col[as.character(lvl)]

    ## Plot
    for (i in seq_along(xy)) {
      ## Get contour
      level <- xy[[i]]

      ## Inverse ILR transform
      tern <- ilr_inv(cbind(level$x, level$y))

      ## Plot ternary lines
      ternary_lines(tern, col = col[[i]], ...)
    }

    invisible(list(levels = lvl, colors = col))
  }
)

#' @export
#' @rdname ternary_density
#' @aliases ternary_density,ANY,missing,missing-method
setMethod(
  f = "ternary_density",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, h = NULL, n = 25, nlevels = 10, levels = NULL,
                        palette = function(i) grDevices::hcl.colors(i, "YlOrRd", rev = TRUE),
                        ...) {
    xyz <- grDevices::xyz.coords(x)
    pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
                               h = h, n = n, nlevels = nlevels, levels = levels,
                               palette = palette, ...)
    invisible(pt)
  }
)

#' Calculate KDE Contour Lines
#'
#' Computes 2D-KDE contours coordinates.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param h A length-one [`numeric`] vector giving the bandwidth.
#' @param n A length-one [`numeric`] specifying the number of grid points.
#' @param nlevels A length-one [`numeric`] vector specifying the number of
#'  contour levels desired. Only used if `levels` is `NULL`.
#' @param levels A [`numeric`] vector of levels at which to draw contour lines.
#' @return
#'  A [`list`] of contours, each itself a list with elements
#'  (see [grDevices::contourLines()]):
#'  \tabular{ll}{
#'   `level` \tab The contour level. \cr
#'   `x`     \tab The ILR x-coordinates of the contour. \cr
#'   `y`     \tab The ILR y-coordinates of the contour. \cr
#'  }
#' @keywords internal
#' @noRd
coordinates_kde <- function(x, y, z, h = NULL, n = 25,
                            nlevels = 10, levels = NULL) {
  ## ILR
  coda <- cbind(x, y, z)
  ratio <- ilr(coda)

  ## Compute KDE
  lims <- expand_range(ratio, mult = 0.2)
  dens <- kde(
    x = ratio[, 1],
    y = ratio[, 2],
    h = h,
    n = n,
    xlim = lims, # x and y range should be same
    ylim = lims
  )

  ## Compute contours
  grDevices::contourLines(
    x = dens$x,
    y = dens$y,
    z = dens$z,
    nlevels = nlevels,
    levels = levels %||% pretty(range(dens$z, na.rm = TRUE), nlevels)
  )
}

## Adapted from MASS::kde2d
kde <- function(x, y, h = NULL, n = 25, gx = NULL, gy = NULL,
                xlim = range(x), ylim = range(y)) {
  n <- rep(n, length.out = 2L)
  if (is.null(gx)) gx <- seq(xlim[1L], xlim[2L], length.out = n[1L])
  if (is.null(gy)) gy <- seq(ylim[1L], ylim[2L], length.out = n[2L])

  if (is.null(h)) {
    h <- c(bandwidth(x), bandwidth(y))
  } else {
    h <- rep(h, length.out = 2L)
  }
  h <- h / 4

  if (any(h <= 0))
    stop("Bandwidths must be strictly positive.", call. = FALSE)

  ax <- outer(gx, x, "-") / h[1L]
  ay <- outer(gy, y, "-") / h[2L]

  nx <- length(x)
  z <- tcrossprod(
    matrix(stats::dnorm(ax), ncol = nx),
    matrix(stats::dnorm(ay), ncol = nx)
  )
  z <- z / (nx * h[1L] * h[2L])

  list(x = gx, y = gy, z = z)
}

# Copied from MASS::bandwidth.nrd()
bandwidth <- function(x) {
  r <- stats::quantile(x, c(0.25, 0.75))
  h <- (r[2L] - r[1L]) / 1.34
  4 * 1.06 * min(sqrt(stats::var(x)), h) * length(x)^(-1 / 5)
}

# TERNARY LABELS
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_labels
#' @aliases ternary_labels,numeric,numeric,numeric-method
setMethod(
  f = "ternary_labels",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, center = FALSE, scale = FALSE,
                        labels = seq_along(x), type = c("text", "shadow"), ...) {
    ## Validation
    type <- match.arg(type, several.ok = FALSE)

    ## Compute label positions
    coords <- coordinates_ternary(x, y, z, center = center, scale = scale)
    labs <- compute_labels(x = coords$x, y = coords$y, labels = labels, ...)

    ## Draw labels
    fun <- switch(
      type,
      text = graphics::text,
      shadow = text_shadow
    )
    fun(labs, labels = labels, ...)

    coords <- utils::modifyList(coords, list(x = x, y = y, z = z))
    invisible(coords)
  }
)

#' @export
#' @rdname ternary_labels
#' @aliases ternary_labels,ANY,missing,missing-method
setMethod(
  f = "ternary_labels",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, center = FALSE, scale = FALSE, labels = seq_along(x$x), ...) {
    x <- grDevices::xyz.coords(x)
    coords <- methods::callGeneric(x = x$x, y = x$y, z = x$z,
                                   center = center, scale = scale,
                                   labels = labels, ...)
    invisible(coords)
  }
)


#' Compute Label Positions
#'
#' @return A [`list`] with elements `x`, `y` and `labels`.
#' @source
#'  This function is modeled after [car::pointLabel()] (originally from the
#'  \pkg{maptools} package).
#' @keywords internal
#' @noRd
compute_labels <- function(x, y, labels, ..., iter = 50,
                           cex = graphics::par("cex"),
                           font = NULL, vfont = NULL) {
  ## Coordinates
  bound <- graphics::par("usr")
  ratio <- graphics::par("pin")[1] / graphics::par("pin")[2] # x/y ratio

  to_unity <- function(x, y) {
    list(x = (x - bound[1]) / (bound[2] - bound[1]) * ratio,
         y = (y - bound[3]) / (bound[4] - bound[3]) / ratio)
  }
  to_usr <- function(x, y) {
    list(x = bound[1] + x / ratio * (bound[2] - bound[1]),
         y = bound[3] + y * ratio * (bound[4] - bound[3]))
  }

  xy <- to_unity(x = x, y = y)
  x <- xy$x
  y <- xy$y
  n <- length(x)

  ## 8 positions: corners and side mid-points of the rectangle
  ## Position 7 (top right) is the most preferred
  width <- graphics::strwidth(labels, units = "figure", cex = cex,
                              font = font, vfont = vfont)
  height <- graphics::strheight(labels, units = "figure", cex = cex,
                                font = font, vfont = vfont)
  width <- (width + 0.02) * ratio
  height <- (height + 0.02) / ratio

  makeoff <- function(pos) {
    c(-1, -1, -1, 0, 0, 1, 1, 1)[pos] * (width / 2) +
      1i * c(-1, 0, 1, -1, 1, -1, 0, 1)[pos] * (height / 2)
  }

  ## Find intersection area of two rectangles
  overlap <- function(xy1, off1, xy2, off2) {
    w <- pmin(Re(xy1 + off1 / 2), Re(xy2 + off2 / 2)) -
      pmax(Re(xy1 - off1 / 2), Re(xy2 - off2 / 2))
    h <- pmin(Im(xy1 + off1 / 2), Im(xy2 + off2 / 2)) -
      pmax(Im(xy1 - off1 / 2), Im(xy2 - off2 / 2))
    w[w <= 0] <- 0
    h[h <= 0] <- 0
    w * h
  }

  objective <- function(gene) {
    offset <- makeoff(gene)

    if (!is.null(rectidx1)) {
      area <- sum(overlap(xy[rectidx1] + offset[rectidx1], rectv[rectidx1],
                          xy[rectidx2] + offset[rectidx2], rectv[rectidx2]))
    } else {
      area <- 0
    }

    ## Penalize labels which go outside the image area
    ## Count points outside of the image
    a <- Re(xy + offset - rectv / 2) < 0 | Re(xy + offset + rectv / 2) > ratio
    b <- Im(xy + offset - rectv / 2) < 0 | Im(xy + offset + rectv / 2) > 1 / ratio
    outside <- sum(a | b)
    res <- 1000 * area + outside
    res
  }

  # Make a list of label rectangles in their reference positions,
  # centered over the map feature; the real labels are displaced
  # from these positions so as not to overlap
  # Note that some labels can be bigger than others
  xy <- x + 1i * y
  rectv <- width + 1i * height

  rectidx1 <- rectidx2 <- array(0, (length(x)^2 - length(x)) / 2)
  k <- 0
  for (i in seq_along(x))
    for (j in seq_len(i - 1)) {
      k <- k + 1
      rectidx1[k] <- i
      rectidx2[k] <- j
    }
  maylap <- overlap(xy[rectidx1], 2 * rectv[rectidx1],
                    xy[rectidx2], 2 * rectv[rectidx2]) > 0
  rectidx1 <- rectidx1[maylap]
  rectidx2 <- rectidx2[maylap]

  ## Simulated annealing
  ## Initial state
  gene <- rep(8, n)
  score <- objective(gene)
  ## Initial "best" solution
  bestgene <- gene
  bestscore <- score
  iter <- seq_len(iter)
  temp <- 2.5
  for (i in iter) {
    k <- 1 # Energy evaluation count
    for (j in iter) {
      newgene <- gene
      newgene[sample(n, 1)] <- sample(8, 1)
      newscore <- objective(newgene)
      if (newscore <= score || stats::runif(1) < exp((score - newscore) / temp)) {
        ## keep the new set if it has the same or better score or
        ## if it's worse randomly based on the annealing criteria
        k <- k + 1
        score <- newscore
        gene <- newgene
      }
      if (score <= bestscore) {
        bestscore <- score
        bestgene <- gene
      }
      if (bestscore == 0 || k == 10) break
    }
    if (bestscore == 0) break
    temp <- 0.9 * temp
  }

  nx <- Re(xy + makeoff(bestgene))
  ny <- Im(xy + makeoff(bestgene))

  xy <- to_usr(x = nx, y = ny)
  xy$labels <- labels
  xy
}

#' Shadow Text
#'
#' @param x,y A [`numeric`] vector. If `y` is `NULL`, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param labels A [`character`] vector specifying the text to be written.
#' @param width Thickness of the shadow, as a fraction of the plotting size.
#' @param theta Angles for plotting the background.
#' @param cex A [`numeric`] character expansion factor.
#' @param col The color to be used for the text.
#' @param bg The color to be used for the shadow.
#' @param font,vfont The font to be used (see [graphics::text()]).
#' @param ... Further parameters to be passed to [graphics::text()].
#' @return
#'  `text_shadow()` is called it for its side-effects: it results in a graphic
#'  being displayed.
#' @author N. Frerebeau
#' @keywords internal
#' @noRd
text_shadow <- function(x, y = NULL, labels = seq_along(x$x),
                        width = 1/10, theta = seq(0, 2 * pi, length.out = 50),
                        cex = graphics::par("cex"), col = graphics::par("fg"),
                        bg = graphics::par("bg"), font = NULL, vfont = NULL, ...) {

  x <- grDevices::xy.coords(x = x, y = y)

  xo <- width * graphics::strwidth("M", units = "user", cex = cex, font = font, vfont = vfont)
  yo <- width * graphics::strheight("X", units = "user", cex = cex, font = font, vfont = vfont)

  for (i in theta) {
    graphics::text(x = x$x + cos(i) * xo, y = x$y + sin(i) * yo, labels = labels,
                   col = bg, cex = cex, font = font, vfont = vfont, ...)
  }

  graphics::text(x = x$x, y = x$y, labels = labels, col = col, cex = cex,
                 font = font, vfont = vfont, ...)

  invisible(NULL)
}

# MEAN
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_mean
#' @aliases ternary_mean,numeric,numeric,numeric-method
setMethod(
  f = "ternary_mean",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, ...) {
    x <- gmean(x)
    y <- gmean(y)
    z <- gmean(z)

    pt <- ternary_points(x = x, y = y, z = z, ...)
    invisible(pt)
  }
)

#' @export
#' @rdname ternary_mean
#' @aliases ternary_mean,ANY,missing,missing-method
setMethod(
  f = "ternary_mean",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, ...) {
    xyz <- grDevices::xyz.coords(x)
    pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z, ...)
    invisible(pt)
  }
)

#' 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 na.rm A [`logical`] scalar: should `NA` values be stripped before the
#'  computation proceeds?
#' @return A [`numeric`] vector.
#' @keywords internal
#' @noRd
gmean <- function(x, trim = 0, na.rm = FALSE) {
  index <- is.finite(x) & x > 0
  exp(mean(log(unclass(x)[index]), trim = trim, na.rm = na.rm))
}

# HELPERS

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

`%||%` <- function(x, y) {
  if (!is.null(x)) x else y
}

map_color <- function(values, palette = NULL, scheme = "viridis",
                      ignore_zero = FALSE) {
  if (isFALSE(palette)) return(values)

  if (is.null(palette)) {
    palette <- function(x) {
      x <- (x - min(x)) / (max(x) - min(x)) # Rescale to [0,1]
      col <- grDevices::hcl.colors(256L, palette = scheme)
      grDevices::rgb(grDevices::colorRamp(col)(x), maxColorValue = 255)
    }
  }

  color <- rep(NA, length(values))
  ok <- is.finite(values) # Remove NA/Inf (if any)
  if (ignore_zero) ok[ok] <- values[ok] > 0
  color[ok] <- palette(values[ok])
  color
}

#' Expand Range
#'
#' @param x A [`numeric`] vector.
#' @param mult A [`numeric`] value giving the multiplicative constant.
#' @param add A [`numeric`] value giving the additive constant.
#' @return A length-two [`numeric`] vector.
#' @keywords internal
#' @noRd
expand_range <- function(x, mult = 0, add = 0) {
  lims <- range(x)
  lims <- lims + c(-1, 1) * (diff(lims) * mult + add)
  lims
}

#' Rotate Around a Point
#'
#' @param x A column vector giving the x and y coordinates of the point to be
#'  rotated.
#' @param theta A length-one [`numeric`] vector specifying the rotation angle
#'  (in radian).
#' @param origin A length-two [`numeric`] vector specifying the coordinates
#'  of the point to rotate around.
#' @return A `matrix` of coordinates.
#' @keywords internal
#' @noRd
rotate <- function(x, theta = 0, origin = c(0.5, sqrt(3) / 6)) {
  ## Translation matrix
  trans <- diag(1, 3, 3)
  trans[, 3] <- c(origin, 1)

  ## Rotation matrix
  rot <- matrix(
    data = c(cos(theta), sin(theta), 0, -sin(theta), cos(theta), 0, 0, 0, 1),
    nrow = 3,
    ncol = 3
  )

  x <- as.matrix(x)
  if (nrow(x) < 3) x <- rbind(x, rep(1, ncol(x)))
  t(trans %*% rot %*% solve(trans) %*% x)
}

#' Check Object Length
#'
#' @param x An object to be checked.
#' @param expected An appropriate expected value.
#' @return
#'  Throws an error, if any, and returns `x` invisibly otherwise.
#' @keywords internal
#' @noRd
assert_length <- function(x, expected) {
  arg <- deparse(substitute(x))
  if (length(x) != expected) {
    str <- tr_("%s must be of length %d; not %d.")
    msg <- sprintf(str, sQuote(arg), expected, length(x))
    stop(msg, call. = FALSE)
  }
  invisible(x)
}

assert_center <- function(x, current = get_center()) {
  centered <- is.numeric(current) && !all(current == 1)
  ok <- isTRUE(x) || is.numeric(x)
  if (!ok && centered) {
    msg <- tr_("The current plot has been centered, but your data doesn't seem to be.")
    message(msg)
  }
  invisible(x)
}

assert_scale <- function(x, current = get_scale()) {
  scaled <- is.numeric(current) && !all(current == 1)
  ok <- isTRUE(x) || is.numeric(x)
  if (!ok && scaled) {
    msg <- tr_("The current plot has been scaled, but your data doesn't seem to be.")
    message(msg)
  }
  invisible(x)
}

assert_package <- function(x) {
  if (!requireNamespace(x, quietly = TRUE)) {
    msg <- tr_("Package %s needed for this function to work. Please install it.")
    stop(sprintf(msg, sQuote(x)), call. = FALSE)
  }
  invisible(NULL)
}

# PAIRS
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_pairs
#' @aliases ternary_pairs,matrix-method
setMethod(
  f = "ternary_pairs",
  signature = c(x = "matrix"),
  definition = function(x, margin = NULL, ...) {

    ## Save and restore graphical parameters
    old_par <- graphics::par(no.readonly = TRUE)
    on.exit(graphics::par(old_par), add = TRUE)

    ## Layout
    n <- ncol(x)
    ni <- seq_len(n)
    parts <- colnames(x) %||% paste0("X", ni)
    zlab <- "*"

    if (!is.null(margin)) {
      margin <- if (is.character(margin)) which(parts == margin) else as.integer(margin)
      star <- x[, margin]
      zlab <- parts[margin]
      parts <- parts[-margin]
      n <- n - 1
      ni <- ni[-margin]
    }

    k <- (n^2 - n) / 2
    lay <- matrix(0, nrow = n, ncol = n)
    lay[lower.tri(lay, diag = FALSE)] <- seq_len(k)
    diag(lay) <- seq(k + 1, k + n)
    lay <- t(lay)

    graphics::layout(lay)
    graphics::par(mar = c(0, 0, 0, 0) + 0.1, oma = c(1, 1, 1, 1))

    ## Ternary plots
    p <- utils::combn(ni, 2, simplify = FALSE)
    for (i in p) {
      if (is.null(margin)) {
        star <- apply(X = x[, -i, drop = FALSE], MARGIN = 1, FUN = gmean)
      }
      z <- cbind(x[, i], star)
      ternary_plot(z, zlab = zlab, axes = FALSE, frame.plot = TRUE, ...)
    }

    ## Graphical parameters
    str.wid <- max(graphics::strwidth(parts, "user"))
    cex.lab <- list(...)$cex.lab %||% max(0.8, min(2, 0.9 / str.wid))
    col.lab <- list(...)$col.lab %||% graphics::par("col.lab")
    font.lab <- list(...)$font.lab %||% graphics::par("font.lab")

    ## Add labels
    for (part in parts) {
      graphics::plot(NULL, xlim = c(0, 1), ylim = c(0, 1),
                     axes = FALSE, ann = FALSE)
      graphics::text(x = 0.5, y = 0.5, labels = part,
                     cex = cex.lab, col = col.lab, font = font.lab)
    }

    invisible(x)
  }
)

#' @export
#' @rdname ternary_pairs
#' @aliases ternary_pairs,data.frame,missing,missing-method
setMethod(
  f = "ternary_pairs",
  signature = c(x = "data.frame"),
  definition = function(x, margin = NULL, ...) {
    x <- data.matrix(x)
    methods::callGeneric(x = x, margin = margin, ...)

    invisible(x)
  }
)

# TERNARY GRID
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_grid
ternary_grid <- function(primary = NULL, secondary = NULL,
                         center = NULL, scale = NULL,
                         col.primary = "darkgray", col.secondary = "lightgray",
                         lty.primary = "dashed", lty.secondary = "dotted",
                         lwd.primary = 1, lwd.secondary = lwd.primary) {
  ## Center and scale
  if (is.null(center)) center <- get_center()
  if (is.null(scale)) scale <- get_scale()

  ## Primary grid
  if (is.null(primary) || (!anyNA(primary) && length(primary) == 1 && primary >= 1)) {
    if (is.null(primary)) primary <- graphics::par("xaxp")[3L]
    i <- seq(from = 0, to = 1, length.out = primary + 1)
    .ternary_grid(i, center = center, scale = scale,
                  col = col.primary, lty = lty.primary, lwd = lwd.primary)
  }

  ## Secondary grid
  if (!is.null(secondary) && !is.na(secondary) && length(secondary) == 1 && secondary > primary) {
    i <- seq(from = 0, to = 1, length.out = secondary + 1)
    .ternary_grid(i, center = center, scale = scale,
                  col = col.secondary, lty = lty.secondary, lwd = lwd.secondary)
  }

  invisible(NULL)
}

.ternary_grid <- function(x, center = NULL, scale = NULL,
                          col = "lightgray", lty = "dotted", lwd = 1, n = 100) {
  ## Reset values if needed
  if (!is.null(center) && all(center == 1)) center <- NULL
  if (!is.null(scale) && scale == 1) scale <- NULL

  x <- x[!(x == 0 | x == 1)]
  if (is.null(scale)) {
    for (i in x) {
      start <- matrix(data = c(i, 0, 1 - i, 1 - i, i, 0, 0, 1 - i, i), ncol = 3)
      end <- matrix(data = c(i, 1 - i, 0, 0, i, 1 - i, 1 - i, 0, i), ncol = 3)

      start <- coordinates_ternary(start, center = center)
      end <- coordinates_ternary(end, center = center)

      graphics::segments(
        x0 = start$x, x1 = end$x,
        y0 = start$y, y1 = end$y,
        lty = lty, lwd = lwd, col = col
      )
    }
  } else {
    for (i in x) {
      start <- matrix(data = c(i, 0, 1 - i, 1 - i, i, 0, 0, 1 - i, i), ncol = 3)
      end <- matrix(data = c(i, 1 - i, 0, 0, i, 1 - i, 1 - i, 0, i), ncol = 3)
      start <- list(x = start[, 2] + start[, 3] / 2, y = start[, 3] * sqrt(3) / 2)
      end <- list(x = end[, 2] + end[, 3] / 2, y = end[, 3] * sqrt(3) / 2)

      mapply(
        FUN = function(x_from, x_to, y_from, y_to, n, center, scale) {
          x <- seq(x_from, x_to, length.out = n)
          y <- seq(y_from, y_to, length.out = n)
          z <- coordinates_cartesian(x, y)
          zz <- coordinates_ternary(z, center = center, scale = scale)
          graphics::lines(
            zz,
            lty = lty, lwd = lwd, col = col
          )
        },
        x_from = start$x, x_to = end$x,
        y_from = start$y, y_to = end$y,
        MoreArgs = list(n = 100, center = center, scale = scale)
      )
    }
  }
}

# TERNARY AXES
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_axis
ternary_axis <- function(side, at = NULL, labels = TRUE, tick = TRUE,
                         center = NULL, scale = NULL,
                         font = NA, lty = "solid",
                         lwd = 1, lwd.ticks = lwd,
                         col = NULL, col.ticks = NULL, ...) {
  ## Center and scale
  if (is.null(center)) center <- get_center()
  if (is.null(scale)) scale <- get_scale()

  ## Graphical parameters
  if (is.na(font)) font <- list(...)$font.axis %||% graphics::par("font.axis")
  if (is.null(col)) col <- list(...)$col.axis %||% graphics::par("col.axis")
  if (is.null(col.ticks)) col.ticks <- col
  cex <- list(...)$cex.axis %||% graphics::par("cex.axis")
  tcl <- list(...)$tcl %||% graphics::par("tcl")

  ## Ticks and labels position
  if (is.null(at)) {
    at <- seq(from = 0, to = 1, length.out = graphics::par("xaxp")[3L] + 1)
    at <- at[!(at == 0 | at == 1)]
  }

  axis_degree <- c(120, 240, 0)[side]
  axis_radian <- c(0, 2 * pi / 3, 4 * pi / 3)[side]

  pos <- matrix(data = 0, nrow = length(at), ncol = 3)
  pos[, side] <- at
  pos[, c(2, 3, 1)[side]] <- 1 - at
  pos <- coordinates_ternary(pos, center = center, scale = scale)

  h <- abs(tcl * graphics::strheight("1", cex = 1))
  dx <- sin(pi / 6) * h
  dy <- cos(pi / 6) * h
  tick_start <- matrix(c(pos$x, pos$y), ncol = 2)
  tick_end <- matrix(c(pos$x + dx * c(1, 1, -2)[side],
                       pos$y + dy * c(-1, 1, 0)[side]), ncol = 2)

  ## Labels
  if (!isFALSE(labels)) {
    if (length(labels) != length(at)) labels <- round(at * 100)
    graphics::text(x = tick_end, label = labels, srt = axis_degree,
                   cex = cex, col = col, font = font, adj = c(1, 0.5))
  }

  ## Ticks
  if (tick) {
    axis_line <- rotate(matrix(c(0, 0, 1, 0), ncol = 2), theta = axis_radian)
    graphics::segments(
      x0 = axis_line[1, 1], x1 = axis_line[2, 1],
      y0 = axis_line[1, 2], y1 = axis_line[2, 2],
      col = col, lwd = lwd, lty = lty
    )
    graphics::segments(
      x0 = tick_start[, 1], x1 = tick_end[, 1],
      y0 = tick_start[, 2], y1 = tick_end[, 2],
      col = col.ticks, lwd = lwd.ticks, lty = lty
    )
  }

  invisible(NULL)
}

# TERNARY SOIL DIAGRAM
#' @include AllGenerics.R
NULL

#' @export
#' @rdname triangle_soil
triangle_soil_hypres <- function(labels = TRUE, symbol = FALSE, ...) {
  .triangle_soil("HYPRES", labels = labels, symbol = symbol, ...)
  invisible(NULL)
}

#' @export
#' @rdname triangle_soil
triangle_soil_folk <- function(labels = TRUE, symbol = FALSE, ...) {
  .triangle_soil("folk1954", labels = labels, symbol = symbol, ...)
  invisible(NULL)
}

#' @export
#' @rdname triangle_soil
triangle_soil_shepard <- function(labels = TRUE, symbol = FALSE, ...) {
  .triangle_soil("shepard1954", labels = labels, symbol = symbol, ...)
  invisible(NULL)
}

#' @export
#' @rdname triangle_soil
triangle_soil_usda <- function(labels = TRUE, symbol = FALSE, ...) {
  .triangle_soil("USDA1951", labels = labels, symbol = symbol, ...)
  invisible(NULL)
}

.triangle_soil <- function(chart, labels = TRUE, symbol = FALSE, ...) {
  ## 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")

  poly <- .soil[[chart]]
  txt_lab <- unique(poly$label)
  txt_symb <- unique(poly$symbol)

  poly$label <- factor(poly$label, levels = unique(poly$label))
  poly <- split(poly, f = poly$label)

  for (i in poly) {
    ternary_polygon(i, ...)
  }
  if (labels) {
    lab <- lapply(X = poly, FUN = function(x) colMeans(x[, c(1, 2, 3, 4)]))
    txt <- if (symbol && !all(txt_symb == "")) txt_symb else txt_lab

    for (i in seq_along(lab)) {
      ternary_text(
        x = lab[[i]][[1]],
        y = lab[[i]][[2]],
        z = lab[[i]][[3]],
        labels = txt[[i]],
        srt = lab[[i]][[4]],
        cex = cex.lab,
        col = col.lab,
        font = font.lab
      )
    }
  }
}

# ELLIPSE
#' @include AllGenerics.R
NULL

# Ellipse ======================================================================
#' @export
#' @rdname ternary_ellipse
#' @aliases ternary_ellipse,numeric,numeric,numeric-method
setMethod(
  f = "ternary_ellipse",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, radius = 1, ...) {
    ## ILR
    coda <- cbind(x, y, z)
    ratio <- ilr(coda)

    ## Compute ellipse
    mu <- colMeans(ratio)
    sigma <- stats::cov(ratio)
    # rob <- robustbase::covMcd(ratio)
    # mu <- rob$center
    # sigma <- rob$cov
    xy <- ellipse(sigma = sigma, mu = mu, radius = radius)

    for (i in seq_along(xy)) {
      tern <- ilr_inv(xy[[i]]) # Inverse transform
      coords <- coordinates_ternary(tern)
      graphics::polygon(x = coords$x, y = coords$y, ...)
    }

    invisible(list(x = x, y = y, z = z))
  }
)

#' @export
#' @rdname ternary_ellipse
#' @aliases ternary_ellipse,ANY,missing,missing-method
setMethod(
  f = "ternary_ellipse",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, radius = 1, ...) {
    xyz <- grDevices::xyz.coords(x)
    pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
                               radius = radius, ...)
    invisible(pt)
  }
)

#' Computes an Ellipse
#'
#' @param sigma A square positive definite \eqn{2 \times 2}{2 x 2} covariance
#'  or correlation `matrix`.
#' @param mu A length-two [`numeric`] vector giving the centre of the ellipse.
#' @param scale If `sigma` is a correlation matrix, then the standard deviations
#'  of each parameter can be given in the scale parameter.
#'  Defaults to `c(1, 1)`, so no rescaling will be done.
#' @param level A length-\eqn{k} [`numeric`] vector giving the confidence level
#'  of a pairwise confidence region.
#' @param radius The size of the ellipse may also be controlled by specifying
#'  the value of a t-statistic on its boundary.
#' @param n A length-one [`numeric`] vector specifying the number of points used
#'  in the ellipse.
#' @note Adapted from [ellipse::ellipse()].
#' @return
#'  A [`list`] of \eqn{k} \eqn{n \times 2}{n x 2} `matrix`, suitable for
#'  plotting.
#' @keywords internal
#' @noRd
ellipse <- function(sigma, mu = c(0, 0), scale = c(1, 1), level = 0.95,
                    radius = sqrt(stats::qchisq(level, 2)), n = 100, ...) {
  r <- sigma[1, 2]

  if (missing(scale)) {
    scale <- sqrt(diag(sigma))
    if (scale[1] > 0) r <- r / scale[1]
    if (scale[2] > 0) r <- r / scale[2]
  }

  r <- min(max(r, -1), 1)  # clamp to -1..1, in case of rounding errors
  d <- acos(r)
  a <- seq(0, 2 * pi, len = n)

  lapply(
    X = radius,
    FUN = function(x) {
      matrix(
        data = c(x * scale[1] * cos(a + d / 2) + mu[1],
                 x * scale[2] * cos(a - d / 2) + mu[2]),
        nrow = n,
        ncol = 2,
        dimnames = list(NULL, c("x", "y"))
      )
    }
  )
}

# Confidence ellipse ===========================================================
#' @export
#' @rdname ternary_ellipse
#' @aliases ternary_confidence,numeric,numeric,numeric-method
setMethod(
  f = "ternary_confidence",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, level = 0.95, ...) {
    df1 <- 2
    df2 <- length(x) - 2
    radius <- sqrt(stats::qf(p = level, df1, df2) * df1 / df2)
    ternary_ellipse(x, y, z, radius = radius, ...)
  }
)

#' @export
#' @rdname ternary_ellipse
#' @aliases ternary_confidence,ANY,missing,missing-method
setMethod(
  f = "ternary_confidence",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, level = 0.95, ...) {
    x <- grDevices::xyz.coords(x)
    methods::callGeneric(x = x$x, y = x$y, z = x$z, level = level, ...)
  }
)

# Probability ellipse ==========================================================
#' @export
#' @rdname ternary_ellipse
#' @aliases ternary_tolerance,numeric,numeric,numeric-method
setMethod(
  f = "ternary_tolerance",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, level = 0.95, ...) {
    radius <- sqrt(stats::qchisq(p = level, df = 2))
    ternary_ellipse(x, y, z, radius = radius, ...)
  }
)

#' @export
#' @rdname ternary_ellipse
#' @aliases ternary_tolerance,ANY,missing,missing-method
setMethod(
  f = "ternary_tolerance",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, level = 0.95, ...) {
    x <- grDevices::xyz.coords(x)
    methods::callGeneric(x = x$x, y = x$y, z = x$z, level = level, ...)
  }
)

# IMAGE
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_image
#' @aliases ternary_image,function-method
setMethod(
  f = "ternary_image",
  signature = c(f = "function"),
  definition = function(f, n = 48, palette = NULL, ...) {

    tri <- .triangle_center(n)
    xyz <- coordinates_cartesian(tri$x, tri$y)
    val <- f(xyz$x, xyz$y, xyz$z, ...)
    col <- map_color(val, palette = palette)

    coords <- .triangle_vertex(tri$x, tri$y, tri$direction, tri$resolution)
    for (i in seq_along(coords)) {
      polygon(coords[[i]][, 1], coords[[i]][, 2], col = col[i], border = NA)
    }

    invisible(list(levels = val, colors = col))
  }
)

#' @export
#' @rdname ternary_tile
#' @aliases tile_bin,numeric,numeric,numeric-method
setMethod(
  f = "tile_bin",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z) {

    total <- x + y + z
    a <- x / total
    b <- y / total
    c <- z / total

    function(x, y, z) {
      tri <- .triangle_center(sqrt(length(x)))
      coords <- .triangle_vertex(tri$x, tri$y, tri$direction, tri$resolution)

      count <- numeric(length(coords))
      for (i in seq_along(coords)) {
        xyz <- coordinates_cartesian(coords[[i]][, 1], coords[[i]][, 2])
        count[[i]] <- sum(min(xyz$x) <= a & a < max(xyz$x) &
                            min(xyz$y) <= b & b < max(xyz$y) &
                            min(xyz$z) <= c & c < max(xyz$z))
      }

      count[count == 0] <- NA
      count
    }
  }
)

#' @export
#' @rdname ternary_tile
#' @aliases tile_bin,ANY,missing,missing-method
setMethod(
  f = "tile_bin",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x) {
    xyz <- grDevices::xyz.coords(x)
    methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z)
  }
)

#' @export
#' @rdname ternary_tile
#' @aliases tile_density,numeric,numeric,numeric-method
setMethod(
  f = "tile_density",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z) {
    ## ILR
    coda <- cbind(x, y, z)
    ratio <- ilr(coda)

    ## Compute KDE
    function(x, y, z) {
      xyz <- cbind(x, y, z)
      xy <- ilr(xyz)
      dens <- kde(
        x = ratio[, 1],
        y = ratio[, 2],
        gx = sort(unique(xy[, 1])),
        gy = sort(unique(xy[, 2]))
      )

      i <- as.numeric(as.factor(rank(xy[, 1])))
      j <- as.numeric(as.factor(rank(xy[, 2])))

      dens$z[cbind(i, j)]
    }
  }
)

#' @export
#' @rdname ternary_tile
#' @aliases tile_density,ANY,missing,missing-method
setMethod(
  f = "tile_density",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x) {
    xyz <- grDevices::xyz.coords(x)
    methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z)
  }
)

#' @export
#' @rdname ternary_tile
#' @aliases tile_interpolate,numeric,numeric,numeric-method
setMethod(
  f = "tile_interpolate",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, value, method = "linear", ...) {
    ## Validation
    assert_package("interp")
    assert_length(value, length(x))

    ## ILR
    coda <- cbind(x, y, z)
    ratio <- ilr(coda)

    ## Interpolate
    function(x, y, z) {
      xyz <- cbind(x, y, z)
      xy <- ilr(xyz)

      interp <- interp::interp(
        x = ratio[, 1],
        y = ratio[, 2],
        z = value,
        xo = sort(unique(xy[, 1])),
        yo = sort(unique(xy[, 2])),
        method = method,
        ...
      )

      i <- as.numeric(as.factor(rank(xy[, 1])))
      j <- as.numeric(as.factor(rank(xy[, 2])))

      interp$z[cbind(i, j)]
    }
  }
)

#' @export
#' @rdname ternary_tile
#' @aliases tile_interpolate,ANY,missing,missing-method
setMethod(
  f = "tile_interpolate",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, value, method = "linear", ...) {
    xyz <- grDevices::xyz.coords(x)
    methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
                         value = value, method = method, ...)
  }
)

#' Tile Center Coordinates
#'
#' Computes tile center cartesian coordinates.
#' @param resolution A length-one [`integer`] vector specifying the maximum
#'  number of tiles on each axis.
#' @return
#'  A [`list`] with the following elements:
#'  \describe{
#'   \item{`x`}{x cartesian coordinates.}
#'   \item{`y`}{y cartesian coordinates.}
#'   \item{`direction`}{`1` means up, `-1` means down.}
#'   \item{`resolution`}{}
#'  }
#' @examples
#' .triangle_center(5)
#' @keywords internal
#' @noRd
.triangle_center <- function(resolution) {

  offset <- 1 / resolution / 2L
  height <- the$top / resolution

  X <- seq(from = offset, to = 1 - offset, by = offset)
  X <- lapply(
    X = seq_len(resolution),
    FUN = function(step) {
      up <- X[seq(from = step, to = (2L * resolution) - step, by = 2L)]
      down <- integer(0)
      if (step != resolution) {
        down <- X[seq(from = step + 1, to = (2L * resolution) - step - 1, by = 2L)]
      }
      c(rbind(up, c(down, NA)))[-(resolution - step + 1) * 2L]
    }
  )

  in_row <- 2L * rev(seq_len(resolution)) - 1L
  is_down <- (1 + unlist(lapply(in_row, seq_len))) %% 2L

  Y <- seq(from = height / 3, to = the$top - (2 * height / 3), length.out = resolution)
  Y <- rep(Y[seq_len(resolution)], in_row)
  Y <- Y + (is_down * height / 3)

  dir <- rep(1, length(is_down))
  dir[is_down == 1] <- -1

  list(
    x = unlist(X),
    y = Y,
    direction = dir,
    resolution = resolution
  )
}

#' Tile Vertex Coordinates
#'
#' Computes tile vertex cartesian coordinates.
#' @param x,y A [`numeric`] vector giving the cartesian coordinates of the
#'  center.
#' @param direction An [`integer`] vector specifying the triangle direction
#'  (`1` means up, `-1` means down).
#' @param resolution A length-one [`integer`] vector specifying the maximum
#'  number of tiles on each axis.
#' @return
#'  A [`list`] of `numeric` [`matrix`].
#' @examples
#' m <- .triangle_center(5)
#' .triangle_vertex(m$x, m$y, m$direction, m$resolution)
#' @keywords internal
#' @noRd
.triangle_vertex <- function(x, y, direction, resolution) {
  n <- length(x)

  width <- 1 / resolution / 2
  height <- the$top / resolution / 3

  tiles <- vector(mode = "list", length = n)
  for (i in seq_len(n)) {
    tiles[[i]] <- matrix(
      data = c(x[i] + c(0, width, -width), y[i] + c(2 * height, -height, -height) * direction[i]),
      ncol = 2, dimnames = list(NULL, c("x", "y"))
    )
  }

  tiles
}

# TERNARY WINDOW
#' @include AllGenerics.R
NULL

#' Set up Ternary Coordinates for Graphics Window
#'
#' @param xlim A length-three [`numeric`] vector giving the `x` limits in the
#'  range \eqn{[0,1]}.
#' @param ylim A length-three [`numeric`] vector giving the `y` limits in the
#'  range \eqn{[0,1]}.
#' @param zlim A length-three [`numeric`] vector giving the `z` limits in the
#'  range \eqn{[0,1]}.
#' @param xlab,ylab,zlab A [`character`] string giving a label for the x, y and
#'  z axes.
#' @return
#'  `ternary_window()` is called it for its side-effects
#'  (see [graphics::plot.window()]).
#' @keywords internal
#' @noRd
ternary_window <- function(xlim = NULL, ylim = NULL, zlim = NULL,
                           xlab = NULL, ylab = NULL, zlab = NULL,
                           cex = 1) {

  n_null <- is.null(xlim) + is.null(ylim) + is.null(zlim)

  if (n_null == 3) {
    dx <- max(graphics::strwidth(c(xlab, ylab, zlab), cex = cex)) / 2
    rx <- expand_range(c(0, 1), add = dx)
    lim <- list(
      x = rx,
      y = the$top / 2 + c(-1, 1) * diff(rx) / 2
    )
  }
  if (n_null == 2) {
    stop(tr_("You must provide at least two coordinates ranges."), call. = FALSE)
  }
  if (n_null < 2) {
    xlims <- if (is.null(xlim)) 1 - ylim - zlim else xlim
    ylims <- if (is.null(ylim)) 1 - xlim - zlim else ylim
    zlims <- if (is.null(zlim)) 1 - xlim - ylim else zlim

    assert_length(xlims, 3)
    assert_length(ylims, 3)
    assert_length(zlims, 3)

    lim <- coordinates_ternary(abs(xlims), abs(ylims), abs(zlims))
    lim$x <- range(lim$x)
    lim$y <- range(lim$y)
  }
  graphics::plot.window(xlim = lim$x, ylim = lim$y, asp = 1)
}

# GENERIC METHODS
#' @include aaa.R
NULL

# Coordinates ==================================================================
#' Ternary Coordinates
#'
#' Computes ternary coordinates.
#' @param x,y,z A [`numeric`] vector giving the x, y and z Cartesian coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param center A [`logical`] scalar or a [`numeric`] vector giving the center.
#' @param scale A [`logical`] scalar or a length-one [`numeric`] vector giving a
#'  scaling factor.
#' @param xlab,ylab,zlab A [`character`] string specifying the names for the x,
#'  y and z variables to be extracted.
#' @param missing A [`logical`] scalar: should [missing values][NA] be replaced
#'  with zeros before the computation proceeds? If `FALSE` (the default),
#'  incomplete cases are removed.
#' @param ... Currently not used.
#' @return
#'  A [`list`] with the components:
#'  \describe{
#'   \item{`x`}{A [`numeric`] vector of x coordinates.}
#'   \item{`y`}{A [`numeric`] vector of y coordinates.}
#'   \item{`center`}{A [`numeric`] vector giving the center.}
#'   \item{`scale`}{A [`numeric`] vector giving the scale factor.}
#'  }
#' @example inst/examples/ex-coordinates.R
#' @author N. Frerebeau
#' @docType methods
#' @family coordinates
#' @aliases coordinates_ternary-method
#' @keywords internal
setGeneric(
  name = "coordinates_ternary",
  def = function(x, y, z, ...) standardGeneric("coordinates_ternary"),
  valueClass = "list"
)

#' Cartesian Coordinates
#'
#' Computes Cartesian coordinates.
#' @param x,y A [`numeric`] vector giving the x and y ternary coordinates of a
#'  set of points. If `y` is missing, an attempt is made to interpret `x` in a
#'  suitable way (see [grDevices::xy.coords()]).
#' @param xlab,ylab A [`character`] string specifying the names for the x and y
#'  variables to be extracted.
#' @param ... Currently not used.
#' @return
#'  A [`list`] with the components:
#'  \describe{
#'   \item{`x`}{A [`numeric`] vector of x coordinates.}
#'   \item{`y`}{A [`numeric`] vector of y coordinates.}
#'   \item{`z`}{A [`numeric`] vector of z coordinates.}
#'  }
#' @example inst/examples/ex-coordinates.R
#' @author N. Frerebeau
#' @docType methods
#' @family coordinates
#' @aliases coordinates_cartesian-method
#' @keywords internal
setGeneric(
  name = "coordinates_cartesian",
  def = function(x, y, ...) standardGeneric("coordinates_cartesian"),
  valueClass = "list"
)

# Plot =========================================================================
#' Ternary Plot
#'
#' Produces a ternary plot.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param center A [`logical`] scalar: should the data be centered?
#' @param scale A [`logical`] scalar: should the data be scaled?
#' @param xlim A length-three [`numeric`] vector giving the `x` limits in the
#'  range \eqn{[0,1]}.
#' @param ylim A length-three [`numeric`] vector giving the `y` limits in the
#'  range \eqn{[0,1]}.
#' @param zlim A length-three [`numeric`] vector giving the `z` limits in the
#'  range \eqn{[0,1]}.
#' @param xlab,ylab,zlab A [`character`] string giving a label for the x, y and
#'  z axes.
#' @param main A [`character`] string giving a main title for the plot.
#' @param sub A [`character`] string giving a subtitle for the plot.
#' @param ann A [`logical`] scalar: should the default annotation (title and x,
#'  y and z axis labels) appear on the plot?
#' @param axes A [`logical`] scalar: should axes be drawn on the plot?
#' @param frame.plot A [`logical`] scalar: should a box be drawn around the
#'  plot?
#' @param panel.first An an `expression` to be evaluated after the plot axes are
#'  set up but before any plotting takes place. This can be useful for drawing
#'  background grids.
#' @param panel.last An `expression` to be evaluated after plotting has taken
#'  place but before the axes, title and box are added.
#' @param ... Other [graphical parameters][graphics::par] may also be passed as
#'  arguments to this function.
#' @return
#'  `ternary_plot()` is called it for its side-effects: it results in a graphic
#'  being displayed. Invisibly returns a [`list`] with the components:
#'  \describe{
#'   \item{`x`}{A [`numeric`] vector of x values.}
#'   \item{`y`}{A [`numeric`] vector of y values.}
#'   \item{`z`}{A [`numeric`] vector of z values.}
#'   \item{`center`}{A [`numeric`] vector giving the center.}
#'   \item{`scale`}{A [`numeric`] vector giving the scale factor.}
#'  }
#' @example inst/examples/ex-plot.R
#' @author N. Frerebeau
#' @docType methods
#' @family graphical elements
#' @aliases ternary_plot-method
setGeneric(
  name = "ternary_plot",
  def = function(x, y, z, ...) standardGeneric("ternary_plot")
)

## Grid ------------------------------------------------------------------------
#' Add Grid to a Ternary Plot
#'
#' Adds a triangular grid to an existing plot.
#' @param primary An [`integer`] specifying the number of cells of the primary
#'  grid in `x`, `y` and `z` direction.
#' @param secondary An [`integer`] specifying the number of cells of the
#'  secondary grid in `x`, `y` and `z` direction.
#' @param center A [`numeric`] vector giving the center. If `NULL` (the
#'  default), `center` is automatically set based on the last call to
#'  [plot()].
#' @param scale A [`numeric`] vector giving the scale factor. If `NULL` (the
#'  default), `scale` is automatically set based on the last call to
#'  [plot()].
#' @param col.primary,col.secondary A [`character`] string specifying the color
#'  of the grid lines.
#' @param lty.primary,lty.secondary A [`character`] string or [`numeric`]
#'  value specifying the line type of the grid lines.
#' @param lwd.primary,lwd.secondary A non-negative [`numeric`] value specifying
#'  the line width of the grid lines.
#' @return
#'  `ternary_grid()` is called it for its side-effects.
#' @example inst/examples/ex-scale.R
#' @author N. Frerebeau
#' @docType methods
#' @family graphical elements
#' @name ternary_grid
NULL

## Axis ------------------------------------------------------------------------
#' Add an Axis to a Ternary Plot
#'
#' Adds an axis to the current plot.
#' @param side An [`integer`] specifying which side of the plot the axis is to
#'  be drawn on. The axis is placed as follows: 1=below, 2=right and 3=left.
#' @param at A [`numeric`] vector giving the points at which tick-marks are to
#'  be drawn.
#' @param labels A [`logical`] scalar specifying whether (numerical) annotations
#'  are to be made at the tickmarks, or a [`character`] vector of labels to be
#'  placed at the tickpoints. If this is not `logical`, `at` should also be
#'  supplied and of the same length.
#' @param tick A [`logical`] scalar: should tickmarks and an axis line be drawn?
#' @param center A [`numeric`] vector giving the center. If `NULL` (the
#'  default), `center` is automatically set based on the last call to
#'  [plot()].
#' @param scale A [`numeric`] vector giving the scale factor. If `NULL` (the
#'  default), `scale` is automatically set based on the last call to
#'  [plot()].
#' @param font font for text. Defaults to `par("font.axis")`.
#' @param lty A [`character`] string or [`numeric`] value specifying the line
#'  type for both the axis line and the tick marks.
#' @param lwd,lwd.ticks A non-negative [`numeric`] value specifying the line
#'  widths for the axis line and the tick marks.
#' @param col,col.ticks Colors for the axis line and the tick marks
#'  respectively. Defaults to `par("col.axis")`.
#' @param ... Other [graphical parameters][graphics::par] may also be passed as
#'  arguments to this function, particularly, `cex.axis`, `col.axis` and
#'  `font.axis` for axis annotation.
#' @return
#'  `ternary_axis()` is called it for its side-effects.
#' @example inst/examples/ex-axis.R
#' @author N. Frerebeau
#' @docType methods
#' @family graphical elements
#' @name ternary_axis
NULL

## Annotation ------------------------------------------------------------------
#' Ternary Plot Annotation
#'
#' @param main A [`character`] string specifying the main title (on top).
#' @param sub A [`character`] string specifying the sub-title (at bottom).
#' @param xlab,ylab,zlab A [`character`] string giving a label for the x, y and
#'  z axes.
#' @param line Specifying a value for `line` overrides the default placement of
#'  labels, and places them this many lines outwards from the plot edge.
#' @param outer A [`logical`] scalar: should the titles be placed in the outer
#'  margins of the plot?
#' @param ... Other [graphical parameters][graphics::par] may also be passed as
#'  arguments to this function, particularly, `font.main`, `cex.main`,
#'  `col.main` and `font.sub`, `cex.sub`, `col.sub` for title annotation;
#'  `font.lab`, `cex.lab` and `col.lab` for axis label.
#' @return
#'  `ternary_title()` is called it for its side-effects.
#' @example inst/examples/ex-title.R
#' @author N. Frerebeau
#' @docType methods
#' @family graphical elements
#' @name ternary_title
NULL

## Box -------------------------------------------------------------------------
#' Draw a Box around a Ternary Plot
#'
#' @param lty A [`character`] string or [`numeric`] value specifying the line
#'  type of the box.
#' @param ... Other [graphical parameters][graphics::par] may also be passed as
#'  arguments to this function, particularly, `col` or `lwd`.
#' @return
#'  `ternary_box()` is called it for its side-effects.
#' @example inst/examples/ex-axis.R
#' @author N. Frerebeau
#' @docType methods
#' @family graphical elements
#' @name ternary_box
NULL

## Pairs -----------------------------------------------------------------------
#' Ternary Plot Matrices
#'
#' Produces a matrix of ternary plots.
#' @param x A [`matrix`] or a [`data.frame`]. Columns are converted to `numeric`
#'  in the same way that [data.matrix()] does.
#' @param margin A [`character`] string or an [`integer`] giving the index of
#'  the column to be used as the third part of the ternary plots. If `NULL`
#'  (the default), marginal compositions will be used (i.e. the geometric mean
#'  of the non-selected parts).
#' @param ... Further [graphical parameters][graphics::par()].
#' @return
#'  `ternary_pairs()` is called it for its side-effects: it results in a graphic
#'  being displayed. Invisibly returns `x`.
#' @example inst/examples/ex-pairs.R
#' @author N. Frerebeau
#' @docType methods
#' @family graphical elements
#' @aliases ternary_pairs-method
setGeneric(
  name = "ternary_pairs",
  def = function(x, ...) standardGeneric("ternary_pairs")
)

# Geometry =====================================================================
#' Add Arrows to a Ternary Plot
#'
#' Draw arrows between pairs of points.
#' @param x0,y0,z0 A [`numeric`] vector giving the x, y and z ternary
#'  coordinates of points from which to draw.
#' @param x1,y1,z1 A [`numeric`] vector giving the x, y and z ternary
#'  coordinates of points to which to draw.
#' @param ... Further arguments to be passed to [graphics::arrows()].
#' @return
#'  `ternary_arrows()` is called it for its side-effects.
#' @seealso [graphics::arrows()]
#' @example inst/examples/ex-arrows.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_arrows-method
setGeneric(
  name = "ternary_arrows",
  def = function(x0, y0, z0, ...) standardGeneric("ternary_arrows")
)

#' Add Line Segments to a Ternary Plot
#'
#' Draw line segments between pairs of points.
#' @param x0,y0,z0 A [`numeric`] vector giving the x, y and z ternary
#'  coordinates of points from which to draw.
#' @param x1,y1,z1 A [`numeric`] vector giving the x, y and z ternary
#'  coordinates of points to which to draw.
#' @param ... Further graphical parameters (see [graphics::par()]) may also be
#'  supplied as arguments, particularly, line type, `lty`, line width, `lwd` and
#'  color, `col`. Also the line characteristics `lend`, `ljoin` and `lmitre`.
#' @return
#'  `ternary_segments()` is called it for its side-effects.
#' @seealso [graphics::segments()]
#' @example inst/examples/ex-segments.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_segments-method
setGeneric(
  name = "ternary_segments",
  def = function(x0, y0, z0, ...) standardGeneric("ternary_segments")
)

#' Add Cross-Hairs to a Ternary Plot
#'
#' Draw lines that intersect at a point.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param x_mark,y_mark,z_mark A [`logical`] scalar: should the `x`, `y` or `z`
#'  axis component be drawn?
#' @param ... Further graphical parameters (see [graphics::par()]) may also be
#'  supplied as arguments, particularly, line type, `lty`, line width, `lwd` and
#'  color, `col`. Also the line characteristics `lend`, `ljoin` and `lmitre`.
#' @return
#'  `ternary_crosshairs()` is called it for its side-effects.
#' @example inst/examples/ex-crosshairs.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_crosshairs-method
setGeneric(
  name = "ternary_crosshairs",
  def = function(x, y, z, ...) standardGeneric("ternary_crosshairs")
)

#' Add Connected Line Segments to a Ternary Plot
#'
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param type A [`character`] string indicating the type of plotting; actually
#'  any of the types as in [graphics::plot.default()].
#' @param ... Further graphical parameters (see [graphics::par()]) may also be
#'  supplied as arguments, particularly, line type, `lty`, line width, `lwd`,
#'  color, `col` and for `type = "b"`, `pch`. Also the line characteristics
#'  `lend`, `ljoin` and `lmitre`.
#' @return
#'  `ternary_lines()` is called it for its side-effects.
#' @seealso [graphics::lines()]
#' @example inst/examples/ex-lines.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_lines-method
setGeneric(
  name = "ternary_lines",
  def = function(x, y, z, ...) standardGeneric("ternary_lines")
)

#' Add Points to a Ternary Plot
#'
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param center A [`logical`] scalar specifying wether the data should be
#'  centered, or a [`numeric`] vector giving the center.
#' @param scale A [`logical`] scalar specifying wether the data should be
#'  scaled, or a [`numeric`] vector giving the scale factor.
#' @param type A [`character`] string indicating the type of plotting; actually
#'  any of the types as in [graphics::plot.default()].
#' @param ... Further graphical parameters (see [graphics::par()]) may also be
#'  supplied as arguments, particularly, plotting character, `pch`, character
#'  expansion, `cex` and color, `col`.
#' @return
#'  `ternary_points()` is called it for its side-effects. Invisibly returns
#'  a [`list`] with the components:
#'  \tabular{ll}{
#'   `x` \tab A [`numeric`] vector of x values. \cr
#'   `y` \tab A [`numeric`] vector of y values. \cr
#'   `z` \tab A [`numeric`] vector of z values. \cr
#'   `center` \tab A [`numeric`] vector giving the center. \cr
#'   `scale` \tab A [`numeric`] vector giving the scale factor. \cr
#'  }
#' @seealso [graphics::points()]
#' @example inst/examples/ex-points.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_points-method
setGeneric(
  name = "ternary_points",
  def = function(x, y, z, ...) standardGeneric("ternary_points")
)

#' Polygon Drawing
#'
#' Draws the polygons whose vertices are given in `x`, `y` and `z`.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param ... Further arguments to be passed to [graphics::polygon()].
#' @return
#'  `ternary_polygon()` is called it for its side-effects.
#' @seealso [graphics::polygon()]
#' @example inst/examples/ex-polygon.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_polygon-method
setGeneric(
  name = "ternary_polygon",
  def = function(x, y, z, ...) standardGeneric("ternary_polygon")
)

#' Add Text to a Ternary Plot
#'
#' Draws the strings given in the vector `labels` at the coordinates given by
#' `x`, `y` and `z`.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param center A [`logical`] scalar specifying wether the data should be
#'  centered, or a [`numeric`] vector giving the center.
#' @param scale A [`logical`] scalar specifying wether the data should be
#'  scaled, or a [`numeric`] vector giving the scale factor.
#' @param labels A [`character`] vector or [`expression`] specifying the text
#'  to be written.
#' @param ... Further arguments to be passed to [graphics::text()].
#' @return
#'  `ternary_text()` is called it for its side-effects.
#' @seealso [graphics::text()]
#' @example inst/examples/ex-text.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_text-method
setGeneric(
  name = "ternary_text",
  def = function(x, y, z, ...) standardGeneric("ternary_text")
)

#' Non-Overlapping Text Labels
#'
#' Optimize the location of text labels to minimize overplotting text.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param center A [`logical`] scalar specifying wether the data should be
#'  centered, or a [`numeric`] vector giving the center.
#' @param scale A [`logical`] scalar specifying wether the data should be
#'  scaled, or a [`numeric`] vector giving the scale factor.
#' @param labels A [`character`] vector or [`expression`] specifying the text
#'  to be written.
#' @param type A [`character`] string specifying the shape of the field.
#'  It must be one of "`text`" or "`shadow`". Any unambiguous substring
#'  can be given.
#' @param ... Further graphical parameters (see [graphics::par()]) may also be
#'  supplied as arguments, particularly, character expansion, `cex` and
#'  color, `col`.
#' @return
#'  `ternary_labels()` is called it for its side-effects.
#' @seealso [graphics::text()]
#' @example inst/examples/ex-labels.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_labels-method
setGeneric(
  name = "ternary_labels",
  def = function(x, y, z, ...) standardGeneric("ternary_labels")
)

## Image -----------------------------------------------------------------------
#' Display a Color Image
#'
#' Creates a grid of colored triangles with colors corresponding to the output
#' of a function.
#' @param f A [`function`] that takes three arguments (x, y and z coordinates)
#'  and returns a `numeric` vector (see [tile_bin()], [tile_density()],
#'  [tile_interpolate()]).
#' @param n A length-one [`integer`] vector specifying the maximum number of
#'  tiles on each axis.
#' @param palette A [`function`] that takes a single `numeric` vector
#'  (the output of `f`) as argument and returns a vector of color.
#'  If `NULL`, the default color scheme will be used. If `FALSE`, the output
#'  of `f` is used as colors.
#' @param ... Further parameters to be passed to `f`.
#' @return
#'  `ternary_image()` is called it for its side-effects.
#'
#'  Invisibly returns a [`list`] with elements `levels` and `colors` that can be
#'  used for a legend.
#' @example inst/examples/ex-image.R
#' @author N. Frerebeau
#' @docType methods
#' @family geometries
#' @aliases ternary_image-method
setGeneric(
  name = "ternary_image",
  def = function(f, ...) standardGeneric("ternary_image")
)

## Tiles -----------------------------------------------------------------------
#' Ternary Tiles
#'
#' Compute tile values.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param value A [`numeric`] vector giving the values to be interpolated.
#' @param method A [`character`] string: specifying the method for interpolation
#'  (see [interp::interp()]).
#' @param ... Further parameters to be passed to internal methods.
#' @return
#'  A [`function`] that takes three [`numeric`] vector as arguments and returns
#'  a `numeric` vector.
#' @example inst/examples/ex-tile.R
#' @seealso [ternary_image()]
#' @author N. Frerebeau
#' @docType methods
#' @family tiles
#' @name ternary_tile
#' @rdname ternary_tile
NULL

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

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

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

# Statistics ===================================================================
## Ellipse ---------------------------------------------------------------------
#' Add an Ellipse to a Ternary Plot
#'
#' Computes and draws a confidence/tolerance ellipse.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param radius A [`numeric`] vector specifying the scaling of the
#'  half-diameters.
#' @param level A [`numeric`] vector specifying the confidence/tolerance level.
#' @param ... Further arguments to be passed to [graphics::polygon()].
#' @details
#'  Ellipse coordinates are computed after an isometric log ratio transformation
#'  of the original data.
#' @return
#'  `ternary_ellipse()` is called it for its side-effects.
#' @seealso [graphics::polygon()]
#' @example inst/examples/ex-ellipse.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases ternary_ellipse-method
setGeneric(
  name = "ternary_ellipse",
  def = function(x, y, z, ...) standardGeneric("ternary_ellipse")
)

#' @rdname ternary_ellipse
setGeneric(
  name = "ternary_confidence",
  def = function(x, y, z, ...) standardGeneric("ternary_confidence")
)

#' @rdname ternary_ellipse
setGeneric(
  name = "ternary_tolerance",
  def = function(x, y, z, ...) standardGeneric("ternary_tolerance")
)

## Convex hull -----------------------------------------------------------------
#' Convex Hull of a Set of Points
#'
#' Computes and draws the convex hull of the set of points specified.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param center A [`logical`] scalar specifying wether the data should be
#'  centered, or a [`numeric`] vector giving the center.
#' @param scale A [`logical`] scalar specifying wether the data should be
#'  scaled, or a [`numeric`] vector giving the scale factor.
#' @param ... Further arguments to be passed to [graphics::polygon()].
#' @return
#'  `ternary_hull()` is called it for its side-effects.
#' @seealso [grDevices::chull()], [graphics::polygon()]
#' @example inst/examples/ex-hull.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases ternary_hull-method
setGeneric(
  name = "ternary_hull",
  def = function(x, y, z, ...) standardGeneric("ternary_hull")
)

## Contour ---------------------------------------------------------------------
#' Contour Lines
#'
#' Computes and draws contour lines.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param value A [`numeric`] vector giving the values to be interpolated.
#' @param n A length-one [`numeric`] specifying the number of grid points.
#' @param nlevels A length-one [`numeric`] vector specifying the number of
#'  contour levels desired. Only used if `levels` is `NULL`.
#' @param levels A [`numeric`] vector of levels at which to draw contour lines.
#' @param palette A color palette [`function`] that takes a single integer
#'  argument (the number of levels) and returns a vector of colors.
#' @param ilr A [`logical`] scalar: should interpolation be computed in ILR
#'  space? If `FALSE`, interpolation is computed in Cartesian space.
#' @param method A [`character`] string: specifying the method for interpolation
#'  (see [interp::interp()]).
#' @param extrapolate A [`logical`] scalar: should extrapolation be used outside
#'  of the convex hull determined by the data points (see [interp::interp()])?
#' @param ... Further arguments to be passed to [ternary_lines()].
#' @details
#'  Contour are computed from a bivariate interpolation onto a grid,
#'  after an isometric log ratio transformation of the original data.
#' @return
#'  `ternary_contour()` is called it for its side-effects.
#'
#'  Invisibly returns a [`list`] with elements `levels` (the contour levels) and
#'  `colors` (the contour colors) that can be used for a legend.
#' @note
#'  The \pkg{interp} package needs to be installed on your machine.
#' @seealso [interp::interp()], [grDevices::contourLines()]
#' @example inst/examples/ex-contour.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases ternary_contour-method
setGeneric(
  name = "ternary_contour",
  def = function(x, y, z, ...) standardGeneric("ternary_contour")
)

## Density ---------------------------------------------------------------------
#' Density Contour Lines
#'
#' Computes and draws density contour lines.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param h A length-one [`numeric`] vector giving the bandwidth.
#' @param n A length-one [`numeric`] specifying the number of grid points.
#' @param nlevels A length-one [`numeric`] vector specifying the number of
#'  contour levels desired. Only used if `levels` is `NULL`.
#' @param levels A [`numeric`] vector of levels at which to draw contour lines.
#' @param palette A color palette [`function`] that takes a single integer
#'  argument (the number of levels) and returns a vector of colors.
#' @param ... Further arguments to be passed to [ternary_lines()].
#' @details
#'  Two-dimensional kernel density estimation with an axis-aligned bivariate
#'  normal kernel. Normal kernel is evaluated on a square grid, after an
#'  isometric log ratio transformation of the original data.
#' @return
#'  `ternary_density()` is called it for its side-effects.
#'
#'  Invisibly returns a [`list`] with elements `levels` (the contour levels) and
#'  `colors` (the contour colors) that can be used for a legend.
#' @note
#'  **This must be considered as experimental and subject to major changes
#'  in a future release.**
#' @source
#'  Two-dimensional kernel density estimation is adapted from [`MASS::kde2d()`].
#' @seealso [grDevices::contourLines()]
#' @example inst/examples/ex-density.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases ternary_density-method
setGeneric(
  name = "ternary_density",
  def = function(x, y, z, ...) standardGeneric("ternary_density")
)

## Mean ------------------------------------------------------------------------
#' Compositional Mean
#'
#' Computes and draws the closed geometric mean of the set of points specified.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param ... Further arguments to be passed to [graphics::points()].
#' @return
#'  `ternary_mean()` is called it for its side-effects.
#' @example inst/examples/ex-mean.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases ternary_mean-method
setGeneric(
  name = "ternary_mean",
  def = function(x, y, z, ...) standardGeneric("ternary_mean")
)

## PCA -------------------------------------------------------------------------
#' Principal Component Analysis
#'
#' Computes and draws principal component.
#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
#'  of a set of points. If `y` and `z` are missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
#' @param axis An [`integer`] specifying the dimension to be plotted.
#' @param ... Further arguments to be passed to [graphics::lines()].
#' @return
#'  `ternary_pca()` is called it for its side-effects.
#' @example inst/examples/ex-pca.R
#' @author N. Frerebeau
#' @docType methods
#' @family statistics
#' @aliases ternary_pca-method
setGeneric(
  name = "ternary_pca",
  def = function(x, y, z, ...) standardGeneric("ternary_pca")
)

# Chart ========================================================================
#' Ceramic Phase Diagram
#'
#' @param labels A [`logical`] scalar: should labels be displayed?
#' @param symbol A [`logical`] scalar: should symbol be used instead of full
#'  labels? Only used if `labels` is `TRUE`.
#' @param mol A [`logical`] scalar: should molarity be used instead of molar
#'  mass?
#' @param ... Further arguments to be passed to [graphics::polygon()].
#' @example inst/examples/ex-phases.R
#' @author N. Frerebeau
#' @docType methods
#' @family charts
#' @name triangle_phase_cas
NULL

#' Soil Texture Triangle
#'
#' @param labels A [`logical`] scalar: should labels be displayed?
#' @param symbol A [`logical`] scalar: should symbol be used instead of full
#'  labels? Only used if `labels` is `TRUE`.
#' @param ... Further arguments to be passed to [graphics::polygon()].
#' @example inst/examples/ex-soil.R
#' @author N. Frerebeau
#' @docType methods
#' @family charts
#' @name triangle_soil
NULL

# TERNARY PHASE DIAGRAM
#' @include AllGenerics.R
NULL

# CAS phase diagram ============================================================
#' @export
#' @rdname triangle_phase_cas
triangle_phase_cas <- function(labels = TRUE, symbol = FALSE,
                               mol = FALSE, ...) {
  oxide_mass <- c(CaO  = 56.0774, Al2O3 = 101.9600, SiO2 = 60.0800)
  .triangle_phases("CAS", oxide_mass = oxide_mass, mol = mol,
                   labels = labels, symbol = symbol, ...)

  invisible(NULL)
}

#' @export
#' @rdname triangle_phase_cas
triangle_phase_ceramic <- function(labels = TRUE, symbol = FALSE,
                                   mol = FALSE, ...) {
  oxide_mass <- c(CaO  = 56.0774, Al2O3 = 101.9600, SiO2 = 60.0800)
  .triangle_phases("ceramic", oxide_mass = oxide_mass, mol = mol,
                   labels = labels, symbol = symbol, ...)

  invisible(NULL)
}

.triangle_phases <- function(chart, oxide_mass, mol = FALSE,
                             labels = TRUE, symbol = FALSE, ...) {
  ## 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")

  poly <- .phases[[chart]]

  if (!mol) {
    oxyde_mol <- as.matrix(poly[, c(1, 2, 3)])

    ## Molar mass (g/mol)
    phase_mass <- oxyde_mol %*% oxide_mass

    ## Oxide weight (%)
    poly[, c(1, 2, 3)] <- t(t(oxyde_mol) * oxide_mass) / as.vector(phase_mass)
  }

  lab <- poly[!duplicated(poly$label), ]
  txt <- if (symbol && !all(lab$symbol == "")) lab$symbol else lab$label

  poly$group <- factor(poly$group, levels = unique(poly$group))
  poly <- split(poly, f = poly$group)

  for (i in poly) {
    ternary_polygon(i, ...)
  }
  if (labels) {
    ternary_points(lab, cex = cex.lab, col = col.lab, ...)
    ternary_text(lab, label = txt, pos = lab$pos,
                 cex = cex.lab, col = col.lab, font = font.lab)
  }
}

# CONVEX HULL
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_hull
#' @aliases ternary_hull,numeric,numeric,numeric-method
setMethod(
  f = "ternary_hull",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, center = FALSE, scale = FALSE, ...) {
    coords <- coordinates_ternary(x, y, z, center = center, scale = scale)
    hull <- grDevices::chull(coords)
    graphics::polygon(x = coords$x[hull], y = coords$y[hull], ...)

    coords <- utils::modifyList(coords, list(x = x, y = y, z = z))
    invisible(coords)
  }
)

#' @export
#' @rdname ternary_hull
#' @aliases ternary_hull,ANY,missing,missing-method
setMethod(
  f = "ternary_hull",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, center = FALSE, scale = FALSE, ...) {
    xyz <- grDevices::xyz.coords(x)
    coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
                                   center = center, scale = scale, ...)
    invisible(coords)
  }
)

# PCA
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_pca
#' @aliases ternary_pca,numeric,numeric,numeric-method
setMethod(
  f = "ternary_pca",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, axis = 1, ...) {
    ## CLR
    coda <- cbind(x, y, z)
    ratio <- clr(coda)

    z <- ratio - rowMeans(ratio) # Center
    m <- colMeans(z)

    ## Get eigenvectors
    eig <- eigen(stats::cov(z))$vectors[, axis[[1L]]] + m

    ## Standard coordinates
    std <- z %*% eig

    lam <- seq(-5, 5, length.out = nrow(ratio))
    axe <- cbind(eig[1L] * lam, eig[2L] * lam, eig[3L] * lam) +
      cbind(m[1L] * (1 - lam), m[2L] * (1 - lam), m[3L] * (1 - lam))

    ## Inverse CLR
    axe <- clr_inv(axe)
    coords <- coordinates_ternary(axe)

    ## Plot
    graphics::lines(coords, ...)

    invisible(list(x = x, y = y, z = z))
  }
)

#' @export
#' @rdname ternary_pca
#' @aliases ternary_pca,ANY,missing,missing-method
setMethod(
  f = "ternary_pca",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, axis = 1, ...) {
    xyz <- grDevices::xyz.coords(x)
    coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z, axis = axis, ...)
    invisible(coords)
  }
)

# TERNARY POINTS
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_points
#' @aliases ternary_points,numeric,numeric,numeric-method
setMethod(
  f = "ternary_points",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, center = FALSE, scale = FALSE, type = "p", ...) {
    pt <- coordinates_ternary(x, y, z, center = center, scale = scale)
    graphics::points(x = pt, type = type, ...)

    pt <- utils::modifyList(pt, list(x = x, y = y, z = z))
    invisible(pt)
  }
)

#' @export
#' @rdname ternary_points
#' @aliases ternary_points,ANY,missing,missing-method
setMethod(
  f = "ternary_points",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, center = FALSE, scale = FALSE, type = "p", ...) {
    xyz <- grDevices::xyz.coords(x)
    pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
                               center = center, scale = scale,
                               type = type, ...)
    invisible(pt)
  }
)

# TERNARY CROSS-HAIRS
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_crosshairs
#' @aliases ternary_crosshairs,numeric,numeric,numeric-method
setMethod(
  f = "ternary_crosshairs",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z,
                        x_mark = TRUE, y_mark = TRUE, z_mark = TRUE, ...) {
    total <- x + y + z
    zero <- rep(0, length(x))

    if (x_mark) {
      ternary_segments(x0 = x, y0 = y, z0 = z,
                       x1 = x / total, y1 = 1 - (x / total), z1 = zero, ...)
    }
    if (y_mark) {
      ternary_segments(x0 = x, y0 = y, z0 = z,
                       x1 = zero, y1 = y / total, z1 = 1 - (y / total), ...)
    }
    if (z_mark) {
      ternary_segments(x0 = x, y0 = y, z0 = z,
                       x1 = 1 - (z / total), y1 = zero, z1 = z / total, ...)
    }

    invisible(list(x = x, y = y, z = z))
  }
)

#' @export
#' @rdname ternary_crosshairs
#' @aliases ternary_crosshairs,ANY,missing,missing-method
setMethod(
  f = "ternary_crosshairs",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, x_mark = TRUE, y_mark = TRUE, z_mark = TRUE, ...) {
    xyz <- grDevices::xyz.coords(x)
    pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
                               x_mark = x_mark, y_mark = y_mark, z_mark = z_mark, ...)
    invisible(pt)
  }
)

# TERNARY ARROWS
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_arrows
#' @aliases ternary_arrows,numeric,numeric,numeric-method
setMethod(
  f = "ternary_arrows",
  signature = c(x0 = "numeric", y0 = "numeric", z0 = "numeric"),
  definition = function(x0, y0, z0, x1 = x0, y1 = y0, z1 = z0, ...) {
    coords0 <- coordinates_ternary(x0, y0, z0)
    coords1 <- coordinates_ternary(x1, y1, z1)
    graphics::arrows(x0 = coords0$x, y0 = coords0$y,
                     x1 = coords1$x, y1 = coords1$y, ...)

    invisible(NULL)
  }
)

# TERNARY PLOT
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_plot
#' @aliases ternary_plot,numeric,numeric,numeric-method
setMethod(
  f = "ternary_plot",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, center = FALSE, scale = FALSE,
                        xlim = NULL, ylim = NULL, zlim = NULL,
                        xlab = NULL, ylab = NULL, zlab = NULL,
                        main = NULL, sub = NULL, ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = axes,
                        panel.first = NULL, panel.last = NULL, ...) {

    ## Save and restore graphical parameters
    ## pty: square plotting region, independent of device size
    old_par <- graphics::par(pty = "s", no.readonly = TRUE)
    on.exit(graphics::par(old_par), add = TRUE)

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

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

    ## Set plotting coordinates
    ternary_window(xlim = xlim, ylim = ylim, zlim = zlim,
                   xlab = xlab, ylab = ylab, zlab = zlab,
                   cex = cex.lab)

    ## Reset center and scale
    assign("center", value = NULL, envir = the)
    assign("scale", value = NULL, envir = the)

    ## Compute ternary coordinates
    pt <- coordinates_ternary(x = x, y = y, z = z, center = center, scale = scale)

    ## Save center and scale for further use, e.g. grid or axes.
    assign("center", value = pt$center, envir = the)
    assign("scale", value = pt$scale, envir = the)

    ## Evaluate pre-plot expressions
    panel.first

    ## Plot
    graphics::points(x = pt, ...)

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

    ## Construct Axis
    if (axes) {
      ternary_axis(side = 1, center = pt$center, scale = pt$scale, col = fg)
      ternary_axis(side = 2, center = pt$center, scale = pt$scale, col = fg)
      ternary_axis(side = 3, center = pt$center, scale = pt$scale, col = fg)
    }

    ## Plot frame
    if (frame.plot) {
      ternary_box(lty = "solid", lwd = 1, col = fg)
    }

    ## Add annotation
    if (ann) {
      ternary_title(main = main, sub = sub, xlab = xlab, ylab = ylab,
                    zlab = zlab, ...)
    }

    invisible(pt)
  }
)

#' @export
#' @rdname ternary_plot
#' @aliases ternary_plot,ANY,missing,missing-method
setMethod(
  f = "ternary_plot",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, xlim = NULL, ylim = NULL, zlim = NULL,
                        xlab = NULL, ylab = NULL, zlab = NULL,
                        main = NULL, sub = NULL, ann = graphics::par("ann"),
                        axes = TRUE, frame.plot = axes,
                        panel.first = NULL, panel.last = NULL, ...) {

    xyz <- grDevices::xyz.coords(x, xlab = xlab, ylab = ylab, zlab = zlab)
    pt <- methods::callGeneric(
      x = xyz$x, y = xyz$y, z = xyz$z,
      xlim = xlim,
      ylim = ylim,
      zlim = zlim,
      xlab = xlab %||% xyz$xlab %||% "x",
      ylab = ylab %||% xyz$ylab %||% "y",
      zlab = zlab %||% xyz$zlab %||% "z",
      main = main, sub = sub, ann = ann,
      axes = axes,
      frame.plot = frame.plot,
      panel.first = panel.first,
      panel.last = panel.last,
      ...
    )

    invisible(pt)
  }
)

# TERNARY BOX
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_box
ternary_box <- function(lty = "solid", ...) {

  ## Graphical parameters
  col <- list(...)$col %||% graphics::par("fg")
  lwd <- list(...)$lwd %||% 1

  graphics::polygon(x = c(0, 0.5, 1), y = c(0, the$top, 0),
                    border = col, lty = lty, lwd = lwd)

  invisible(NULL)
}

# TERNARY TEXT
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_text
#' @aliases ternary_text,numeric,numeric,numeric-method
setMethod(
  f = "ternary_text",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, center = FALSE, scale = FALSE,
                        labels = seq_along(x), ...) {
    coords <- coordinates_ternary(x, y, z, center = center, scale = scale)
    graphics::text(x = coords, labels = labels, ...)

    invisible(list(x = x, y = y, z = z))
  }
)

#' @export
#' @rdname ternary_text
#' @aliases ternary_text,ANY,missing,missing-method
setMethod(
  f = "ternary_text",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, center = FALSE, scale = FALSE, labels = seq_along(x$x), ...) {
    x <- grDevices::xyz.coords(x)
    force(labels)
    coords <- methods::callGeneric(x = x$x, y = x$y, z = x$z,
                                   center = center, scale = scale,
                                   labels = labels, ...)
    invisible(coords)
  }
)

# TERNARY SEGMENTS
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_segments
#' @aliases ternary_segments,numeric,numeric,numeric-method
setMethod(
  f = "ternary_segments",
  signature = c(x0 = "numeric", y0 = "numeric", z0 = "numeric"),
  definition = function(x0, y0, z0, x1 = x0, y1 = y0, z1 = z0, ...) {
    coords0 <- coordinates_ternary(x0, y0, z0)
    coords1 <- coordinates_ternary(x1, y1, z1)
    graphics::segments(x0 = coords0$x, y0 = coords0$y,
                       x1 = coords1$x, y1 = coords1$y, ...)

    invisible(NULL)
  }
)

# TERNARY POLYGON
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_polygon
#' @aliases ternary_polygon,numeric,numeric,numeric-method
setMethod(
  f = "ternary_polygon",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, ...) {
    coords <- coordinates_ternary(x, y, z)
    graphics::polygon(x = coords, ...)

    invisible(list(x = x, y = y, z = z))
  }
)

#' @export
#' @rdname ternary_polygon
#' @aliases ternary_polygon,ANY,missing,missing-method
setMethod(
  f = "ternary_polygon",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, ...) {
    xyz <- grDevices::xyz.coords(x)
    coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z, ...)
    invisible(coords)
  }
)

# TERNARY LINES
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ternary_lines
#' @aliases ternary_lines,numeric,numeric,numeric-method
setMethod(
  f = "ternary_lines",
  signature = c(x = "numeric", y = "numeric", z = "numeric"),
  definition = function(x, y, z, type = "l", ...) {
    coords <- coordinates_ternary(x, y, z)
    graphics::lines(x = coords, type = type, ...)

    invisible(list(x = x, y = y, z = z))
  }
)

#' @export
#' @rdname ternary_lines
#' @aliases ternary_lines,ANY,missing,missing-method
setMethod(
  f = "ternary_lines",
  signature = c(x = "ANY", y = "missing", z = "missing"),
  definition = function(x, type = "l", ...) {
    xyz <- grDevices::xyz.coords(x)
    coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
                                   type = type, ...)
    invisible(coords)
  }
)

# INTERNAL ENVIRONMENT

the <- new.env(parent = emptyenv())
the$top <- sqrt(3) / 2
the$center <- NULL # Updated by plot()
the$scale <- NULL  # Updated by plot()

get_center <- function(...) {
  get("center", envir = the)
}

get_scale <- function(...) {
  get("scale", envir = the)
}

1		# TERNARY TITLE
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_title
7		ternary_title <- function(main = NULL, sub = NULL, xlab = NULL, ylab = NULL,
8		zlab = NULL, line = NA, outer = FALSE, ...) {
9		## Graphical parameters
10	86x	cex.lab <- list(...)$cex.lab %\|\|% graphics::par("cex.lab")
11	86x	col.lab <- list(...)$col.lab %\|\|% graphics::par("col.lab")
12	86x	font.lab <- list(...)$font.lab %\|\|% graphics::par("font.lab")
13
14		## Axes labels
15	86x	xlab <- grDevices::as.graphicsAnnot(xlab)
16	86x	ylab <- grDevices::as.graphicsAnnot(ylab)
17	86x	zlab <- grDevices::as.graphicsAnnot(zlab)
18
19	86x	if (!is.null(xlab)) {
20	86x	graphics::text(x = 0, y = 0, label = xlab, pos = 1,
21	86x	col = col.lab, cex = cex.lab, font = font.lab)
22		}
23	86x	if (!is.null(ylab)) {
24	86x	graphics::text(x = 1, y = 0, label = ylab, pos = 1,
25	86x	col = col.lab, cex = cex.lab, font = font.lab)
26		}
27	86x	if (!is.null(zlab)) {
28	86x	graphics::text(x = 0.5, y = the$top, label = zlab, pos = 3,
29	86x	col = col.lab, cex = cex.lab, font = font.lab)
30		}
31
32		## Title
33	86x	graphics::title(main = main, sub = sub, xlab = NULL, ylab = NULL,
34	86x	line = line, outer = outer)
35
36	86x	invisible(NULL)
37		}

1		# CONTOUR
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_contour
7		#' @aliases ternary_contour,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_contour",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, value, n = 50, nlevels = 10,
12		levels = pretty(range(value, na.rm = TRUE), nlevels),
13		ilr = TRUE, method = "linear", extrapolate = FALSE,
14		palette = function(i) grDevices::hcl.colors(i, "YlOrRd", rev = TRUE),
15		...) {
16		## Calculate contour lines
17	2x	xy <- coordinates_contour(x = x, y = y, z = z, value = value, n = n,
18	2x	nlevels = nlevels, levels = levels,
19	2x	ilr = ilr, method = method,
20	2x	extrapolate = extrapolate)
21
22		## Get contour levels
23	2x	lvl <- vapply(X = xy, FUN = getElement, FUN.VALUE = numeric(1),
24	2x	name = "level")
25
26		## Colors
27		## (number of levels may differ from nlevels due to interp())
28	2x	col <- palette(length(unique(lvl)))
29	2x	names(col) <- unique(lvl)
30	2x	col <- col[as.character(lvl)]
31
32		## Plot
33	2x	for (i in seq_along(xy)) {
34		## Get contour
35	133x	level <- xy[[i]]
36
37		## Inverse ILR transform
38	133x	tern <- cbind(level$x, level$y)
39	133x	tern <- if (ilr) ilr_inv(tern) else coordinates_cartesian(tern)
40
41		## Plot ternary lines
42	133x	ternary_lines(tern, col = col[[i]], ...)
43		}
44
45	2x	invisible(list(levels = lvl, colors = col))
46		}
47		)
48
49		#' @export
50		#' @rdname ternary_contour
51		#' @aliases ternary_contour,ANY,missing,missing-method
52		setMethod(
53		f = "ternary_contour",
54		signature = c(x = "ANY", y = "missing", z = "missing"),
55		definition = function(x, value, n = 50, nlevels = 10,
56		levels = pretty(range(value, na.rm = TRUE), nlevels),
57		ilr = TRUE, method = "linear", extrapolate = FALSE,
58		palette = function(i) grDevices::hcl.colors(i, "YlOrRd", rev = TRUE),
59		...) {
60	2x	xyz <- grDevices::xyz.coords(x)
61	2x	coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z, value = value,
62	2x	n = n, nlevels = nlevels, levels = levels,
63	2x	ilr = ilr, method = method,
64	2x	extrapolate = extrapolate,
65	2x	palette = palette, ...)
66	2x	invisible(coords)
67		}
68		)
69
70		#' Calculate Contour Lines
71		#'
72		#' Computes contours coordinates.
73		#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
74		#' of a set of points. If `y` and `z` are missing, an attempt is made to
75		#' interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
76		#' @param value A [`numeric`] [`matrix`] containing the values to be plotted.
77		#' @param n A length-one [`numeric`] specifying the number of grid points.
78		#' @param nlevels A length-one [`numeric`] vector specifying the number of
79		#' contour levels desired. Only used if `levels` is `NULL`.
80		#' @param levels A [`numeric`] vector of levels at which to draw contour lines.
81		#' @param ilr A [`logical`] scalar: should interpolation be computed in ILR
82		#' space? If `FALSE`, interpolation is computed in Cartesian space.
83		#' @param method A [`character`] string: specifying the method for interpolation
84		#' (see [interp::interp()]).
85		#' @param extrapolate A [`logical`] scalar: should extrapolation be used outside
86		#' of the convex hull determined by the data points (see [interp::interp()])?
87		#' @param ... Further parameters to be passed to [interp::interp()].
88		#' @return
89		#' A [`list`] of contours, each itself a list with elements
90		#' (see [grDevices::contourLines()]):
91		#' \tabular{ll}{
92		#' `level` \tab The contour level. \cr
93		#' `x` \tab The (ILR) x-coordinates of the contour. \cr
94		#' `y` \tab The (ILR) y-coordinates of the contour. \cr
95		#' }
96		#' @keywords internal
97		#' @noRd
98		coordinates_contour <- function(x, y, z, value, n = 50, nlevels = 10,
99		levels = pretty(range(value, na.rm = TRUE), nlevels),
100		ilr = TRUE, method = "linear", extrapolate = FALSE,
101		...) {
102		## Validation
103	2x	assert_package("interp")
104	2x	assert_length(value, length(x))
105
106		## ILR vs Cartesian
107	2x	coda <- cbind(x, y, z)
108	2x	ratio <- if (ilr) ilr(coda) else do.call(cbind, coordinates_ternary(coda))
109
110		## Remove NA/Inf (if any)
111	2x	ok <- apply(X = ratio, MARGIN = 1, FUN = function(x) all(is.finite(x)))
112	2x	ratio <- ratio[ok, , drop = FALSE]
113	2x	value <- value[ok]
114
115		## Interpolate
116	2x	xlim <- expand_range(ratio[, 1], mult = 0.2)
117	2x	ylim <- expand_range(ratio[, 2], mult = 0.2)
118	2x	interp <- interp::interp(
119	2x	x = ratio[, 1],
120	2x	y = ratio[, 2],
121	2x	z = value,
122	2x	xo = seq(xlim[1L], xlim[2L], length.out = n),
123	2x	yo = seq(ylim[1L], ylim[2L], length.out = n),
124	2x	method = method,
125	2x	extrap = extrapolate,
126		...
127		)
128
129		## Compute contours
130	2x	grDevices::contourLines(
131	2x	x = interp$x,
132	2x	y = interp$y,
133	2x	z = interp$z,
134	2x	nlevels = nlevels,
135	2x	levels = levels
136		)
137		}

1		# DENSITY
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_density
7		#' @aliases ternary_density,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_density",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, h = NULL, n = 25, nlevels = 10, levels = NULL,
12		palette = function(i) grDevices::hcl.colors(i, "YlOrRd", rev = TRUE),
13		...) {
14		## Calculate density contour lines
15	1x	xy <- coordinates_kde(x = x, y = y, z = z, h = h, n = n,
16	1x	nlevels = nlevels, levels = levels)
17
18		## Get contour levels
19	1x	lvl <- vapply(X = xy, FUN = getElement, FUN.VALUE = numeric(1),
20	1x	name = "level")
21
22		## Colors
23		## (number of levels may differ from nlevels due to pretty())
24	1x	col <- palette(length(unique(lvl)))
25	1x	names(col) <- unique(lvl)
26	1x	col <- col[as.character(lvl)]
27
28		## Plot
29	1x	for (i in seq_along(xy)) {
30		## Get contour
31	9x	level <- xy[[i]]
32
33		## Inverse ILR transform
34	9x	tern <- ilr_inv(cbind(level$x, level$y))
35
36		## Plot ternary lines
37	9x	ternary_lines(tern, col = col[[i]], ...)
38		}
39
40	1x	invisible(list(levels = lvl, colors = col))
41		}
42		)
43
44		#' @export
45		#' @rdname ternary_density
46		#' @aliases ternary_density,ANY,missing,missing-method
47		setMethod(
48		f = "ternary_density",
49		signature = c(x = "ANY", y = "missing", z = "missing"),
50		definition = function(x, h = NULL, n = 25, nlevels = 10, levels = NULL,
51		palette = function(i) grDevices::hcl.colors(i, "YlOrRd", rev = TRUE),
52		...) {
53	1x	xyz <- grDevices::xyz.coords(x)
54	1x	pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
55	1x	h = h, n = n, nlevels = nlevels, levels = levels,
56	1x	palette = palette, ...)
57	1x	invisible(pt)
58		}
59		)
60
61		#' Calculate KDE Contour Lines
62		#'
63		#' Computes 2D-KDE contours coordinates.
64		#' @param x,y,z A [`numeric`] vector giving the x, y and z ternary coordinates
65		#' of a set of points. If `y` and `z` are missing, an attempt is made to
66		#' interpret `x` in a suitable way (see [grDevices::xyz.coords()]).
67		#' @param h A length-one [`numeric`] vector giving the bandwidth.
68		#' @param n A length-one [`numeric`] specifying the number of grid points.
69		#' @param nlevels A length-one [`numeric`] vector specifying the number of
70		#' contour levels desired. Only used if `levels` is `NULL`.
71		#' @param levels A [`numeric`] vector of levels at which to draw contour lines.
72		#' @return
73		#' A [`list`] of contours, each itself a list with elements
74		#' (see [grDevices::contourLines()]):
75		#' \tabular{ll}{
76		#' `level` \tab The contour level. \cr
77		#' `x` \tab The ILR x-coordinates of the contour. \cr
78		#' `y` \tab The ILR y-coordinates of the contour. \cr
79		#' }
80		#' @keywords internal
81		#' @noRd
82		coordinates_kde <- function(x, y, z, h = NULL, n = 25,
83		nlevels = 10, levels = NULL) {
84		## ILR
85	1x	coda <- cbind(x, y, z)
86	1x	ratio <- ilr(coda)
87
88		## Compute KDE
89	1x	lims <- expand_range(ratio, mult = 0.2)
90	1x	dens <- kde(
91	1x	x = ratio[, 1],
92	1x	y = ratio[, 2],
93	1x	h = h,
94	1x	n = n,
95	1x	xlim = lims, # x and y range should be same
96	1x	ylim = lims
97		)
98
99		## Compute contours
100	1x	grDevices::contourLines(
101	1x	x = dens$x,
102	1x	y = dens$y,
103	1x	z = dens$z,
104	1x	nlevels = nlevels,
105	1x	levels = levels %\|\|% pretty(range(dens$z, na.rm = TRUE), nlevels)
106		)
107		}
108
109		## Adapted from MASS::kde2d
110		kde <- function(x, y, h = NULL, n = 25, gx = NULL, gy = NULL,
111		xlim = range(x), ylim = range(y)) {
112	2x	n <- rep(n, length.out = 2L)
113	1x	if (is.null(gx)) gx <- seq(xlim[1L], xlim[2L], length.out = n[1L])
114	1x	if (is.null(gy)) gy <- seq(ylim[1L], ylim[2L], length.out = n[2L])
115
116	2x	if (is.null(h)) {
117	2x	h <- c(bandwidth(x), bandwidth(y))
118		} else {
119	!	h <- rep(h, length.out = 2L)
120		}
121	2x	h <- h / 4
122
123	2x	if (any(h <= 0))
124	!	stop("Bandwidths must be strictly positive.", call. = FALSE)
125
126	2x	ax <- outer(gx, x, "-") / h[1L]
127	2x	ay <- outer(gy, y, "-") / h[2L]
128
129	2x	nx <- length(x)
130	2x	z <- tcrossprod(
131	2x	matrix(stats::dnorm(ax), ncol = nx),
132	2x	matrix(stats::dnorm(ay), ncol = nx)
133		)
134	2x	z <- z / (nx * h[1L] * h[2L])
135
136	2x	list(x = gx, y = gy, z = z)
137		}
138
139		# Copied from MASS::bandwidth.nrd()
140		bandwidth <- function(x) {
141	4x	r <- stats::quantile(x, c(0.25, 0.75))
142	4x	h <- (r[2L] - r[1L]) / 1.34
143	4x	4 * 1.06 * min(sqrt(stats::var(x)), h) * length(x)^(-1 / 5)
144		}

1		# HELPERS
2
3		## https://michaelchirico.github.io/potools/articles/developers.html
4		tr_ <- function(...) {
5	9x	enc2utf8(gettext(paste0(...), domain = "R-isopleuros"))
6		}
7
8		`%\|\|%` <- function(x, y) {
9	441x	if (!is.null(x)) x else y
10		}
11
12		map_color <- function(values, palette = NULL, scheme = "viridis",
13		ignore_zero = FALSE) {
14	1x	if (isFALSE(palette)) return(values)
15
16	3x	if (is.null(palette)) {
17	3x	palette <- function(x) {
18	3x	x <- (x - min(x)) / (max(x) - min(x)) # Rescale to [0,1]
19	3x	col <- grDevices::hcl.colors(256L, palette = scheme)
20	3x	grDevices::rgb(grDevices::colorRamp(col)(x), maxColorValue = 255)
21		}
22		}
23
24	3x	color <- rep(NA, length(values))
25	3x	ok <- is.finite(values) # Remove NA/Inf (if any)
26	!	if (ignore_zero) ok[ok] <- values[ok] > 0
27	3x	color[ok] <- palette(values[ok])
28	3x	color
29		}
30
31		#' Expand Range
32		#'
33		#' @param x A [`numeric`] vector.
34		#' @param mult A [`numeric`] value giving the multiplicative constant.
35		#' @param add A [`numeric`] value giving the additive constant.
36		#' @return A length-two [`numeric`] vector.
37		#' @keywords internal
38		#' @noRd
39		expand_range <- function(x, mult = 0, add = 0) {
40	92x	lims <- range(x)
41	92x	lims <- lims + c(-1, 1) * (diff(lims) * mult + add)
42	92x	lims
43		}
44
45		#' Rotate Around a Point
46		#'
47		#' @param x A column vector giving the x and y coordinates of the point to be
48		#' rotated.
49		#' @param theta A length-one [`numeric`] vector specifying the rotation angle
50		#' (in radian).
51		#' @param origin A length-two [`numeric`] vector specifying the coordinates
52		#' of the point to rotate around.
53		#' @return A `matrix` of coordinates.
54		#' @keywords internal
55		#' @noRd
56		rotate <- function(x, theta = 0, origin = c(0.5, sqrt(3) / 6)) {
57		## Translation matrix
58	138x	trans <- diag(1, 3, 3)
59	138x	trans[, 3] <- c(origin, 1)
60
61		## Rotation matrix
62	138x	rot <- matrix(
63	138x	data = c(cos(theta), sin(theta), 0, -sin(theta), cos(theta), 0, 0, 0, 1),
64	138x	nrow = 3,
65	138x	ncol = 3
66		)
67
68	138x	x <- as.matrix(x)
69	138x	if (nrow(x) < 3) x <- rbind(x, rep(1, ncol(x)))
70	138x	t(trans %% rot %% solve(trans) %*% x)
71		}
72
73		#' Check Object Length
74		#'
75		#' @param x An object to be checked.
76		#' @param expected An appropriate expected value.
77		#' @return
78		#' Throws an error, if any, and returns `x` invisibly otherwise.
79		#' @keywords internal
80		#' @noRd
81		assert_length <- function(x, expected) {
82	2763x	arg <- deparse(substitute(x))
83	2763x	if (length(x) != expected) {
84	3x	str <- tr_("%s must be of length %d; not %d.")
85	3x	msg <- sprintf(str, sQuote(arg), expected, length(x))
86	3x	stop(msg, call. = FALSE)
87		}
88	2760x	invisible(x)
89		}
90
91		assert_center <- function(x, current = get_center()) {
92	884x	centered <- is.numeric(current) && !all(current == 1)
93	884x	ok <- isTRUE(x) \|\| is.numeric(x)
94	884x	if (!ok && centered) {
95	2x	msg <- tr_("The current plot has been centered, but your data doesn't seem to be.")
96	2x	message(msg)
97		}
98	884x	invisible(x)
99		}
100
101		assert_scale <- function(x, current = get_scale()) {
102	884x	scaled <- is.numeric(current) && !all(current == 1)
103	884x	ok <- isTRUE(x) \|\| is.numeric(x)
104	884x	if (!ok && scaled) {
105	2x	msg <- tr_("The current plot has been scaled, but your data doesn't seem to be.")
106	2x	message(msg)
107		}
108	884x	invisible(x)
109		}
110
111		assert_package <- function(x) {
112	3x	if (!requireNamespace(x, quietly = TRUE)) {
113	!	msg <- tr_("Package %s needed for this function to work. Please install it.")
114	!	stop(sprintf(msg, sQuote(x)), call. = FALSE)
115		}
116	3x	invisible(NULL)
117		}

1		# PAIRS
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_pairs
7		#' @aliases ternary_pairs,matrix-method
8		setMethod(
9		f = "ternary_pairs",
10		signature = c(x = "matrix"),
11		definition = function(x, margin = NULL, ...) {
12
13		## Save and restore graphical parameters
14	3x	old_par <- graphics::par(no.readonly = TRUE)
15	3x	on.exit(graphics::par(old_par), add = TRUE)
16
17		## Layout
18	3x	n <- ncol(x)
19	3x	ni <- seq_len(n)
20	3x	parts <- colnames(x) %\|\|% paste0("X", ni)
21	3x	zlab <- "*"
22
23	3x	if (!is.null(margin)) {
24	1x	margin <- if (is.character(margin)) which(parts == margin) else as.integer(margin)
25	1x	star <- x[, margin]
26	1x	zlab <- parts[margin]
27	1x	parts <- parts[-margin]
28	1x	n <- n - 1
29	1x	ni <- ni[-margin]
30		}
31
32	3x	k <- (n^2 - n) / 2
33	3x	lay <- matrix(0, nrow = n, ncol = n)
34	3x	lay[lower.tri(lay, diag = FALSE)] <- seq_len(k)
35	3x	diag(lay) <- seq(k + 1, k + n)
36	3x	lay <- t(lay)
37
38	3x	graphics::layout(lay)
39	3x	graphics::par(mar = c(0, 0, 0, 0) + 0.1, oma = c(1, 1, 1, 1))
40
41		## Ternary plots
42	3x	p <- utils::combn(ni, 2, simplify = FALSE)
43	3x	for (i in p) {
44	39x	if (is.null(margin)) {
45	24x	star <- apply(X = x[, -i, drop = FALSE], MARGIN = 1, FUN = gmean)
46		}
47	39x	z <- cbind(x[, i], star)
48	39x	ternary_plot(z, zlab = zlab, axes = FALSE, frame.plot = TRUE, ...)
49		}
50
51		## Graphical parameters
52	3x	str.wid <- max(graphics::strwidth(parts, "user"))
53	3x	cex.lab <- list(...)$cex.lab %\|\|% max(0.8, min(2, 0.9 / str.wid))
54	3x	col.lab <- list(...)$col.lab %\|\|% graphics::par("col.lab")
55	3x	font.lab <- list(...)$font.lab %\|\|% graphics::par("font.lab")
56
57		## Add labels
58	3x	for (part in parts) {
59	16x	graphics::plot(NULL, xlim = c(0, 1), ylim = c(0, 1),
60	16x	axes = FALSE, ann = FALSE)
61	16x	graphics::text(x = 0.5, y = 0.5, labels = part,
62	16x	cex = cex.lab, col = col.lab, font = font.lab)
63		}
64
65	3x	invisible(x)
66		}
67		)
68
69		#' @export
70		#' @rdname ternary_pairs
71		#' @aliases ternary_pairs,data.frame,missing,missing-method
72		setMethod(
73		f = "ternary_pairs",
74		signature = c(x = "data.frame"),
75		definition = function(x, margin = NULL, ...) {
76	3x	x <- data.matrix(x)
77	3x	methods::callGeneric(x = x, margin = margin, ...)
78
79	3x	invisible(x)
80		}
81		)

1		# TERNARY COORDINATES
2		#' @include AllGenerics.R
3		NULL
4
5		# Ternary to cartesian =========================================================
6		#' @export
7		#' @rdname coordinates_ternary
8		#' @aliases coordinates_ternary,numeric,numeric,numeric-method
9		setMethod(
10		f = "coordinates_ternary",
11		signature = c(x = "numeric", y = "numeric", z = "numeric"),
12		definition = function(x, y, z, center = FALSE, scale = FALSE,
13		missing = getOption("isopleuros.missing")) {
14		## Validation
15	886x	n <- length(x)
16	886x	assert_length(y, n)
17	885x	assert_length(z, n)
18	884x	assert_center(center)
19	884x	assert_scale(scale)
20
21		## Missing values
22	884x	if (missing) {
23	!	x[is.na(x)] <- 0
24	!	y[is.na(y)] <- 0
25	!	z[is.na(z)] <- 0
26		}
27
28	884x	total <- x + y + z
29	884x	na <- is.na(x) \| is.na(y) \| is.na(z)
30	884x	zero <- total == 0
31
32	884x	x <- x[!na & !zero]
33	884x	y <- y[!na & !zero]
34	884x	z <- z[!na & !zero]
35	884x	total <- total[!na & !zero]
36
37		## Validation
38	884x	if (any(x < 0 \| y < 0 \| z < 0)) {
39	1x	stop(tr_("Positive values are expected."), call. = FALSE)
40		}
41
42	883x	coord <- matrix(data = c(x, y, z), ncol = 3) / total
43	883x	coord <- scale(coord, center = center, scale = scale)
44
45	883x	list(
46	883x	x = coord$y + coord$z / 2,
47	883x	y = coord$z * sqrt(3) / 2,
48	883x	center = coord$center,
49	883x	scale = coord$scale
50		)
51		}
52		)
53
54		#' @export
55		#' @rdname coordinates_ternary
56		#' @aliases coordinates_ternary,ANY,missing,missing-method
57		setMethod(
58		f = "coordinates_ternary",
59		signature = c(x = "ANY", y = "missing", z = "missing"),
60		definition = function(x, xlab = NULL, ylab = NULL, zlab = NULL,
61		center = FALSE, scale = FALSE,
62		missing = getOption("isopleuros.missing")) {
63	504x	xyz <- grDevices::xyz.coords(x, xlab = xlab, ylab = ylab, zlab = zlab)
64	504x	methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
65	504x	center = center, scale = scale, missing = missing)
66		}
67		)
68
69		# Cartesian to ternary =========================================================
70		#' @export
71		#' @rdname coordinates_cartesian
72		#' @aliases coordinates_cartesian,numeric,numeric-method
73		setMethod(
74		f = "coordinates_cartesian",
75		signature = c(x = "numeric", y = "numeric"),
76		definition = function(x, y) {
77		## Validation
78	652x	n <- length(x)
79	652x	assert_length(y, n)
80
81	651x	k <- y * 2 / sqrt(3)
82	651x	j <- x - k / 2
83	651x	i <- 1 - (j + k)
84
85	651x	list(
86	651x	x = i,
87	651x	y = j,
88	651x	z = k
89		)
90		}
91		)
92
93		#' @export
94		#' @rdname coordinates_cartesian
95		#' @aliases coordinates_cartesian,ANY,missing-method
96		setMethod(
97		f = "coordinates_cartesian",
98		signature = c(x = "ANY", y = "missing"),
99		definition = function(x, xlab = NULL, ylab = NULL) {
100	47x	xy <- grDevices::xy.coords(x, xlab = xlab, ylab = ylab)
101	47x	methods::callGeneric(x = xy$x, y = xy$y)
102		}
103		)
104
105		# Scale ========================================================================
106		#' Center and Scale
107		#'
108		#' @param x,y,z The x, y and z coordinates of a set of points. Both y and z can
109		#' be left at `NULL`. In this case, an attempt is made to interpret x in a way
110		#' suitable for plotting.
111		#' @param center A [`logical`] scalar or a [`numeric`] vector giving the center.
112		#' @param scale A [`logical`] scalar or a length-one [`numeric`] vector giving a
113		#' scaling factor.
114		#' @return
115		#' A [`list`] with the components:
116		#' \tabular{ll}{
117		#' `x` \tab A [`numeric`] vector of x values. \cr
118		#' `y` \tab A [`numeric`] vector of y values. \cr
119		#' `z` \tab A [`numeric`] vector of z values. \cr
120		#' `center` \tab A [`numeric`] vector giving the center. \cr
121		#' `scale` \tab A [`numeric`] vector giving the scale factor. \cr
122		#' }
123		#' @keywords internal
124		#' @noRd
125		scale <- function(x, y = NULL, z = NULL, center = TRUE, scale = TRUE) {
126	883x	xyz <- grDevices::xyz.coords(x = x, y = y, z = z)
127	883x	xyz <- matrix(data = c(xyz$x, xyz$y, xyz$z), ncol = 3)
128
129	883x	y <- xyz
130	883x	if (!isFALSE(center) && !is.null(center)) {
131	162x	if (isTRUE(center)) {
132	4x	center <- apply(X = xyz, MARGIN = 2, FUN = gmean, simplify = TRUE)
133	4x	center <- center / sum(center)
134		}
135	162x	assert_length(center, NCOL(xyz))
136
137	162x	y <- t(t(y) / center)
138	162x	y <- y / rowSums(y)
139		} else {
140	721x	center <- rep(1, NCOL(xyz))
141		}
142
143	883x	if (!isFALSE(scale) && !is.null(scale)) {
144	166x	if (isTRUE(scale)) {
145	4x	scale <- sqrt(mean(diag(stats::cov(clr(xyz)))))
146		}
147	166x	assert_length(scale, 1)
148
149	166x	y <- y^(1 / scale)
150	166x	y <- y / rowSums(y)
151		} else {
152	717x	scale <- 1
153		}
154
155	883x	list(
156	883x	x = y[, 1],
157	883x	y = y[, 2],
158	883x	z = y[, 3],
159	883x	center = center,
160	883x	scale = scale
161		)
162		}
163
164		#' Geometric Mean
165		#'
166		#' @param x A [`numeric`] vector.
167		#' @param trim A length-one [`numeric`] vector specifying the fraction (0 to 0.5)
168		#' of observations to be trimmed from each end of `x` before the mean is
169		#' computed.
170		#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
171		#' computation proceeds?
172		#' @return A [`numeric`] vector.
173		#' @keywords internal
174		#' @noRd
175		gmean <- function(x, trim = 0, na.rm = FALSE) {
176		if (na.rm) x <- x[is.finite(x)]
177		x <- x[x > 0]
178		exp(mean(log(x), trim = trim))
179		}
180
181		# Centered Log-Ratios ==========================================================
182		#' Centered Log-Ratios (CLR)
183		#'
184		#' Computes CLR transformation.
185		#' @param x A [`numeric`] `matrix`.
186		#' @keywords internal
187		#' @noRd
188		clr <- function(x) {
189	6x	J <- ncol(x)
190	6x	clr <- log(x, base = exp(1)) %*% diag(J)
191
192	6x	clr
193		}
194
195		#' Inverse Centered Log-Ratios Transformation
196		#'
197		#' Computes inverse CLR transformation.
198		#' @param x A [`numeric`] `matrix` of log ratios.
199		#' @keywords internal
200		#' @noRd
201		clr_inv <- function(x) {
202	2x	y <- exp(x)
203	2x	y <- y / rowSums(y)
204
205	2x	y
206		}
207
208		# Isometric Log-Ratios =========================================================
209		#' Isometric Log-Ratios (ILR)
210		#'
211		#' Computes ILR transformation.
212		#' @param x A [`numeric`] `matrix`.
213		#' @keywords internal
214		#' @noRd
215		ilr <- function(x) {
216	10x	D <- ncol(x)
217	10x	H <- ilr_base(D)
218
219		## Rotated and centered values
220	10x	y <- log(x, base = exp(1))
221	10x	ilr <- y %*% H
222
223	10x	ilr
224		}
225
226		#' Canonical Basis for Isometric Log-Ratio transformation
227		#'
228		#' Computes the canonical basis in the CLR plane used for ILR transformation.
229		#' @param a A [`numeric`] value giving the number of parts of the simplex.
230		#' @keywords internal
231		#' @noRd
232		ilr_base <- function(n) {
233	112x	seq_parts <- seq_len(n - 1)
234
235		## Helmert matrix (rotation matrix)
236	112x	H <- stats::contr.helmert(n) # n x n-1
237	112x	H <- t(H) / sqrt((seq_parts + 1) * seq_parts) # n-1 x n
238
239		## Center
240	112x	m <- diag(x = 1, nrow = n) - matrix(data = 1 / n, nrow = n, ncol = n)
241	112x	H <- tcrossprod(m, H)
242
243	112x	H
244		}
245
246		#' Inverse Isometric Log-Ratio Transformation
247		#'
248		#' Computes inverse ILR transformation.
249		#' @param x A [`numeric`] `matrix` of log ratios.
250		#' @keywords internal
251		#' @noRd
252		ilr_inv <- function(x) {
253	102x	D <- ncol(x) + 1
254	102x	H <- ilr_base(D)
255
256	102x	y <- tcrossprod(x, H)
257	102x	y <- exp(y)
258	102x	y <- y / rowSums(y)
259
260	102x	y
261		}

1		# TERNARY LABELS
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_labels
7		#' @aliases ternary_labels,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_labels",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, center = FALSE, scale = FALSE,
12		labels = seq_along(x), type = c("text", "shadow"), ...) {
13		## Validation
14	!	type <- match.arg(type, several.ok = FALSE)
15
16		## Compute label positions
17	!	coords <- coordinates_ternary(x, y, z, center = center, scale = scale)
18	!	labs <- compute_labels(x = coords$x, y = coords$y, labels = labels, ...)
19
20		## Draw labels
21	!	fun <- switch(
22	!	type,
23	!	text = graphics::text,
24	!	shadow = text_shadow
25		)
26	!	fun(labs, labels = labels, ...)
27
28	!	coords <- utils::modifyList(coords, list(x = x, y = y, z = z))
29	!	invisible(coords)
30		}
31		)
32
33		#' @export
34		#' @rdname ternary_labels
35		#' @aliases ternary_labels,ANY,missing,missing-method
36		setMethod(
37		f = "ternary_labels",
38		signature = c(x = "ANY", y = "missing", z = "missing"),
39		definition = function(x, center = FALSE, scale = FALSE, labels = seq_along(x$x), ...) {
40	!	x <- grDevices::xyz.coords(x)
41	!	coords <- methods::callGeneric(x = x$x, y = x$y, z = x$z,
42	!	center = center, scale = scale,
43	!	labels = labels, ...)
44	!	invisible(coords)
45		}
46		)
47
48
49		#' Compute Label Positions
50		#'
51		#' @return A [`list`] with elements `x`, `y` and `labels`.
52		#' @source
53		#' This function is modeled after [car::pointLabel()] (originally from the
54		#' \pkg{maptools} package).
55		#' @keywords internal
56		#' @noRd
57		compute_labels <- function(x, y, labels, ..., iter = 50,
58		cex = graphics::par("cex"),
59		font = NULL, vfont = NULL) {
60		## Coordinates
61	!	bound <- graphics::par("usr")
62	!	ratio <- graphics::par("pin")[1] / graphics::par("pin")[2] # x/y ratio
63
64	!	to_unity <- function(x, y) {
65	!	list(x = (x - bound[1]) / (bound[2] - bound[1]) * ratio,
66	!	y = (y - bound[3]) / (bound[4] - bound[3]) / ratio)
67		}
68	!	to_usr <- function(x, y) {
69	!	list(x = bound[1] + x / ratio * (bound[2] - bound[1]),
70	!	y = bound[3] + y * ratio * (bound[4] - bound[3]))
71		}
72
73	!	xy <- to_unity(x = x, y = y)
74	!	x <- xy$x
75	!	y <- xy$y
76	!	n <- length(x)
77
78		## 8 positions: corners and side mid-points of the rectangle
79		## Position 7 (top right) is the most preferred
80	!	width <- graphics::strwidth(labels, units = "figure", cex = cex,
81	!	font = font, vfont = vfont)
82	!	height <- graphics::strheight(labels, units = "figure", cex = cex,
83	!	font = font, vfont = vfont)
84	!	width <- (width + 0.02) * ratio
85	!	height <- (height + 0.02) / ratio
86
87	!	makeoff <- function(pos) {
88	!	c(-1, -1, -1, 0, 0, 1, 1, 1)[pos] * (width / 2) +
89	!	1i * c(-1, 0, 1, -1, 1, -1, 0, 1)[pos] * (height / 2)
90		}
91
92		## Find intersection area of two rectangles
93	!	overlap <- function(xy1, off1, xy2, off2) {
94	!	w <- pmin(Re(xy1 + off1 / 2), Re(xy2 + off2 / 2)) -
95	!	pmax(Re(xy1 - off1 / 2), Re(xy2 - off2 / 2))
96	!	h <- pmin(Im(xy1 + off1 / 2), Im(xy2 + off2 / 2)) -
97	!	pmax(Im(xy1 - off1 / 2), Im(xy2 - off2 / 2))
98	!	w[w <= 0] <- 0
99	!	h[h <= 0] <- 0
100	!	w * h
101		}
102
103	!	objective <- function(gene) {
104	!	offset <- makeoff(gene)
105
106	!	if (!is.null(rectidx1)) {
107	!	area <- sum(overlap(xy[rectidx1] + offset[rectidx1], rectv[rectidx1],
108	!	xy[rectidx2] + offset[rectidx2], rectv[rectidx2]))
109		} else {
110	!	area <- 0
111		}
112
113		## Penalize labels which go outside the image area
114		## Count points outside of the image
115	!	a <- Re(xy + offset - rectv / 2) < 0 \| Re(xy + offset + rectv / 2) > ratio
116	!	b <- Im(xy + offset - rectv / 2) < 0 \| Im(xy + offset + rectv / 2) > 1 / ratio
117	!	outside <- sum(a \| b)
118	!	res <- 1000 * area + outside
119	!	res
120		}
121
122		# Make a list of label rectangles in their reference positions,
123		# centered over the map feature; the real labels are displaced
124		# from these positions so as not to overlap
125		# Note that some labels can be bigger than others
126	!	xy <- x + 1i * y
127	!	rectv <- width + 1i * height
128
129	!	rectidx1 <- rectidx2 <- array(0, (length(x)^2 - length(x)) / 2)
130	!	k <- 0
131	!	for (i in seq_along(x))
132	!	for (j in seq_len(i - 1)) {
133	!	k <- k + 1
134	!	rectidx1[k] <- i
135	!	rectidx2[k] <- j
136		}
137	!	maylap <- overlap(xy[rectidx1], 2 * rectv[rectidx1],
138	!	xy[rectidx2], 2 * rectv[rectidx2]) > 0
139	!	rectidx1 <- rectidx1[maylap]
140	!	rectidx2 <- rectidx2[maylap]
141
142		## Simulated annealing
143		## Initial state
144	!	gene <- rep(8, n)
145	!	score <- objective(gene)
146		## Initial "best" solution
147	!	bestgene <- gene
148	!	bestscore <- score
149	!	iter <- seq_len(iter)
150	!	temp <- 2.5
151	!	for (i in iter) {
152	!	k <- 1 # Energy evaluation count
153	!	for (j in iter) {
154	!	newgene <- gene
155	!	newgene[sample(n, 1)] <- sample(8, 1)
156	!	newscore <- objective(newgene)
157	!	if (newscore <= score \|\| stats::runif(1) < exp((score - newscore) / temp)) {
158		## keep the new set if it has the same or better score or
159		## if it's worse randomly based on the annealing criteria
160	!	k <- k + 1
161	!	score <- newscore
162	!	gene <- newgene
163		}
164	!	if (score <= bestscore) {
165	!	bestscore <- score
166	!	bestgene <- gene
167		}
168	!	if (bestscore == 0 \|\| k == 10) break
169		}
170	!	if (bestscore == 0) break
171	!	temp <- 0.9 * temp
172		}
173
174	!	nx <- Re(xy + makeoff(bestgene))
175	!	ny <- Im(xy + makeoff(bestgene))
176
177	!	xy <- to_usr(x = nx, y = ny)
178	!	xy$labels <- labels
179	!	xy
180		}
181
182		#' Shadow Text
183		#'
184		#' @param x,y A [`numeric`] vector. If `y` is `NULL`, an attempt is made to
185		#' interpret `x` in a suitable way (see [grDevices::xy.coords()]).
186		#' @param labels A [`character`] vector specifying the text to be written.
187		#' @param width Thickness of the shadow, as a fraction of the plotting size.
188		#' @param theta Angles for plotting the background.
189		#' @param cex A [`numeric`] character expansion factor.
190		#' @param col The color to be used for the text.
191		#' @param bg The color to be used for the shadow.
192		#' @param font,vfont The font to be used (see [graphics::text()]).
193		#' @param ... Further parameters to be passed to [graphics::text()].
194		#' @return
195		#' `text_shadow()` is called it for its side-effects: it results in a graphic
196		#' being displayed.
197		#' @author N. Frerebeau
198		#' @keywords internal
199		#' @noRd
200		text_shadow <- function(x, y = NULL, labels = seq_along(x$x),
201		width = 1/10, theta = seq(0, 2 * pi, length.out = 50),
202		cex = graphics::par("cex"), col = graphics::par("fg"),
203		bg = graphics::par("bg"), font = NULL, vfont = NULL, ...) {
204
205	!	x <- grDevices::xy.coords(x = x, y = y)
206
207	!	xo <- width * graphics::strwidth("M", units = "user", cex = cex, font = font, vfont = vfont)
208	!	yo <- width * graphics::strheight("X", units = "user", cex = cex, font = font, vfont = vfont)
209
210	!	for (i in theta) {
211	!	graphics::text(x = x$x + cos(i) * xo, y = x$y + sin(i) * yo, labels = labels,
212	!	col = bg, cex = cex, font = font, vfont = vfont, ...)
213		}
214
215	!	graphics::text(x = x$x, y = x$y, labels = labels, col = col, cex = cex,
216	!	font = font, vfont = vfont, ...)
217
218	!	invisible(NULL)
219		}

1		# MEAN
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_mean
7		#' @aliases ternary_mean,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_mean",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, ...) {
12	1x	x <- gmean(x)
13	1x	y <- gmean(y)
14	1x	z <- gmean(z)
15
16	1x	pt <- ternary_points(x = x, y = y, z = z, ...)
17	1x	invisible(pt)
18		}
19		)
20
21		#' @export
22		#' @rdname ternary_mean
23		#' @aliases ternary_mean,ANY,missing,missing-method
24		setMethod(
25		f = "ternary_mean",
26		signature = c(x = "ANY", y = "missing", z = "missing"),
27		definition = function(x, ...) {
28	1x	xyz <- grDevices::xyz.coords(x)
29	1x	pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z, ...)
30	1x	invisible(pt)
31		}
32		)
33
34		#' Geometric Mean
35		#'
36		#' @param x A [`numeric`] vector.
37		#' @param trim A length-one [`numeric`] vector specifying the fraction (0 to 0.5)
38		#' of observations to be trimmed from each end of `x` before the mean is
39		#' computed.
40		#' @param na.rm A [`logical`] scalar: should `NA` values be stripped before the
41		#' computation proceeds?
42		#' @return A [`numeric`] vector.
43		#' @keywords internal
44		#' @noRd
45		gmean <- function(x, trim = 0, na.rm = FALSE) {
46	504x	index <- is.finite(x) & x > 0
47	504x	exp(mean(log(unclass(x)[index]), trim = trim, na.rm = na.rm))
48		}

1		# TERNARY GRID
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_grid
7		ternary_grid <- function(primary = NULL, secondary = NULL,
8		center = NULL, scale = NULL,
9		col.primary = "darkgray", col.secondary = "lightgray",
10		lty.primary = "dashed", lty.secondary = "dotted",
11		lwd.primary = 1, lwd.secondary = lwd.primary) {
12		## Center and scale
13	21x	if (is.null(center)) center <- get_center()
14	21x	if (is.null(scale)) scale <- get_scale()
15
16		## Primary grid
17	21x	if (is.null(primary) \|\| (!anyNA(primary) && length(primary) == 1 && primary >= 1)) {
18	5x	if (is.null(primary)) primary <- graphics::par("xaxp")[3L]
19	21x	i <- seq(from = 0, to = 1, length.out = primary + 1)
20	21x	.ternary_grid(i, center = center, scale = scale,
21	21x	col = col.primary, lty = lty.primary, lwd = lwd.primary)
22		}
23
24		## Secondary grid
25	21x	if (!is.null(secondary) && !is.na(secondary) && length(secondary) == 1 && secondary > primary) {
26	10x	i <- seq(from = 0, to = 1, length.out = secondary + 1)
27	10x	.ternary_grid(i, center = center, scale = scale,
28	10x	col = col.secondary, lty = lty.secondary, lwd = lwd.secondary)
29		}
30
31	21x	invisible(NULL)
32		}
33
34		.ternary_grid <- function(x, center = NULL, scale = NULL,
35		col = "lightgray", lty = "dotted", lwd = 1, n = 100) {
36		## Reset values if needed
37	29x	if (!is.null(center) && all(center == 1)) center <- NULL
38	29x	if (!is.null(scale) && scale == 1) scale <- NULL
39
40	31x	x <- x[!(x == 0 \| x == 1)]
41	31x	if (is.null(scale)) {
42	29x	for (i in x) {
43	167x	start <- matrix(data = c(i, 0, 1 - i, 1 - i, i, 0, 0, 1 - i, i), ncol = 3)
44	167x	end <- matrix(data = c(i, 1 - i, 0, 0, i, 1 - i, 1 - i, 0, i), ncol = 3)
45
46	167x	start <- coordinates_ternary(start, center = center)
47	167x	end <- coordinates_ternary(end, center = center)
48
49	167x	graphics::segments(
50	167x	x0 = start$x, x1 = end$x,
51	167x	y0 = start$y, y1 = end$y,
52	167x	lty = lty, lwd = lwd, col = col
53		)
54		}
55		} else {
56	2x	for (i in x) {
57	8x	start <- matrix(data = c(i, 0, 1 - i, 1 - i, i, 0, 0, 1 - i, i), ncol = 3)
58	8x	end <- matrix(data = c(i, 1 - i, 0, 0, i, 1 - i, 1 - i, 0, i), ncol = 3)
59	8x	start <- list(x = start[, 2] + start[, 3] / 2, y = start[, 3] * sqrt(3) / 2)
60	8x	end <- list(x = end[, 2] + end[, 3] / 2, y = end[, 3] * sqrt(3) / 2)
61
62	8x	mapply(
63	8x	FUN = function(x_from, x_to, y_from, y_to, n, center, scale) {
64	24x	x <- seq(x_from, x_to, length.out = n)
65	24x	y <- seq(y_from, y_to, length.out = n)
66	24x	z <- coordinates_cartesian(x, y)
67	24x	zz <- coordinates_ternary(z, center = center, scale = scale)
68	24x	graphics::lines(
69	24x	zz,
70	24x	lty = lty, lwd = lwd, col = col
71		)
72		},
73	8x	x_from = start$x, x_to = end$x,
74	8x	y_from = start$y, y_to = end$y,
75	8x	MoreArgs = list(n = 100, center = center, scale = scale)
76		)
77		}
78		}
79		}

1		# TERNARY AXES
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_axis
7		ternary_axis <- function(side, at = NULL, labels = TRUE, tick = TRUE,
8		center = NULL, scale = NULL,
9		font = NA, lty = "solid",
10		lwd = 1, lwd.ticks = lwd,
11		col = NULL, col.ticks = NULL, ...) {
12		## Center and scale
13	3x	if (is.null(center)) center <- get_center()
14	3x	if (is.null(scale)) scale <- get_scale()
15
16		## Graphical parameters
17	138x	if (is.na(font)) font <- list(...)$font.axis %\|\|% graphics::par("font.axis")
18	!	if (is.null(col)) col <- list(...)$col.axis %\|\|% graphics::par("col.axis")
19	138x	if (is.null(col.ticks)) col.ticks <- col
20	138x	cex <- list(...)$cex.axis %\|\|% graphics::par("cex.axis")
21	138x	tcl <- list(...)$tcl %\|\|% graphics::par("tcl")
22
23		## Ticks and labels position
24	138x	if (is.null(at)) {
25	138x	at <- seq(from = 0, to = 1, length.out = graphics::par("xaxp")[3L] + 1)
26	138x	at <- at[!(at == 0 \| at == 1)]
27		}
28
29	138x	axis_degree <- c(120, 240, 0)[side]
30	138x	axis_radian <- c(0, 2 * pi / 3, 4 * pi / 3)[side]
31
32	138x	pos <- matrix(data = 0, nrow = length(at), ncol = 3)
33	138x	pos[, side] <- at
34	138x	pos[, c(2, 3, 1)[side]] <- 1 - at
35	138x	pos <- coordinates_ternary(pos, center = center, scale = scale)
36
37	138x	h <- abs(tcl * graphics::strheight("1", cex = 1))
38	138x	dx <- sin(pi / 6) * h
39	138x	dy <- cos(pi / 6) * h
40	138x	tick_start <- matrix(c(pos$x, pos$y), ncol = 2)
41	138x	tick_end <- matrix(c(pos$x + dx * c(1, 1, -2)[side],
42	138x	pos$y + dy * c(-1, 1, 0)[side]), ncol = 2)
43
44		## Labels
45	138x	if (!isFALSE(labels)) {
46	138x	if (length(labels) != length(at)) labels <- round(at * 100)
47	138x	graphics::text(x = tick_end, label = labels, srt = axis_degree,
48	138x	cex = cex, col = col, font = font, adj = c(1, 0.5))
49		}
50
51		## Ticks
52	138x	if (tick) {
53	138x	axis_line <- rotate(matrix(c(0, 0, 1, 0), ncol = 2), theta = axis_radian)
54	138x	graphics::segments(
55	138x	x0 = axis_line[1, 1], x1 = axis_line[2, 1],
56	138x	y0 = axis_line[1, 2], y1 = axis_line[2, 2],
57	138x	col = col, lwd = lwd, lty = lty
58		)
59	138x	graphics::segments(
60	138x	x0 = tick_start[, 1], x1 = tick_end[, 1],
61	138x	y0 = tick_start[, 2], y1 = tick_end[, 2],
62	138x	col = col.ticks, lwd = lwd.ticks, lty = lty
63		)
64		}
65
66	138x	invisible(NULL)
67		}

1		# TERNARY SOIL DIAGRAM
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname triangle_soil
7		triangle_soil_hypres <- function(labels = TRUE, symbol = FALSE, ...) {
8	1x	.triangle_soil("HYPRES", labels = labels, symbol = symbol, ...)
9	1x	invisible(NULL)
10		}
11
12		#' @export
13		#' @rdname triangle_soil
14		triangle_soil_folk <- function(labels = TRUE, symbol = FALSE, ...) {
15	1x	.triangle_soil("folk1954", labels = labels, symbol = symbol, ...)
16	1x	invisible(NULL)
17		}
18
19		#' @export
20		#' @rdname triangle_soil
21		triangle_soil_shepard <- function(labels = TRUE, symbol = FALSE, ...) {
22	1x	.triangle_soil("shepard1954", labels = labels, symbol = symbol, ...)
23	1x	invisible(NULL)
24		}
25
26		#' @export
27		#' @rdname triangle_soil
28		triangle_soil_usda <- function(labels = TRUE, symbol = FALSE, ...) {
29	1x	.triangle_soil("USDA1951", labels = labels, symbol = symbol, ...)
30	1x	invisible(NULL)
31		}
32
33		.triangle_soil <- function(chart, labels = TRUE, symbol = FALSE, ...) {
34		## Graphical parameters
35	4x	cex.lab <- list(...)$cex.lab %\|\|% graphics::par("cex.lab")
36	4x	col.lab <- list(...)$col.lab %\|\|% graphics::par("col.lab")
37	4x	font.lab <- list(...)$font.lab %\|\|% graphics::par("font.lab")
38
39	4x	poly <- .soil[[chart]]
40	4x	txt_lab <- unique(poly$label)
41	4x	txt_symb <- unique(poly$symbol)
42
43	4x	poly$label <- factor(poly$label, levels = unique(poly$label))
44	4x	poly <- split(poly, f = poly$label)
45
46	4x	for (i in poly) {
47	37x	ternary_polygon(i, ...)
48		}
49	4x	if (labels) {
50	4x	lab <- lapply(X = poly, FUN = function(x) colMeans(x[, c(1, 2, 3, 4)]))
51	4x	txt <- if (symbol && !all(txt_symb == "")) txt_symb else txt_lab
52
53	4x	for (i in seq_along(lab)) {
54	37x	ternary_text(
55	37x	x = lab[[i]][[1]],
56	37x	y = lab[[i]][[2]],
57	37x	z = lab[[i]][[3]],
58	37x	labels = txt[[i]],
59	37x	srt = lab[[i]][[4]],
60	37x	cex = cex.lab,
61	37x	col = col.lab,
62	37x	font = font.lab
63		)
64		}
65		}
66		}

1		# ELLIPSE
2		#' @include AllGenerics.R
3		NULL
4
5		# Ellipse ======================================================================
6		#' @export
7		#' @rdname ternary_ellipse
8		#' @aliases ternary_ellipse,numeric,numeric,numeric-method
9		setMethod(
10		f = "ternary_ellipse",
11		signature = c(x = "numeric", y = "numeric", z = "numeric"),
12		definition = function(x, y, z, radius = 1, ...) {
13		## ILR
14	3x	coda <- cbind(x, y, z)
15	3x	ratio <- ilr(coda)
16
17		## Compute ellipse
18	3x	mu <- colMeans(ratio)
19	3x	sigma <- stats::cov(ratio)
20		# rob <- robustbase::covMcd(ratio)
21		# mu <- rob$center
22		# sigma <- rob$cov
23	3x	xy <- ellipse(sigma = sigma, mu = mu, radius = radius)
24
25	3x	for (i in seq_along(xy)) {
26	5x	tern <- ilr_inv(xy[[i]]) # Inverse transform
27	5x	coords <- coordinates_ternary(tern)
28	5x	graphics::polygon(x = coords$x, y = coords$y, ...)
29		}
30
31	3x	invisible(list(x = x, y = y, z = z))
32		}
33		)
34
35		#' @export
36		#' @rdname ternary_ellipse
37		#' @aliases ternary_ellipse,ANY,missing,missing-method
38		setMethod(
39		f = "ternary_ellipse",
40		signature = c(x = "ANY", y = "missing", z = "missing"),
41		definition = function(x, radius = 1, ...) {
42	1x	xyz <- grDevices::xyz.coords(x)
43	1x	pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
44	1x	radius = radius, ...)
45	1x	invisible(pt)
46		}
47		)
48
49		#' Computes an Ellipse
50		#'
51		#' @param sigma A square positive definite \eqn{2 \times 2}{2 x 2} covariance
52		#' or correlation `matrix`.
53		#' @param mu A length-two [`numeric`] vector giving the centre of the ellipse.
54		#' @param scale If `sigma` is a correlation matrix, then the standard deviations
55		#' of each parameter can be given in the scale parameter.
56		#' Defaults to `c(1, 1)`, so no rescaling will be done.
57		#' @param level A length-\eqn{k} [`numeric`] vector giving the confidence level
58		#' of a pairwise confidence region.
59		#' @param radius The size of the ellipse may also be controlled by specifying
60		#' the value of a t-statistic on its boundary.
61		#' @param n A length-one [`numeric`] vector specifying the number of points used
62		#' in the ellipse.
63		#' @note Adapted from [ellipse::ellipse()].
64		#' @return
65		#' A [`list`] of \eqn{k} \eqn{n \times 2}{n x 2} `matrix`, suitable for
66		#' plotting.
67		#' @keywords internal
68		#' @noRd
69		ellipse <- function(sigma, mu = c(0, 0), scale = c(1, 1), level = 0.95,
70		radius = sqrt(stats::qchisq(level, 2)), n = 100, ...) {
71	3x	r <- sigma[1, 2]
72
73	3x	if (missing(scale)) {
74	3x	scale <- sqrt(diag(sigma))
75	3x	if (scale[1] > 0) r <- r / scale[1]
76	3x	if (scale[2] > 0) r <- r / scale[2]
77		}
78
79	3x	r <- min(max(r, -1), 1) # clamp to -1..1, in case of rounding errors
80	3x	d <- acos(r)
81	3x	a <- seq(0, 2 * pi, len = n)
82
83	3x	lapply(
84	3x	X = radius,
85	3x	FUN = function(x) {
86	5x	matrix(
87	5x	data = c(x * scale[1] * cos(a + d / 2) + mu[1],
88	5x	x * scale[2] * cos(a - d / 2) + mu[2]),
89	5x	nrow = n,
90	5x	ncol = 2,
91	5x	dimnames = list(NULL, c("x", "y"))
92		)
93		}
94		)
95		}
96
97		# Confidence ellipse ===========================================================
98		#' @export
99		#' @rdname ternary_ellipse
100		#' @aliases ternary_confidence,numeric,numeric,numeric-method
101		setMethod(
102		f = "ternary_confidence",
103		signature = c(x = "numeric", y = "numeric", z = "numeric"),
104		definition = function(x, y, z, level = 0.95, ...) {
105	1x	df1 <- 2
106	1x	df2 <- length(x) - 2
107	1x	radius <- sqrt(stats::qf(p = level, df1, df2) * df1 / df2)
108	1x	ternary_ellipse(x, y, z, radius = radius, ...)
109		}
110		)
111
112		#' @export
113		#' @rdname ternary_ellipse
114		#' @aliases ternary_confidence,ANY,missing,missing-method
115		setMethod(
116		f = "ternary_confidence",
117		signature = c(x = "ANY", y = "missing", z = "missing"),
118		definition = function(x, level = 0.95, ...) {
119	1x	x <- grDevices::xyz.coords(x)
120	1x	methods::callGeneric(x = x$x, y = x$y, z = x$z, level = level, ...)
121		}
122		)
123
124		# Probability ellipse ==========================================================
125		#' @export
126		#' @rdname ternary_ellipse
127		#' @aliases ternary_tolerance,numeric,numeric,numeric-method
128		setMethod(
129		f = "ternary_tolerance",
130		signature = c(x = "numeric", y = "numeric", z = "numeric"),
131		definition = function(x, y, z, level = 0.95, ...) {
132	1x	radius <- sqrt(stats::qchisq(p = level, df = 2))
133	1x	ternary_ellipse(x, y, z, radius = radius, ...)
134		}
135		)
136
137		#' @export
138		#' @rdname ternary_ellipse
139		#' @aliases ternary_tolerance,ANY,missing,missing-method
140		setMethod(
141		f = "ternary_tolerance",
142		signature = c(x = "ANY", y = "missing", z = "missing"),
143		definition = function(x, level = 0.95, ...) {
144	1x	x <- grDevices::xyz.coords(x)
145	1x	methods::callGeneric(x = x$x, y = x$y, z = x$z, level = level, ...)
146		}
147		)

1		# IMAGE
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_image
7		#' @aliases ternary_image,function-method
8		setMethod(
9		f = "ternary_image",
10		signature = c(f = "function"),
11		definition = function(f, n = 48, palette = NULL, ...) {
12
13	4x	tri <- .triangle_center(n)
14	4x	xyz <- coordinates_cartesian(tri$x, tri$y)
15	4x	val <- f(xyz$x, xyz$y, xyz$z, ...)
16	4x	col <- map_color(val, palette = palette)
17
18	4x	coords <- .triangle_vertex(tri$x, tri$y, tri$direction, tri$resolution)
19	4x	for (i in seq_along(coords)) {
20	2128x	polygon(coords[[i]][, 1], coords[[i]][, 2], col = col[i], border = NA)
21		}
22
23	4x	invisible(list(levels = val, colors = col))
24		}
25		)
26
27		#' @export
28		#' @rdname ternary_tile
29		#' @aliases tile_bin,numeric,numeric,numeric-method
30		setMethod(
31		f = "tile_bin",
32		signature = c(x = "numeric", y = "numeric", z = "numeric"),
33		definition = function(x, y, z) {
34
35	1x	total <- x + y + z
36	1x	a <- x / total
37	1x	b <- y / total
38	1x	c <- z / total
39
40	1x	function(x, y, z) {
41	1x	tri <- .triangle_center(sqrt(length(x)))
42	1x	coords <- .triangle_vertex(tri$x, tri$y, tri$direction, tri$resolution)
43
44	1x	count <- numeric(length(coords))
45	1x	for (i in seq_along(coords)) {
46	576x	xyz <- coordinates_cartesian(coords[[i]][, 1], coords[[i]][, 2])
47	576x	count[[i]] <- sum(min(xyz$x) <= a & a < max(xyz$x) &
48	576x	min(xyz$y) <= b & b < max(xyz$y) &
49	576x	min(xyz$z) <= c & c < max(xyz$z))
50		}
51
52	1x	count[count == 0] <- NA
53	1x	count
54		}
55		}
56		)
57
58		#' @export
59		#' @rdname ternary_tile
60		#' @aliases tile_bin,ANY,missing,missing-method
61		setMethod(
62		f = "tile_bin",
63		signature = c(x = "ANY", y = "missing", z = "missing"),
64		definition = function(x) {
65	1x	xyz <- grDevices::xyz.coords(x)
66	1x	methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z)
67		}
68		)
69
70		#' @export
71		#' @rdname ternary_tile
72		#' @aliases tile_density,numeric,numeric,numeric-method
73		setMethod(
74		f = "tile_density",
75		signature = c(x = "numeric", y = "numeric", z = "numeric"),
76		definition = function(x, y, z) {
77		## ILR
78	1x	coda <- cbind(x, y, z)
79	1x	ratio <- ilr(coda)
80
81		## Compute KDE
82	1x	function(x, y, z) {
83	1x	xyz <- cbind(x, y, z)
84	1x	xy <- ilr(xyz)
85	1x	dens <- kde(
86	1x	x = ratio[, 1],
87	1x	y = ratio[, 2],
88	1x	gx = sort(unique(xy[, 1])),
89	1x	gy = sort(unique(xy[, 2]))
90		)
91
92	1x	i <- as.numeric(as.factor(rank(xy[, 1])))
93	1x	j <- as.numeric(as.factor(rank(xy[, 2])))
94
95	1x	dens$z[cbind(i, j)]
96		}
97		}
98		)
99
100		#' @export
101		#' @rdname ternary_tile
102		#' @aliases tile_density,ANY,missing,missing-method
103		setMethod(
104		f = "tile_density",
105		signature = c(x = "ANY", y = "missing", z = "missing"),
106		definition = function(x) {
107	1x	xyz <- grDevices::xyz.coords(x)
108	1x	methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z)
109		}
110		)
111
112		#' @export
113		#' @rdname ternary_tile
114		#' @aliases tile_interpolate,numeric,numeric,numeric-method
115		setMethod(
116		f = "tile_interpolate",
117		signature = c(x = "numeric", y = "numeric", z = "numeric"),
118		definition = function(x, y, z, value, method = "linear", ...) {
119		## Validation
120	1x	assert_package("interp")
121	1x	assert_length(value, length(x))
122
123		## ILR
124	1x	coda <- cbind(x, y, z)
125	1x	ratio <- ilr(coda)
126
127		## Interpolate
128	1x	function(x, y, z) {
129	1x	xyz <- cbind(x, y, z)
130	1x	xy <- ilr(xyz)
131
132	1x	interp <- interp::interp(
133	1x	x = ratio[, 1],
134	1x	y = ratio[, 2],
135	1x	z = value,
136	1x	xo = sort(unique(xy[, 1])),
137	1x	yo = sort(unique(xy[, 2])),
138	1x	method = method,
139		...
140		)
141
142	1x	i <- as.numeric(as.factor(rank(xy[, 1])))
143	1x	j <- as.numeric(as.factor(rank(xy[, 2])))
144
145	1x	interp$z[cbind(i, j)]
146		}
147		}
148		)
149
150		#' @export
151		#' @rdname ternary_tile
152		#' @aliases tile_interpolate,ANY,missing,missing-method
153		setMethod(
154		f = "tile_interpolate",
155		signature = c(x = "ANY", y = "missing", z = "missing"),
156		definition = function(x, value, method = "linear", ...) {
157	1x	xyz <- grDevices::xyz.coords(x)
158	1x	methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
159	1x	value = value, method = method, ...)
160		}
161		)
162
163		#' Tile Center Coordinates
164		#'
165		#' Computes tile center cartesian coordinates.
166		#' @param resolution A length-one [`integer`] vector specifying the maximum
167		#' number of tiles on each axis.
168		#' @return
169		#' A [`list`] with the following elements:
170		#' \describe{
171		#' \item{`x`}{x cartesian coordinates.}
172		#' \item{`y`}{y cartesian coordinates.}
173		#' \item{`direction`}{`1` means up, `-1` means down.}
174		#' \item{`resolution`}{}
175		#' }
176		#' @examples
177		#' .triangle_center(5)
178		#' @keywords internal
179		#' @noRd
180		.triangle_center <- function(resolution) {
181
182	5x	offset <- 1 / resolution / 2L
183	5x	height <- the$top / resolution
184
185	5x	X <- seq(from = offset, to = 1 - offset, by = offset)
186	5x	X <- lapply(
187	5x	X = seq_len(resolution),
188	5x	FUN = function(step) {
189	116x	up <- X[seq(from = step, to = (2L * resolution) - step, by = 2L)]
190	116x	down <- integer(0)
191	116x	if (step != resolution) {
192	111x	down <- X[seq(from = step + 1, to = (2L * resolution) - step - 1, by = 2L)]
193		}
194	116x	c(rbind(up, c(down, NA)))[-(resolution - step + 1) * 2L]
195		}
196		)
197
198	5x	in_row <- 2L * rev(seq_len(resolution)) - 1L
199	5x	is_down <- (1 + unlist(lapply(in_row, seq_len))) %% 2L
200
201	5x	Y <- seq(from = height / 3, to = the$top - (2 * height / 3), length.out = resolution)
202	5x	Y <- rep(Y[seq_len(resolution)], in_row)
203	5x	Y <- Y + (is_down * height / 3)
204
205	5x	dir <- rep(1, length(is_down))
206	5x	dir[is_down == 1] <- -1
207
208	5x	list(
209	5x	x = unlist(X),
210	5x	y = Y,
211	5x	direction = dir,
212	5x	resolution = resolution
213		)
214		}
215
216		#' Tile Vertex Coordinates
217		#'
218		#' Computes tile vertex cartesian coordinates.
219		#' @param x,y A [`numeric`] vector giving the cartesian coordinates of the
220		#' center.
221		#' @param direction An [`integer`] vector specifying the triangle direction
222		#' (`1` means up, `-1` means down).
223		#' @param resolution A length-one [`integer`] vector specifying the maximum
224		#' number of tiles on each axis.
225		#' @return
226		#' A [`list`] of `numeric` [`matrix`].
227		#' @examples
228		#' m <- .triangle_center(5)
229		#' .triangle_vertex(m$x, m$y, m$direction, m$resolution)
230		#' @keywords internal
231		#' @noRd
232		.triangle_vertex <- function(x, y, direction, resolution) {
233	5x	n <- length(x)
234
235	5x	width <- 1 / resolution / 2
236	5x	height <- the$top / resolution / 3
237
238	5x	tiles <- vector(mode = "list", length = n)
239	5x	for (i in seq_len(n)) {
240	2704x	tiles[[i]] <- matrix(
241	2704x	data = c(x[i] + c(0, width, -width), y[i] + c(2 * height, -height, -height) * direction[i]),
242	2704x	ncol = 2, dimnames = list(NULL, c("x", "y"))
243		)
244		}
245
246	5x	tiles
247		}

1		# TERNARY WINDOW
2		#' @include AllGenerics.R
3		NULL
4
5		#' Set up Ternary Coordinates for Graphics Window
6		#'
7		#' @param xlim A length-three [`numeric`] vector giving the `x` limits in the
8		#' range \eqn{[0,1]}.
9		#' @param ylim A length-three [`numeric`] vector giving the `y` limits in the
10		#' range \eqn{[0,1]}.
11		#' @param zlim A length-three [`numeric`] vector giving the `z` limits in the
12		#' range \eqn{[0,1]}.
13		#' @param xlab,ylab,zlab A [`character`] string giving a label for the x, y and
14		#' z axes.
15		#' @return
16		#' `ternary_window()` is called it for its side-effects
17		#' (see [graphics::plot.window()]).
18		#' @keywords internal
19		#' @noRd
20		ternary_window <- function(xlim = NULL, ylim = NULL, zlim = NULL,
21		xlab = NULL, ylab = NULL, zlab = NULL,
22		cex = 1) {
23
24	91x	n_null <- is.null(xlim) + is.null(ylim) + is.null(zlim)
25
26	91x	if (n_null == 3) {
27	87x	dx <- max(graphics::strwidth(c(xlab, ylab, zlab), cex = cex)) / 2
28	87x	rx <- expand_range(c(0, 1), add = dx)
29	87x	lim <- list(
30	87x	x = rx,
31	87x	y = the$top / 2 + c(-1, 1) * diff(rx) / 2
32		)
33		}
34	91x	if (n_null == 2) {
35	1x	stop(tr_("You must provide at least two coordinates ranges."), call. = FALSE)
36		}
37	90x	if (n_null < 2) {
38	3x	xlims <- if (is.null(xlim)) 1 - ylim - zlim else xlim
39	3x	ylims <- if (is.null(ylim)) 1 - xlim - zlim else ylim
40	3x	zlims <- if (is.null(zlim)) 1 - xlim - ylim else zlim
41
42	3x	assert_length(xlims, 3)
43	3x	assert_length(ylims, 3)
44	3x	assert_length(zlims, 3)
45
46	3x	lim <- coordinates_ternary(abs(xlims), abs(ylims), abs(zlims))
47	3x	lim$x <- range(lim$x)
48	3x	lim$y <- range(lim$y)
49		}
50	90x	graphics::plot.window(xlim = lim$x, ylim = lim$y, asp = 1)
51		}

1		# TERNARY PHASE DIAGRAM
2		#' @include AllGenerics.R
3		NULL
4
5		# CAS phase diagram ============================================================
6		#' @export
7		#' @rdname triangle_phase_cas
8		triangle_phase_cas <- function(labels = TRUE, symbol = FALSE,
9		mol = FALSE, ...) {
10	2x	oxide_mass <- c(CaO = 56.0774, Al2O3 = 101.9600, SiO2 = 60.0800)
11	2x	.triangle_phases("CAS", oxide_mass = oxide_mass, mol = mol,
12	2x	labels = labels, symbol = symbol, ...)
13
14	2x	invisible(NULL)
15		}
16
17		#' @export
18		#' @rdname triangle_phase_cas
19		triangle_phase_ceramic <- function(labels = TRUE, symbol = FALSE,
20		mol = FALSE, ...) {
21	2x	oxide_mass <- c(CaO = 56.0774, Al2O3 = 101.9600, SiO2 = 60.0800)
22	2x	.triangle_phases("ceramic", oxide_mass = oxide_mass, mol = mol,
23	2x	labels = labels, symbol = symbol, ...)
24
25	2x	invisible(NULL)
26		}
27
28		.triangle_phases <- function(chart, oxide_mass, mol = FALSE,
29		labels = TRUE, symbol = FALSE, ...) {
30		## Graphical parameters
31	4x	cex.lab <- list(...)$cex.lab %\|\|% graphics::par("cex.lab")
32	4x	col.lab <- list(...)$col.lab %\|\|% graphics::par("col.lab")
33	4x	font.lab <- list(...)$font.lab %\|\|% graphics::par("font.lab")
34
35	4x	poly <- .phases[[chart]]
36
37	4x	if (!mol) {
38	2x	oxyde_mol <- as.matrix(poly[, c(1, 2, 3)])
39
40		## Molar mass (g/mol)
41	2x	phase_mass <- oxyde_mol %*% oxide_mass
42
43		## Oxide weight (%)
44	2x	poly[, c(1, 2, 3)] <- t(t(oxyde_mol) * oxide_mass) / as.vector(phase_mass)
45		}
46
47	4x	lab <- poly[!duplicated(poly$label), ]
48	4x	txt <- if (symbol && !all(lab$symbol == "")) lab$symbol else lab$label
49
50	4x	poly$group <- factor(poly$group, levels = unique(poly$group))
51	4x	poly <- split(poly, f = poly$group)
52
53	4x	for (i in poly) {
54	40x	ternary_polygon(i, ...)
55		}
56	4x	if (labels) {
57	4x	ternary_points(lab, cex = cex.lab, col = col.lab, ...)
58	4x	ternary_text(lab, label = txt, pos = lab$pos,
59	4x	cex = cex.lab, col = col.lab, font = font.lab)
60		}
61		}

1		# CONVEX HULL
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_hull
7		#' @aliases ternary_hull,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_hull",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, center = FALSE, scale = FALSE, ...) {
12	1x	coords <- coordinates_ternary(x, y, z, center = center, scale = scale)
13	1x	hull <- grDevices::chull(coords)
14	1x	graphics::polygon(x = coords$x[hull], y = coords$y[hull], ...)
15
16	1x	coords <- utils::modifyList(coords, list(x = x, y = y, z = z))
17	1x	invisible(coords)
18		}
19		)
20
21		#' @export
22		#' @rdname ternary_hull
23		#' @aliases ternary_hull,ANY,missing,missing-method
24		setMethod(
25		f = "ternary_hull",
26		signature = c(x = "ANY", y = "missing", z = "missing"),
27		definition = function(x, center = FALSE, scale = FALSE, ...) {
28	1x	xyz <- grDevices::xyz.coords(x)
29	1x	coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
30	1x	center = center, scale = scale, ...)
31	1x	invisible(coords)
32		}
33		)

1		# PCA
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_pca
7		#' @aliases ternary_pca,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_pca",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, axis = 1, ...) {
12		## CLR
13	1x	coda <- cbind(x, y, z)
14	1x	ratio <- clr(coda)
15
16	1x	z <- ratio - rowMeans(ratio) # Center
17	1x	m <- colMeans(z)
18
19		## Get eigenvectors
20	1x	eig <- eigen(stats::cov(z))$vectors[, axis[[1L]]] + m
21
22		## Standard coordinates
23	1x	std <- z %*% eig
24
25	1x	lam <- seq(-5, 5, length.out = nrow(ratio))
26	1x	axe <- cbind(eig[1L] * lam, eig[2L] * lam, eig[3L] * lam) +
27	1x	cbind(m[1L] * (1 - lam), m[2L] * (1 - lam), m[3L] * (1 - lam))
28
29		## Inverse CLR
30	1x	axe <- clr_inv(axe)
31	1x	coords <- coordinates_ternary(axe)
32
33		## Plot
34	1x	graphics::lines(coords, ...)
35
36	1x	invisible(list(x = x, y = y, z = z))
37		}
38		)
39
40		#' @export
41		#' @rdname ternary_pca
42		#' @aliases ternary_pca,ANY,missing,missing-method
43		setMethod(
44		f = "ternary_pca",
45		signature = c(x = "ANY", y = "missing", z = "missing"),
46		definition = function(x, axis = 1, ...) {
47	1x	xyz <- grDevices::xyz.coords(x)
48	1x	coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z, axis = axis, ...)
49	1x	invisible(coords)
50		}
51		)

1		# TERNARY POINTS
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_points
7		#' @aliases ternary_points,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_points",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, center = FALSE, scale = FALSE, type = "p", ...) {
12	12x	pt <- coordinates_ternary(x, y, z, center = center, scale = scale)
13	12x	graphics::points(x = pt, type = type, ...)
14
15	12x	pt <- utils::modifyList(pt, list(x = x, y = y, z = z))
16	12x	invisible(pt)
17		}
18		)
19
20		#' @export
21		#' @rdname ternary_points
22		#' @aliases ternary_points,ANY,missing,missing-method
23		setMethod(
24		f = "ternary_points",
25		signature = c(x = "ANY", y = "missing", z = "missing"),
26		definition = function(x, center = FALSE, scale = FALSE, type = "p", ...) {
27	11x	xyz <- grDevices::xyz.coords(x)
28	11x	pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
29	11x	center = center, scale = scale,
30	11x	type = type, ...)
31	11x	invisible(pt)
32		}
33		)

1		# TERNARY CROSS-HAIRS
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_crosshairs
7		#' @aliases ternary_crosshairs,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_crosshairs",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z,
12		x_mark = TRUE, y_mark = TRUE, z_mark = TRUE, ...) {
13	3x	total <- x + y + z
14	3x	zero <- rep(0, length(x))
15
16	3x	if (x_mark) {
17	1x	ternary_segments(x0 = x, y0 = y, z0 = z,
18	1x	x1 = x / total, y1 = 1 - (x / total), z1 = zero, ...)
19		}
20	3x	if (y_mark) {
21	1x	ternary_segments(x0 = x, y0 = y, z0 = z,
22	1x	x1 = zero, y1 = y / total, z1 = 1 - (y / total), ...)
23		}
24	3x	if (z_mark) {
25	1x	ternary_segments(x0 = x, y0 = y, z0 = z,
26	1x	x1 = 1 - (z / total), y1 = zero, z1 = z / total, ...)
27		}
28
29	3x	invisible(list(x = x, y = y, z = z))
30		}
31		)
32
33		#' @export
34		#' @rdname ternary_crosshairs
35		#' @aliases ternary_crosshairs,ANY,missing,missing-method
36		setMethod(
37		f = "ternary_crosshairs",
38		signature = c(x = "ANY", y = "missing", z = "missing"),
39		definition = function(x, x_mark = TRUE, y_mark = TRUE, z_mark = TRUE, ...) {
40	3x	xyz <- grDevices::xyz.coords(x)
41	3x	pt <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
42	3x	x_mark = x_mark, y_mark = y_mark, z_mark = z_mark, ...)
43	3x	invisible(pt)
44		}
45		)

1		# TERNARY ARROWS
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_arrows
7		#' @aliases ternary_arrows,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_arrows",
10		signature = c(x0 = "numeric", y0 = "numeric", z0 = "numeric"),
11		definition = function(x0, y0, z0, x1 = x0, y1 = y0, z1 = z0, ...) {
12	1x	coords0 <- coordinates_ternary(x0, y0, z0)
13	1x	coords1 <- coordinates_ternary(x1, y1, z1)
14	1x	graphics::arrows(x0 = coords0$x, y0 = coords0$y,
15	1x	x1 = coords1$x, y1 = coords1$y, ...)
16
17	1x	invisible(NULL)
18		}
19		)

1		# TERNARY PLOT
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_plot
7		#' @aliases ternary_plot,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_plot",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, center = FALSE, scale = FALSE,
12		xlim = NULL, ylim = NULL, zlim = NULL,
13		xlab = NULL, ylab = NULL, zlab = NULL,
14		main = NULL, sub = NULL, ann = graphics::par("ann"),
15		axes = TRUE, frame.plot = axes,
16		panel.first = NULL, panel.last = NULL, ...) {
17
18		## Save and restore graphical parameters
19		## pty: square plotting region, independent of device size
20	91x	old_par <- graphics::par(pty = "s", no.readonly = TRUE)
21	91x	on.exit(graphics::par(old_par), add = TRUE)
22
23		## Graphical parameters
24	91x	fg <- list(...)$fg %\|\|% graphics::par("fg")
25	91x	cex.lab <- list(...)$cex.lab %\|\|% graphics::par("cex.lab")
26
27		## Open new window
28	91x	grDevices::dev.hold()
29	91x	on.exit(grDevices::dev.flush(), add = TRUE)
30	91x	graphics::plot.new()
31
32		## Set plotting coordinates
33	91x	ternary_window(xlim = xlim, ylim = ylim, zlim = zlim,
34	91x	xlab = xlab, ylab = ylab, zlab = zlab,
35	91x	cex = cex.lab)
36
37		## Reset center and scale
38	90x	assign("center", value = NULL, envir = the)
39	90x	assign("scale", value = NULL, envir = the)
40
41		## Compute ternary coordinates
42	90x	pt <- coordinates_ternary(x = x, y = y, z = z, center = center, scale = scale)
43
44		## Save center and scale for further use, e.g. grid or axes.
45	90x	assign("center", value = pt$center, envir = the)
46	90x	assign("scale", value = pt$scale, envir = the)
47
48		## Evaluate pre-plot expressions
49	90x	panel.first
50
51		## Plot
52	90x	graphics::points(x = pt, ...)
53
54		## Evaluate post-plot and pre-axis expressions
55	90x	panel.last
56
57		## Construct Axis
58	90x	if (axes) {
59	45x	ternary_axis(side = 1, center = pt$center, scale = pt$scale, col = fg)
60	45x	ternary_axis(side = 2, center = pt$center, scale = pt$scale, col = fg)
61	45x	ternary_axis(side = 3, center = pt$center, scale = pt$scale, col = fg)
62		}
63
64		## Plot frame
65	90x	if (frame.plot) {
66	88x	ternary_box(lty = "solid", lwd = 1, col = fg)
67		}
68
69		## Add annotation
70	90x	if (ann) {
71	86x	ternary_title(main = main, sub = sub, xlab = xlab, ylab = ylab,
72	86x	zlab = zlab, ...)
73		}
74
75	90x	invisible(pt)
76		}
77		)
78
79		#' @export
80		#' @rdname ternary_plot
81		#' @aliases ternary_plot,ANY,missing,missing-method
82		setMethod(
83		f = "ternary_plot",
84		signature = c(x = "ANY", y = "missing", z = "missing"),
85		definition = function(x, xlim = NULL, ylim = NULL, zlim = NULL,
86		xlab = NULL, ylab = NULL, zlab = NULL,
87		main = NULL, sub = NULL, ann = graphics::par("ann"),
88		axes = TRUE, frame.plot = axes,
89		panel.first = NULL, panel.last = NULL, ...) {
90
91	91x	xyz <- grDevices::xyz.coords(x, xlab = xlab, ylab = ylab, zlab = zlab)
92	91x	pt <- methods::callGeneric(
93	91x	x = xyz$x, y = xyz$y, z = xyz$z,
94	91x	xlim = xlim,
95	91x	ylim = ylim,
96	91x	zlim = zlim,
97	91x	xlab = xlab %\|\|% xyz$xlab %\|\|% "x",
98	91x	ylab = ylab %\|\|% xyz$ylab %\|\|% "y",
99	91x	zlab = zlab %\|\|% xyz$zlab %\|\|% "z",
100	91x	main = main, sub = sub, ann = ann,
101	91x	axes = axes,
102	91x	frame.plot = frame.plot,
103	91x	panel.first = panel.first,
104	91x	panel.last = panel.last,
105		...
106		)
107
108	90x	invisible(pt)
109		}
110		)

1		# TERNARY BOX
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_box
7		ternary_box <- function(lty = "solid", ...) {
8
9		## Graphical parameters
10	88x	col <- list(...)$col %\|\|% graphics::par("fg")
11	88x	lwd <- list(...)$lwd %\|\|% 1
12
13	88x	graphics::polygon(x = c(0, 0.5, 1), y = c(0, the$top, 0),
14	88x	border = col, lty = lty, lwd = lwd)
15
16	88x	invisible(NULL)
17		}

1		# TERNARY TEXT
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_text
7		#' @aliases ternary_text,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_text",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, center = FALSE, scale = FALSE,
12		labels = seq_along(x), ...) {
13	42x	coords <- coordinates_ternary(x, y, z, center = center, scale = scale)
14	42x	graphics::text(x = coords, labels = labels, ...)
15
16	42x	invisible(list(x = x, y = y, z = z))
17		}
18		)
19
20		#' @export
21		#' @rdname ternary_text
22		#' @aliases ternary_text,ANY,missing,missing-method
23		setMethod(
24		f = "ternary_text",
25		signature = c(x = "ANY", y = "missing", z = "missing"),
26		definition = function(x, center = FALSE, scale = FALSE, labels = seq_along(x$x), ...) {
27	5x	x <- grDevices::xyz.coords(x)
28	5x	force(labels)
29	5x	coords <- methods::callGeneric(x = x$x, y = x$y, z = x$z,
30	5x	center = center, scale = scale,
31	5x	labels = labels, ...)
32	5x	invisible(coords)
33		}
34		)

1		# TERNARY SEGMENTS
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_segments
7		#' @aliases ternary_segments,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_segments",
10		signature = c(x0 = "numeric", y0 = "numeric", z0 = "numeric"),
11		definition = function(x0, y0, z0, x1 = x0, y1 = y0, z1 = z0, ...) {
12	4x	coords0 <- coordinates_ternary(x0, y0, z0)
13	4x	coords1 <- coordinates_ternary(x1, y1, z1)
14	4x	graphics::segments(x0 = coords0$x, y0 = coords0$y,
15	4x	x1 = coords1$x, y1 = coords1$y, ...)
16
17	4x	invisible(NULL)
18		}
19		)

1		# TERNARY POLYGON
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_polygon
7		#' @aliases ternary_polygon,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_polygon",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, ...) {
12	78x	coords <- coordinates_ternary(x, y, z)
13	78x	graphics::polygon(x = coords, ...)
14
15	78x	invisible(list(x = x, y = y, z = z))
16		}
17		)
18
19		#' @export
20		#' @rdname ternary_polygon
21		#' @aliases ternary_polygon,ANY,missing,missing-method
22		setMethod(
23		f = "ternary_polygon",
24		signature = c(x = "ANY", y = "missing", z = "missing"),
25		definition = function(x, ...) {
26	78x	xyz <- grDevices::xyz.coords(x)
27	78x	coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z, ...)
28	78x	invisible(coords)
29		}
30		)

1		# TERNARY LINES
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ternary_lines
7		#' @aliases ternary_lines,numeric,numeric,numeric-method
8		setMethod(
9		f = "ternary_lines",
10		signature = c(x = "numeric", y = "numeric", z = "numeric"),
11		definition = function(x, y, z, type = "l", ...) {
12	143x	coords <- coordinates_ternary(x, y, z)
13	143x	graphics::lines(x = coords, type = type, ...)
14
15	143x	invisible(list(x = x, y = y, z = z))
16		}
17		)
18
19		#' @export
20		#' @rdname ternary_lines
21		#' @aliases ternary_lines,ANY,missing,missing-method
22		setMethod(
23		f = "ternary_lines",
24		signature = c(x = "ANY", y = "missing", z = "missing"),
25		definition = function(x, type = "l", ...) {
26	143x	xyz <- grDevices::xyz.coords(x)
27	143x	coords <- methods::callGeneric(x = xyz$x, y = xyz$y, z = xyz$z,
28	143x	type = type, ...)
29	143x	invisible(coords)
30		}
31		)

1		# INTERNAL ENVIRONMENT
2
3		the <- new.env(parent = emptyenv())
4		the$top <- sqrt(3) / 2
5		the$center <- NULL # Updated by plot()
6		the$scale <- NULL # Updated by plot()
7
8		get_center <- function(...) {
9	908x	get("center", envir = the)
10		}
11
12		get_scale <- function(...) {
13	908x	get("scale", envir = the)
14		}