1library("matrixStats")
2
3# - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4# Consistency checks
5# - - - - - - - - - - - - - - - - - - - - - - - - - - - -
6set.seed(1)
7
8mean2_R <- function(x, na.rm = FALSE, idxs = NULL) {
9  if (is.null(idxs)) {
10    mean(x, na.rm = na.rm)
11  } else {
12    mean(x[idxs], na.rm = na.rm)
13  }
14} # mean2_R()
15
16
17cat("Consistency checks:\n")
18for (kk in 1:20) {
19  cat("Random test #", kk, "\n", sep = "")
20
21  # Simulate data in a matrix of any shape
22  n <- sample(100L, size = 1L)
23  x <- rnorm(n, sd = 100)
24
25  # Add NAs?
26  if ((kk %% 4) %in% c(3, 0)) {
27    cat("Adding NAs\n")
28    nna <- sample(n, size = 1L)
29    na_values <- c(NA_real_, NaN)
30    t <- sample(na_values, size = nna, replace = TRUE)
31    x[sample(length(x), size = nna)] <- t
32  }
33
34  # Integer or double?
35  if ((kk %% 4) %in% c(2, 0)) {
36    cat("Coercing to integers\n")
37    storage.mode(x) <- "integer"
38  }
39
40  na.rm <- sample(c(TRUE, FALSE), size = 1L)
41
42  # Sum over all
43  y0 <- mean2_R(x, na.rm = na.rm)
44  y1 <- mean2(x, na.rm = na.rm)
45  stopifnot(all.equal(y1, y0))
46
47  # Sum over subset
48  nidxs <- sample(n, size = 1L)
49  idxs <- sample(n, size = nidxs)
50  y0 <- mean2_R(x, na.rm = na.rm, idxs = idxs)
51  y1 <- mean2(x, na.rm = na.rm, idxs = idxs)
52  stopifnot(all.equal(y1, y0))
53
54  if (storage.mode(x) == "integer") {
55    storage.mode(x) <- "logical"
56
57    y0 <- mean2_R(x, na.rm = na.rm)
58    y1 <- mean2(x, na.rm = na.rm)
59    stopifnot(all.equal(y1, y0))
60
61    y0 <- mean2_R(x, na.rm = na.rm, idxs = idxs)
62    y1 <- mean2(x, na.rm = na.rm, idxs = idxs)
63    stopifnot(all.equal(y1, y0))
64  }
65} # for (kk ...)
66
67
68# - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69# Special cases
70# - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71for (na.rm in c(FALSE, TRUE)) {
72  # Averaging over zero elements (integers)
73  x <- integer(0)
74  s1 <- mean(x, na.rm = na.rm)
75  s2 <- mean2(x, na.rm = na.rm)
76  stopifnot(identical(s1, s2))
77
78  x <- 1:5
79  idxs <- integer(0)
80  s1 <- mean(x[idxs], na.rm = na.rm)
81  s2 <- mean2(x, idxs = idxs, na.rm = na.rm)
82  stopifnot(identical(s1, s2))
83
84  # Averaging over NA_integer_:s
85  x <- rep(NA_integer_, times = 5L)
86  s1 <- mean(x, na.rm = na.rm)
87  s2 <- mean2(x, na.rm = na.rm)
88  stopifnot(identical(s1, s2))
89
90  x <- rep(NA_integer_, times = 5L)
91  idxs <- 1:3
92  s1 <- mean(x[idxs], na.rm = na.rm)
93  s2 <- mean2(x, idxs = idxs, na.rm = na.rm)
94  stopifnot(identical(s1, s2))
95
96
97  # Averaging over zero elements (doubles)
98  x <- double(0)
99  s1 <- mean(x)
100  s2 <- mean2(x)
101  stopifnot(identical(s1, s2))
102
103  x <- as.double(1:10)
104  idxs <- integer(0)
105  s1 <- mean(x[idxs])
106  s2 <- mean2(x, idxs = idxs)
107  stopifnot(identical(s1, s2))
108
109  # Averaging over NA_real_:s
110  x <- rep(NA_real_, times = 5L)
111  s1 <- mean(x, na.rm = na.rm)
112  s2 <- mean2(x, na.rm = na.rm)
113  stopifnot(identical(s1, s2))
114
115  x <- rep(NA_real_, times = 5L)
116  idxs <- 1:3
117  s1 <- mean(x[idxs], na.rm = na.rm)
118  s2 <- mean2(x, idxs = idxs, na.rm = na.rm)
119  stopifnot(identical(s1, s2))
120
121  # Averaging over -Inf:s
122  x <- rep(-Inf, times = 3L)
123  s1 <- mean(x, na.rm = na.rm)
124  s2 <- mean2(x, na.rm = na.rm)
125  stopifnot(identical(s1, s2))
126
127  # Averaging over +Inf:s
128  x <- rep(+Inf, times = 3L)
129  s1 <- mean(x, na.rm = na.rm)
130  s2 <- mean2(x, na.rm = na.rm)
131  stopifnot(identical(s1, s2))
132
133  # Averaging over mix of -Inf:s and +Inf:s
134  x <- rep(c(-Inf, +Inf), times = 3L)
135  s1 <- mean(x, na.rm = na.rm)
136  s2 <- mean2(x, na.rm = na.rm)
137  stopifnot(identical(s1, s2))
138
139  # Averaging over mix of -Inf:s and +Inf:s and numerics
140  x <- rep(c(-Inf, +Inf, 3.14), times = 2L)
141  s1 <- mean(x, na.rm = na.rm)
142  s2 <- mean2(x, na.rm = na.rm)
143  stopifnot(identical(s1, s2))
144
145  # Averaging over mix of NaN, NA, +Inf, and numerics
146  x <- c(NaN, NA, +Inf, 3.14)
147  s1 <- mean(x, na.rm = na.rm)
148  s2 <- mean2(x, na.rm = na.rm)
149  if (na.rm) {
150    stopifnot(identical(s2, s1))
151  } else {
152    stopifnot(is.na(s1), is.na(s2))
153    ## NOTE, due to compiler optimization, it is not guaranteed that NA is
154    ## returned here (as one would expect). NaN might very well be returned,
155    ## when both NA and NaN are involved.  This is an accepted feature in R,
156    ## which is documented in help("is.nan").  See also
157    ## https://stat.ethz.ch/pipermail/r-devel/2017-April/074009.html.
158    ## Thus, we cannot guarantee that s1 is identical to s0.
159  }
160
161  # Averaging over mix of NaN, NA_real_, +Inf, and numerics
162  x <- c(NA_real_, NaN, +Inf, 3.14)
163  s1 <- mean(x, na.rm = na.rm)
164  s2 <- mean2(x, na.rm = na.rm)
165  if (na.rm) {
166    stopifnot(identical(s2, s1))
167  } else {
168    stopifnot(is.na(s1), is.na(s2))
169    ## NOTE, due to compiler optimization, it is not guaranteed that NA is
170    ## returned here (as one would expect). NaN might very well be returned,
171    ## when both NA and NaN are involved.  This is an accepted feature in R,
172    ## which is documented in help("is.nan").  See also
173    ## https://stat.ethz.ch/pipermail/r-devel/2017-April/074009.html.
174    ## Thus, we cannot guarantee that s1 is identical to s0.
175  }
176}
177
178
179# - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180# Argument 'idxs'
181# - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182x <- 1:5
183idxs_list <- list(
184  integer = 1:3,
185  double  = as.double(1:3),
186  logical = (x <= 3)
187)
188
189for (idxs in idxs_list) {
190  cat("idxs:\n")
191  str(idxs)
192  s1 <- mean(x[idxs], na.rm = TRUE)
193  s2 <- mean2(x, idxs = idxs, na.rm = TRUE)
194  stopifnot(identical(s1, s2))
195}
196