context("fmean")

if(!is.null(attributes(identical(FALSE, TRUE)))) stop("OECD label issue")

bmean <- base::mean
bsum <- base::sum

# rm(list = ls())
set.seed(101)
x <- rnorm(100)
w <- abs(100*rnorm(100))
wdat <- abs(100*rnorm(32))
xNA <- x
wNA <- w
wdatNA <- wdat
xNA[sample.int(100,20)] <- NA
wNA[sample.int(100,20)] <- NA
wdatNA[sample.int(32, 5)] <- NA
f <- as.factor(sample.int(10, 100, TRUE))
g <- GRP(mtcars, ~ cyl + vs + am)
gf <- as_factor_GRP(g)
mtcNA <- na_insert(mtcars)
mtcNA[27,1] <- NA # single group NA !!
m <- as.matrix(mtcars)
mNA <- as.matrix(mtcNA)
mNAc <- mNA
storage.mode(mNAc) <- "character"

wmean <- function(x, w, na.rm = FALSE) {
  if(na.rm) {
    cc <- complete.cases(x, w)
    x <- x[cc]
    w <- w[cc]
  }
  bsum(x*w)/bsum(w)
}


for (nth in 1:2) {

  if(nth == 2L) {
    if(Sys.getenv("OMP") == "TRUE") {
      fmean <- function(x, ...) collapse::fmean(x, ..., nthreads = 2L)
    } else break
  }

test_that("fmean performs like base::mean", {
  expect_equal(fmean(NA), bmean(NA))
  expect_equal(fmean(NA, na.rm = FALSE), bmean(NA))
  expect_equal(fmean(1), bmean(1, na.rm = TRUE))
  expect_equal(fmean(1:3), bmean(1:3, na.rm = TRUE))
  expect_equal(fmean(-1:1), bmean(-1:1, na.rm = TRUE))
  expect_equal(fmean(1, na.rm = FALSE), bmean(1))
  expect_equal(fmean(1:3, na.rm = FALSE), bmean(1:3))
  expect_equal(fmean(-1:1, na.rm = FALSE), bmean(-1:1))
  expect_equal(fmean(x), bmean(x, na.rm = TRUE))
  expect_equal(fmean(x, na.rm = FALSE), bmean(x))
  expect_equal(fmean(xNA, na.rm = FALSE), bmean(xNA))
  expect_equal(fmean(xNA), bmean(xNA, na.rm = TRUE))
  expect_equal(fmean(mtcars), fmean(m))
  expect_equal(fmean(m), dapply(m, bmean, na.rm = TRUE))
  expect_equal(fmean(m, na.rm = FALSE), dapply(m, bmean))
  expect_equal(fmean(mNA, na.rm = FALSE), dapply(mNA, bmean))
  expect_equal(fmean(mNA), dapply(mNA, bmean, na.rm = TRUE))
  expect_equal(fmean(mtcars), dapply(mtcars, bmean, na.rm = TRUE))
  expect_equal(fmean(mtcars, na.rm = FALSE), dapply(mtcars, bmean))
  expect_equal(fmean(mtcNA, na.rm = FALSE), dapply(mtcNA, bmean))
  expect_equal(fmean(mtcNA), dapply(mtcNA, bmean, na.rm = TRUE))
  expect_equal(fmean(x, f), BY(x, f, bmean, na.rm = TRUE))
  expect_equal(fmean(x, f, na.rm = FALSE), BY(x, f, bmean))
  expect_equal(fmean(xNA, f, na.rm = FALSE), BY(xNA, f, bmean))
  expect_equal(fmean(xNA, f), BY(xNA, f, bmean, na.rm = TRUE))
  expect_equal(fmean(m, g), BY(m, g, bmean, na.rm = TRUE))
  expect_equal(fmean(m, g, na.rm = FALSE), BY(m, g, bmean))
  expect_equal(fmean(mNA, g, na.rm = FALSE), BY(mNA, g, bmean))
  expect_equal(fmean(mNA, g), BY(mNA, g, bmean, na.rm = TRUE))
  expect_equal(fmean(mtcars, g), BY(mtcars, g, bmean, na.rm = TRUE))
  expect_equal(fmean(mtcars, g, na.rm = FALSE), BY(mtcars, g, bmean))
  expect_equal(fmean(mtcNA, g, na.rm = FALSE), BY(mtcNA, g, bmean))
  expect_equal(fmean(mtcNA, g), BY(mtcNA, g, bmean, na.rm = TRUE))
})

test_that("fmean with weights performs as intended (unbiased)", {
  expect_equal(fmean(c(2,2,4,5,5,5)), fmean(c(2,4,5), w = c(2,1,3)))
  expect_equal(fmean(c(2,2,4,5,5,5), na.rm = FALSE), fmean(c(2,4,5), w = c(2,1,3), na.rm = FALSE))
  expect_equal(fmean(c(2.456,2.456,4.123,5.009,5.009,5.009)), fmean(c(2.456,4.123,5.009), w = c(2,1,3)))
  expect_equal(fmean(c(2.456,2.456,4.123,5.009,5.009,5.009), na.rm = FALSE), fmean(c(2.456,4.123,5.009), w = c(2,1,3), na.rm = FALSE))
  expect_equal(fmean(c(2,2,NA,5,5,5)), fmean(c(2,NA,5), w = c(2,1,3)))
  expect_equal(fmean(c(2,2,NA,5,5,5), na.rm = FALSE), fmean(c(2,NA,5), w = c(2,1,3), na.rm = FALSE))
  expect_equal(fmean(c(2,2,NA,5,5,5)), fmean(c(2,4,5), w = c(2,NA,3)))
  expect_equal(fmean(c(2,2,NA,5,5,5), na.rm = FALSE), fmean(c(2,4,5), w = c(2,NA,3), na.rm = FALSE))
  expect_equal(fmean(c(NA,NA,4.123,5.009,5.009,5.009)), fmean(c(NA,4.123,5.009), w = c(2,1,3)))
  expect_equal(fmean(c(NA,NA,4.123,5.009,5.009,5.009), na.rm = FALSE), fmean(c(NA,4.123,5.009), w = c(2,1,3), na.rm = FALSE))
  expect_equal(fmean(c(NA,NA,4.123,5.009,5.009,5.009)), fmean(c(2.456,4.123,5.009), w = c(NA,1,3)))
  expect_equal(fmean(c(NA,NA,4.123,5.009,5.009,5.009), na.rm = FALSE), fmean(c(2.456,4.123,5.009), w = c(NA,1,3), na.rm = FALSE))
  f <- as.factor(rep(1:2, each = 6)); fs <- as.factor(rep(1:2, each = 3))
  v <- c(2,2,4,5,5,5,2,2,4,5,5,5); vs <- c(2,4,5,2,4,5); w <- c(2,1,3,2,1,3)
  v2 <- c(2.456,2.456,4.123,5.009,5.009,5.009,2.456,2.456,4.123,5.009,5.009,5.009); v2s <- c(2.456,4.123,5.009,2.456,4.123,5.009)
  expect_equal(fmean(v, f), fmean(vs, fs, w))
  expect_equal(fmean(v, f, na.rm = FALSE), fmean(vs, fs, w, na.rm = FALSE))
  expect_equal(fmean(v2, f), fmean(v2s, fs, w))
  expect_equal(fmean(v2, f, na.rm = FALSE), fmean(v2s, fs, w, na.rm = FALSE))
  v[c(3,9)] <- NA; vs[c(2,5)] <- NA
  expect_equal(fmean(v, f), fmean(vs, fs, w))
  expect_equal(fmean(v, f, na.rm = FALSE), fmean(vs, fs, w, na.rm = FALSE))
  vs[c(2,5)] <- 4; w[c(2,5)] <- NA
  expect_equal(fmean(v, f), fmean(vs, fs, w))
  expect_equal(fmean(v, f, na.rm = FALSE), fmean(vs, fs, w, na.rm = FALSE))
  w[c(2,5)] <- 1; v2[c(1:2,7:8)] <- NA; v2s[c(1,4)] <- NA
  expect_equal(fmean(v2, f), fmean(v2s, fs, w))
  expect_equal(fmean(v2, f, na.rm = FALSE), fmean(v2s, fs, w, na.rm = FALSE))
  v2s[c(1,4)] <- 2.456; w[c(1,4)] <- NA
  expect_equal(fmean(v2, f), fmean(v2s, fs, w))
  expect_equal(fmean(v2, f, na.rm = FALSE), fmean(v2s, fs, w, na.rm = FALSE))
})

test_that("fmean performs like fmean with weights all equal", {
  expect_equal(fmean(NA), fmean(NA, w = 0.99999999))
  expect_equal(fmean(NA, na.rm = FALSE), fmean(NA, w = 2.946, na.rm = FALSE))
  expect_equal(fmean(1), fmean(1, w = 3))
  expect_equal(fmean(1:3), fmean(1:3, w = rep(0.999,3)))
  expect_equal(fmean(-1:1), fmean(-1:1, w = rep(4.2,3)))
  expect_equal(fmean(1, na.rm = FALSE), fmean(1, w = 5, na.rm = FALSE))
  expect_equal(fmean(1:3, na.rm = FALSE), fmean(1:3, w = rep(1.44565, 3), na.rm = FALSE))
  expect_equal(fmean(-1:1, na.rm = FALSE), fmean(-1:1, w = rep(1.44565, 3), na.rm = FALSE))
  expect_equal(fmean(x), fmean(x, w = rep(1,100)))
  expect_equal(fmean(x, na.rm = FALSE), fmean(x, w = rep(1.44565, 100), na.rm = FALSE))
  expect_equal(fmean(xNA, na.rm = FALSE), fmean(xNA, w = rep(4.676587, 100), na.rm = FALSE))
  expect_equal(fmean(xNA), fmean(xNA, w = rep(4.676587, 100)))
  expect_equal(fmean(m), fmean(m, w = rep(6587.3454, 32)))
  expect_equal(fmean(m, na.rm = FALSE), fmean(m, w = rep(6587.3454, 32), na.rm = FALSE))
  expect_equal(fmean(mNA, na.rm = FALSE), fmean(mNA, w = rep(6587.3454, 32), na.rm = FALSE))
  expect_equal(fmean(mNA), fmean(mNA, w = rep(6587.3454, 32)))
  expect_equal(fmean(mtcars), fmean(mtcars, w = rep(6787.3454, 32)))
  expect_equal(fmean(mtcars, na.rm = FALSE), fmean(mtcars, w = rep(6787.3454, 32), na.rm = FALSE))
  expect_equal(fmean(mtcNA, na.rm = FALSE), fmean(mtcNA, w = rep(6787.3454, 32), na.rm = FALSE))
  expect_equal(fmean(mtcNA), fmean(mtcNA, w = rep(6787.3454, 32)))
  expect_equal(fmean(x, f), fmean(x, f, rep(546.78,100)))
  expect_equal(fmean(x, f, na.rm = FALSE), fmean(x, f, rep(5.88,100), na.rm = FALSE))
  expect_equal(fmean(xNA, f, na.rm = FALSE), fmean(xNA, f, rep(52.7,100), na.rm = FALSE))
  expect_equal(fmean(xNA, f), fmean(xNA, f, rep(5997456,100)))
  expect_equal(fmean(m, g), fmean(m, g, rep(546.78,32)))
  expect_equal(fmean(m, g, na.rm = FALSE), fmean(m, g, rep(0.0001,32), na.rm = FALSE))
  expect_equal(fmean(mNA, g, na.rm = FALSE), fmean(mNA, g, rep(5.7,32), na.rm = FALSE))
  expect_equal(fmean(mNA, g), fmean(mNA, g, rep(1.1,32)))
  expect_equal(fmean(mtcars, g), fmean(mtcars, g, rep(53,32)))
  expect_equal(fmean(mtcars, g, na.rm = FALSE), fmean(mtcars, g, rep(546.78,32), na.rm = FALSE))
  expect_equal(fmean(mtcNA, g, na.rm = FALSE), fmean(mtcNA, g, rep(0.999999,32), na.rm = FALSE))
  expect_equal(fmean(mtcNA, g), fmean(mtcNA, g, rep(999.9999,32)))
})

test_that("fmean with weights performs like wmean (defined above)", {
  # complete weights
  expect_equal(fmean(NA, w = 1), wmean(NA, 1))
  expect_equal(fmean(NA, w = 1, na.rm = FALSE), wmean(NA, 1))
  expect_equal(fmean(1, w = 1), wmean(1, w = 1))
  expect_equal(fmean(1:3, w = 1:3), wmean(1:3, 1:3))
  expect_equal(fmean(-1:1, w = 1:3), wmean(-1:1, 1:3))
  expect_equal(fmean(1, w = 1, na.rm = FALSE), wmean(1, 1))
  expect_equal(fmean(1:3, w = c(0.99,3454,1.111), na.rm = FALSE), wmean(1:3, c(0.99,3454,1.111)))
  expect_equal(fmean(-1:1, w = 1:3, na.rm = FALSE), wmean(-1:1, 1:3))
  expect_equal(fmean(x, w = w), wmean(x, w))
  expect_equal(fmean(x, w = w, na.rm = FALSE), wmean(x, w))
  expect_equal(fmean(xNA, w = w, na.rm = FALSE), wmean(xNA, w))
  expect_equal(fmean(xNA, w = w), wmean(xNA, w, na.rm = TRUE))
  expect_equal(fmean(mtcars, w = wdat), fmean(m, w = wdat))
  expect_equal(fmean(m, w = wdat), dapply(m, wmean, wdat, na.rm = TRUE))
  expect_equal(fmean(m, w = wdat, na.rm = FALSE), dapply(m, wmean, wdat))
  expect_equal(fmean(mNA, w = wdat, na.rm = FALSE), dapply(mNA, wmean, wdat))
  expect_equal(fmean(mNA, w = wdat), dapply(mNA, wmean, wdat, na.rm = TRUE))
  expect_equal(fmean(mtcars, w = wdat), dapply(mtcars, wmean, wdat, na.rm = TRUE))
  expect_equal(fmean(mtcars, w = wdat, na.rm = FALSE), dapply(mtcars, wmean, wdat))
  expect_equal(fmean(mtcNA, w = wdat, na.rm = FALSE), dapply(mtcNA, wmean, wdat))
  expect_equal(fmean(mtcNA, w = wdat), dapply(mtcNA, wmean, wdat, na.rm = TRUE))
  expect_equal(fmean(x, f, w), BY(x, f, wmean, w))
  expect_equal(fmean(x, f, w, na.rm = FALSE), BY(x, f, wmean, w))
  expect_equal(fmean(xNA, f, w, na.rm = FALSE), BY(xNA, f, wmean, w))
  expect_equal(fmean(xNA, f, w), BY(xNA, f, wmean, w, na.rm = TRUE))
  expect_equal(fmean(m, g, wdat), BY(m, gf, wmean, wdat))
  expect_equal(fmean(m, g, wdat, na.rm = FALSE), BY(m, gf, wmean, wdat))
  expect_equal(fmean(mNA, g, wdat, na.rm = FALSE),  BY(mNA, gf, wmean, wdat))
  expect_equal(fmean(mNA, g, wdat), BY(mNA, gf, wmean, wdat, na.rm = TRUE))
  expect_equal(fmean(mtcars, g, wdat), BY(mtcars, gf, wmean, wdat))
  expect_equal(fmean(mtcars, g, wdat, na.rm = FALSE), BY(mtcars, gf, wmean, wdat))
  expect_equal(fmean(mtcNA, g, wdat, na.rm = FALSE), BY(mtcNA, gf, wmean, wdat))
  expect_equal(fmean(mtcNA, g, wdat), BY(mtcNA, gf, wmean, wdat, na.rm = TRUE))
  # missing weights
  expect_equal(fmean(NA, w = NA), wmean(NA, NA))
  expect_equal(fmean(NA, w = NA, na.rm = FALSE), wmean(NA, NA))
  expect_equal(fmean(1, w = NA), wmean(1, w = NA))
  expect_equal(fmean(1:3, w = c(NA,1:2)), wmean(1:3, c(NA,1:2), na.rm = TRUE))
  expect_equal(fmean(-1:1, w = c(NA,1:2)), wmean(-1:1, c(NA,1:2), na.rm = TRUE))
  expect_equal(fmean(1, w = NA, na.rm = FALSE), wmean(1, NA))
  expect_equal(fmean(1:3, w = c(NA,1:2), na.rm = FALSE), wmean(1:3, c(NA,1:2)))
  expect_equal(fmean(-1:1, w = c(NA,1:2), na.rm = FALSE), wmean(-1:1, c(NA,1:2)))
  expect_equal(fmean(x, w = wNA), wmean(x, wNA, na.rm = TRUE))
  expect_equal(fmean(x, w = wNA, na.rm = FALSE), wmean(x, wNA))
  expect_equal(fmean(xNA, w = wNA, na.rm = FALSE), wmean(xNA, wNA))
  expect_equal(fmean(xNA, w = wNA), wmean(xNA, wNA, na.rm = TRUE))
  expect_equal(fmean(mtcars, w = wdatNA), fmean(m, w = wdatNA))
  expect_equal(fmean(m, w = wdatNA), dapply(m, wmean, wdatNA, na.rm = TRUE))
  expect_equal(fmean(m, w = wdatNA, na.rm = FALSE), dapply(m, wmean, wdatNA))
  expect_equal(fmean(mNA, w = wdatNA, na.rm = FALSE), dapply(mNA, wmean, wdatNA))
  expect_equal(fmean(mNA, w = wdatNA), dapply(mNA, wmean, wdatNA, na.rm = TRUE))
  expect_equal(fmean(mtcars, w = wdatNA), dapply(mtcars, wmean, wdatNA, na.rm = TRUE))
  expect_equal(fmean(mtcars, w = wdatNA, na.rm = FALSE), dapply(mtcars, wmean, wdatNA))
  expect_equal(fmean(mtcNA, w = wdatNA, na.rm = FALSE), dapply(mtcNA, wmean, wdatNA))
  expect_equal(fmean(mtcNA, w = wdatNA), dapply(mtcNA, wmean, wdatNA, na.rm = TRUE))
  expect_equal(fmean(x, f, wNA), BY(x, f, wmean, wNA, na.rm = TRUE))
  expect_equal(fmean(x, f, wNA, na.rm = FALSE), BY(x, f, wmean, wNA))
  expect_equal(fmean(xNA, f, wNA, na.rm = FALSE), BY(xNA, f, wmean, wNA))
  expect_equal(fmean(xNA, f, wNA), BY(xNA, f, wmean, wNA, na.rm = TRUE))
  expect_equal(fmean(m, g, wdatNA), BY(m, gf, wmean, wdatNA, na.rm = TRUE))
  expect_equal(fmean(m, g, wdatNA, na.rm = FALSE), BY(m, gf, wmean, wdatNA))
  expect_equal(fmean(mNA, g, wdatNA, na.rm = FALSE),  BY(mNA, gf, wmean, wdatNA))
  expect_equal(fmean(mNA, g, wdatNA), BY(mNA, gf, wmean, wdatNA, na.rm = TRUE))
  expect_equal(fmean(mtcars, g, wdatNA), BY(mtcars, gf, wmean, wdatNA, na.rm = TRUE))
  expect_equal(fmean(mtcars, g, wdatNA, na.rm = FALSE), BY(mtcars, gf, wmean, wdatNA))
  expect_equal(fmean(mtcNA, g, wdatNA, na.rm = FALSE), BY(mtcNA, gf, wmean, wdatNA))
  expect_equal(fmean(mtcNA, g, wdatNA), BY(mtcNA, gf, wmean, wdatNA, na.rm = TRUE))
})

test_that("fmean performs numerically stable", {
  expect_true(all_obj_equal(replicate(50, fmean(1), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(NA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(NA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, f), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, f, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, f, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, f), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, g), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, g, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, g, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, g), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, g), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, g, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, g, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, g), simplify = FALSE)))
})

