1#' Estimates quantiles for each row (column) in a matrix
2#'
3#' Estimates quantiles for each row (column) in a matrix.
4#'
5#' @inheritParams rowAlls
6#'
7#' @param x An \code{\link[base]{integer}}, \code{\link[base]{numeric}} or
8#' \code{\link[base]{logical}} NxK \code{\link[base]{matrix}} with N >= 0.
9#'
10#' @param probs A \code{\link[base]{numeric}} \code{\link[base]{vector}} of J
11#' probabilities in [0, 1].
12#'
13#' @param type An \code{\link[base]{integer}} specify the type of estimator.
14#' See \code{\link[stats]{quantile}} for more details.
15#'
16#' @param ... Additional arguments passed to \code{\link[stats]{quantile}}.
17#'
18#' @param drop If TRUE, singleton dimensions in the result are dropped,
19#' otherwise not.
20#'
21#' @return Returns a NxJ (KxJ) \code{\link[base]{matrix}}, where N (K) is the
22#' number of rows (columns) for which the J quantiles are calculated.
23#' The return type is either integer or numeric depending on \code{type}.
24#'
25#' @example incl/rowQuantiles.R
26#'
27#' @author Henrik Bengtsson
28#' @seealso \code{\link[stats]{quantile}}.
29#' @keywords array iteration robust univar
30#'
31#' @importFrom stats quantile
32#' @export
33rowQuantiles <- function(x, rows = NULL, cols = NULL,
34                         probs = seq(from = 0, to = 1, by = 0.25),
35                         na.rm = FALSE, type = 7L, ..., useNames = NA, drop = TRUE) {
36  # Argument 'x':
37  if (!is.matrix(x)) defunctShouldBeMatrix(x)
38  if (!is.numeric(x) && !is.integer(x) && !is.logical(x)) {
39    .Defunct(msg = sprintf("Argument 'x' is of type %s. Only 'integer', 'numeric', and 'logical' is supported.", sQuote(storage.mode(x))))  #nolint
40  }
41
42  # Argument 'probs':
43  if (anyMissing(probs)) {
44    stop(sprintf("Argument '%s' must not contain missing values", "probs"))
45  }
46  eps <- 100 * .Machine$double.eps
47  if (any((probs < -eps | probs > 1 + eps))) {
48    stop(sprintf("Argument '%s' is out of range [0-eps, 1+eps]: %g", "probs", probs))
49  }
50
51  # Apply subset
52  if (!is.null(rows) && !is.null(cols)) x <- x[rows, cols, drop = FALSE]
53  else if (!is.null(rows)) x <- x[rows, , drop = FALSE]
54  else if (!is.null(cols)) x <- x[, cols, drop = FALSE]
55
56  # Coerce?
57#  if (is.logical(x)) {
58#    storage.mode(x) <- "integer"
59#  }
60
61  # Argument 'x':
62  nrow <- nrow(x)
63  ncol <- ncol(x)
64
65  # Allocate result
66  na_value <- NA_real_
67  if (type != 7L) storage.mode(na_value) <- storage.mode(x)
68  q <- matrix(na_value, nrow = nrow, ncol = length(probs))
69
70  if (nrow > 0L && ncol > 0L) {
71    na_rows <- rowAnyMissings(x)
72    has_na <- any(na_rows)
73    if (!has_na) na.rm <- FALSE
74
75    if (!has_na && type == 7L) {
76      n <- ncol
77      idxs <- 1 + (n - 1) * probs
78      idxs_lo <- floor(idxs)
79      idxs_hi <- ceiling(idxs)
80      partial <- sort.int(unique(c(idxs_lo, idxs_hi)))
81
82      # Adjust?
83      idxs_adj <- which(idxs > idxs_lo)
84      adj <- (length(idxs_adj) > 0L)
85      # Adjust
86      if (adj) {
87        idxs_hi <- idxs_hi[idxs_adj]
88        w <- (idxs - idxs_lo)[idxs_adj]
89        q_lo <- matrix(na_value, nrow = length(idxs_adj), ncol = nrow)
90        q_hi <- matrix(na_value, nrow = length(idxs_adj), ncol = nrow)
91        cols <- seq_len(ncol)
92        for (rr in seq_len(nrow)) {
93          x_rr <- .subset(x, rr, cols, drop = TRUE)
94          xp_rr <- sort.int(x_rr, partial = partial)
95	  q_rr <- .subset(xp_rr, idxs_lo)
96          q[rr,] <- q_rr
97          q_hi[,rr] <- .subset(xp_rr, idxs_hi)
98          q_lo[,rr] <- .subset(q_rr, idxs_adj)
99          # Not needed anymore
100          x_rr <- xp_rr <- NULL
101        }
102	q_adj <- (1 - w) * q_lo + w * q_hi
103	for (cc in seq_along(idxs_adj)) {
104	  q[, idxs_adj[cc]] <- q_adj[cc, , drop = TRUE]
105	}
106        # Not needed anymore
107        q_adj <- q_lo <- q_hi <- NULL
108      } else {
109        cols <- seq_len(ncol)
110        for (rr in seq_len(nrow)) {
111          x_rr <- .subset(x, rr, cols, drop = TRUE)
112          xp_rr <- sort.int(x_rr, partial = partial)
113          q[rr,] <- .subset(xp_rr, idxs_lo)
114          # Not needed anymore
115          x_rr <- xp_rr <- NULL
116        }
117      }
118
119      storage.mode(q) <- "numeric"
120    } else {
121      # For each row...
122      rows <- seq_len(nrow)
123
124      # Rows with NAs should return all NAs (so skip those)
125      if (has_na && !na.rm) rows <- rows[!na_rows]
126
127      for (kk in rows) {
128        xkk <- x[kk, ]
129        if (na.rm) xkk <- xkk[!is.na(xkk)]
130        q[kk, ] <- quantile(xkk, probs = probs, na.rm = FALSE, type = type, ...)
131      }
132    } # if (type ...)
133  }
134
135  # Add dim names
136  digits <- max(2L, getOption("digits"))
137  colnames(q) <- sprintf("%.*g%%", digits, 100 * probs)
138
139  # Preserve names attribute?
140  if (is.na(useNames) || useNames) {
141    rownames(q) <- rownames(x)
142  } else {
143    rownames(q) <- NULL
144  }
145
146  # Drop singleton dimensions?
147  if (drop) {
148    q <- drop(q)
149  }
150
151  q
152}
153
154#' @importFrom stats quantile
155#' @rdname rowQuantiles
156#' @export
157colQuantiles <- function(x, rows = NULL, cols = NULL,
158                         probs = seq(from = 0, to = 1, by = 0.25),
159                         na.rm = FALSE, type = 7L, ..., useNames = NA, drop = TRUE) {
160  # Argument 'x':
161  if (!is.matrix(x)) defunctShouldBeMatrix(x)
162  if (!is.numeric(x) && !is.integer(x) && !is.logical(x)) {
163    .Defunct(msg = sprintf("Argument 'x' is of type %s. Only 'integer', 'numeric', and 'logical' is supported.", sQuote(storage.mode(x))))  #nolint
164  }
165
166  # Argument 'probs':
167  if (anyMissing(probs)) {
168    stop(sprintf("Argument '%s' must not contain missing values", "probs"))
169  }
170  eps <- 100 * .Machine$double.eps
171  if (any((probs < -eps | probs > 1 + eps))) {
172    stop(sprintf("Argument '%s' is out of range [0-eps, 1+eps]: %g", "probs", probs))
173  }
174
175  # Apply subset
176  if (!is.null(rows) && !is.null(cols)) x <- x[rows, cols, drop = FALSE]
177  else if (!is.null(rows)) x <- x[rows, , drop = FALSE]
178  else if (!is.null(cols)) x <- x[, cols, drop = FALSE]
179
180  # Coerce?
181#  if (is.logical(x)) {
182#    storage.mode(x) <- "integer"
183#  }
184
185  # Argument 'x':
186  nrow <- nrow(x)
187  ncol <- ncol(x)
188
189  # Allocate result
190  na_value <- NA_real_
191  if (type != 7L) storage.mode(na_value) <- storage.mode(x)
192  q <- matrix(na_value, nrow = ncol, ncol = length(probs))
193
194  if (nrow > 0L && ncol > 0L) {
195    na_cols <- colAnyMissings(x)
196    has_na <- any(na_cols)
197    if (!has_na) na.rm <- FALSE
198
199    if (!has_na && type == 7L) {
200      n <- nrow
201      idxs <- 1 + (n - 1) * probs
202      idxs_lo <- floor(idxs)
203      idxs_hi <- ceiling(idxs)
204      partial <- sort.int(unique(c(idxs_lo, idxs_hi)))
205
206      # Adjust?
207      idxs_adj <- which(idxs > idxs_lo)
208      adj <- (length(idxs_adj) > 0L)
209      if (adj) {
210        idxs_hi <- idxs_hi[idxs_adj]
211        w <- (idxs - idxs_lo)[idxs_adj]
212        q_lo <- matrix(na_value, nrow = length(idxs_adj), ncol = ncol)
213        q_hi <- matrix(na_value, nrow = length(idxs_adj), ncol = ncol)
214        rows <- seq_len(nrow)
215        for (cc in seq_len(ncol)) {
216          x_cc <- .subset(x, rows, cc, drop = TRUE)
217          xp_cc <- sort.int(x_cc, partial = partial)
218	  q_cc <- .subset(xp_cc, idxs_lo)
219          q[cc,] <- q_cc
220          q_hi[,cc] <- .subset(xp_cc, idxs_hi)
221          q_lo[,cc] <- .subset(q_cc, idxs_adj)
222          # Not needed anymore
223          x_cc <- xp_cc <- NULL
224        }
225	q_adj <- (1 - w) * q_lo + w * q_hi
226	for (cc in seq_along(idxs_adj)) {
227	  q[, idxs_adj[cc]] <- q_adj[cc, , drop = TRUE]
228	}
229        # Not needed anymore
230        q_adj <- q_lo <- q_hi <- NULL
231      } else {
232        rows <- seq_len(nrow)
233        for (cc in seq_len(ncol)) {
234          x_cc <- .subset(x, rows, cc, drop = TRUE)
235          xp_cc <- sort.int(x_cc, partial = partial)
236          q[cc,] <- .subset(xp_cc, idxs_lo)
237          # Not needed anymore
238          x_cc <- xp_cc <- NULL
239        }
240      }
241
242    } else {
243      # For each column...
244      cols <- seq_len(ncol)
245
246      # Columns with NAs should return all NAs (so skip those)
247      if (has_na && !na.rm) cols <- cols[!na_cols]
248
249      for (kk in cols) {
250        xkk <- x[, kk]
251        if (na.rm) xkk <- xkk[!is.na(xkk)]
252        q[kk, ] <- quantile(xkk, probs = probs, na.rm = FALSE, type = type, ...)
253      }
254    } # if (type ...)
255  }
256
257  # Add dim names
258  digits <- max(2L, getOption("digits"))
259  colnames(q) <- sprintf("%.*g%%", digits, 100 * probs)
260
261  # Preserve names attribute?
262  if (is.na(useNames) || useNames) {
263    rownames(q) <- colnames(x)
264  } else {
265    rownames(q) <- NULL
266  }
267
268  # Drop singleton dimensions?
269  if (drop) {
270    q <- drop(q)
271  }
272
273  q
274}
275