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