test_that("fmean with complete weights performs numerically stable", {
  expect_true(all_obj_equal(replicate(50, fmean(1, w = 1), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(NA, w = 1), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(NA, w = 1, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, w = w), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, w = w, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, w = w, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, w = w), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, w = wdat), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, w = wdat, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, w = wdat, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, w = wdat), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, w = wdat), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, w = wdat, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, w = wdat, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, w = wdat), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, f, w), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, f, w, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, f, w, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, f, w), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, g, wdat), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, g, wdat, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, g, wdat, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, g, wdat), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, g, wdat), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, g, wdat, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, g, wdat, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, g, wdat), simplify = FALSE)))
})


test_that("fmean with missing weights performs numerically stable", {
  expect_true(all_obj_equal(replicate(50, fmean(1, w = NA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(NA, w = NA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(NA, w = NA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, w = wNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, w = wNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, w = wNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, w = wNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, w = wdatNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, w = wdatNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, w = wdatNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, w = wdatNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, w = wdatNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, w = wdatNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, w = wdatNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, w = wdatNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, f, wNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(x, f, wNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, f, wNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(xNA, f, wNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, g, wdatNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(m, g, wdatNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, g, wdatNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mNA, g, wdatNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, g, wdatNA), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcars, g, wdatNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, g, wdatNA, na.rm = FALSE), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fmean(mtcNA, g, wdatNA), simplify = FALSE)))
})

test_that("fmean handles special values in the right way", {
  expect_equal(fmean(NA), NA_real_)
  expect_equal(fmean(NaN), NaN)
  expect_equal(fmean(Inf), Inf)
  expect_equal(fmean(-Inf), -Inf)
  expect_equal(fmean(TRUE), 1)
  expect_equal(fmean(FALSE), 0)
  expect_equal(fmean(NA, na.rm = FALSE), NA_real_)
  expect_equal(fmean(NaN, na.rm = FALSE), NaN)
  expect_equal(fmean(Inf, na.rm = FALSE), Inf)
  expect_equal(fmean(-Inf, na.rm = FALSE), -Inf)
  expect_equal(fmean(TRUE, na.rm = FALSE), 1)
  expect_equal(fmean(FALSE, na.rm = FALSE), 0)
  expect_equal(fmean(c(1,NA)), 1)
  expect_equal(fmean(c(1,NaN)), 1)
  expect_equal(fmean(c(1,Inf)), Inf)
  expect_equal(fmean(c(1,-Inf)), -Inf)
  expect_equal(fmean(c(FALSE,TRUE)), 0.5)
  expect_equal(fmean(c(FALSE,FALSE)), 0)
  expect_equal(fmean(c(1,Inf), na.rm = FALSE), Inf)
  expect_equal(fmean(c(1,-Inf), na.rm = FALSE), -Inf)
  expect_equal(fmean(c(FALSE,TRUE), na.rm = FALSE), 0.5)
  expect_equal(fmean(c(FALSE,FALSE), na.rm = FALSE), 0)
})

test_that("fmean with weights handles special values in the right way", {
  expect_equal(fmean(NA, w = 1), NA_real_)
  expect_equal(fmean(NaN, w = 1), NaN)
  expect_equal(fmean(Inf, w = 1), Inf)
  expect_equal(fmean(-Inf, w = 1), -Inf)
  expect_equal(fmean(TRUE, w = 1), 1)
  expect_equal(fmean(FALSE, w = 1), 0)
  expect_equal(fmean(NA, w = 1, na.rm = FALSE), NA_real_)
  expect_equal(fmean(NaN, w = 1, na.rm = FALSE), NaN)
  expect_equal(fmean(Inf, w = 1, na.rm = FALSE), Inf)
  expect_equal(fmean(-Inf, w = 1, na.rm = FALSE), -Inf)
  expect_equal(fmean(TRUE, w = 1, na.rm = FALSE), 1)
  expect_equal(fmean(FALSE, w = 1, na.rm = FALSE), 0)
  expect_equal(fmean(NA, w = NA), NA_real_)
  expect_equal(fmean(NaN, w = NA), NA_real_)
  expect_equal(fmean(Inf, w = NA), NA_real_)
  expect_equal(fmean(-Inf, w = NA), NA_real_)
  expect_equal(fmean(TRUE, w = NA), NA_real_)
  expect_equal(fmean(FALSE, w = NA), NA_real_)
  expect_equal(fmean(NA, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fmean(NaN, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fmean(Inf, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fmean(-Inf, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fmean(TRUE, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fmean(FALSE, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fmean(1:3, w = c(1,Inf,3)), NaN)
  expect_equal(fmean(1:3, w = c(1,-Inf,3)), NaN)
  expect_equal(fmean(1:3, w = c(1,Inf,3), na.rm = FALSE), NaN)
  expect_equal(fmean(1:3, w = c(1,-Inf,3), na.rm = FALSE), NaN)
})

test_that("fmean produces errors for wrong input", {
  expect_error(fmean("a"))
  expect_error(fmean(NA_character_))
  expect_error(fmean(mNAc))
  expect_error(fmean(mNAc, f))
  expect_error(fmean(1:2,1:3))
  expect_error(fmean(m,1:31))
  expect_error(fmean(mtcars,1:31))
  expect_error(fmean(mtcars, w = 1:31))
  expect_error(fmean("a", w = 1))
  expect_error(fmean(1:2, w = 1:3))
  expect_error(fmean(NA_character_, w = 1))
  expect_error(fmean(mNAc, w = wdat))
  expect_error(fmean(mNAc, f, wdat))
  expect_error(fmean(mNA, w = 1:33))
  expect_error(fmean(1:2,1:2, 1:3))
  expect_error(fmean(m,1:32,1:20))
  expect_error(fmean(mtcars,1:32,1:10))
  expect_error(fmean(1:2, w = c("a","b")))
  expect_error(fmean(wlddev))
  expect_error(fmean(wlddev, w = wlddev$year))
  expect_error(fmean(wlddev, wlddev$iso3c))
  expect_error(fmean(wlddev, wlddev$iso3c, wlddev$year))
})

}