context("fbetween / B and fwithin / W")

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

# 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(rep(1:10, each = 10))
g <- as.factor(rep(c(1,2,2,3,3,3,4,4,4,4,5,5,5,5,5,7,7,7,7,7,7,7,10,10,10,10,10,10,10,10,10,10)))
mtcNA <- na_insert(mtcars)
mtcNA[1,1] <- NA # single group NA !!
m <- as.matrix(mtcars)
mNA <- as.matrix(mtcNA)
mNAc <- mNA
storage.mode(mNAc) <- "character"

# x = rnorm(1e7)
# xNA = x
# xNA[sample.int(1e7,1e6)] <- NA
# w = abs(100*rnorm(1e7))
# wNA = w
# wNA[sample.int(1e7,1e6)] <- NA

# microbenchmark(fwithin(xNA), fbetween(xNA), fbetween(xNA, w = w), fwithin(xNA, w = w), fbetween(xNA, w = wNA), fwithin(xNA, w = wNA))
# Unit: milliseconds
# expr      min       lq      mean   median       uq      max neval  cld
# fwithin(xNA) 59.89809 61.45215  81.20188 63.21997 65.99563 303.5464   100 a
# fbetween(xNA) 71.32829 73.00953  86.06850 74.51227 77.79108 275.6274   100 ab
# fbetween(xNA, w = w) 81.95167 84.85050 106.61714 86.65870 90.92104 314.8245   100   cd
# fwithin(xNA, w = w) 71.24841 73.72264  88.08572 75.32935 80.46232 279.5597   100 a c
# fbetween(xNA, w = wNA) 90.99712 93.71455 107.38818 95.91545 98.16989 328.8951   100    d
# fwithin(xNA, w = wNA) 80.13678 83.62511 103.55614 86.22361 93.18352 301.7070   100  bcd

bsum <- base::sum

between <- function(x, na.rm = FALSE) {
  if(!na.rm) return(ave(x))
  cc <- !is.na(x)
  x[cc] <- ave(x[cc])
  return(x)
}
within <- function(x, na.rm = FALSE, mean = 0) {
  if(!na.rm) return(x - ave(x) + mean)
  cc <- !is.na(x)
  m <- bsum(x[cc]) / bsum(cc)
  return(x - m + mean)
}

# NOTE: This is what fbetween and fwithin currently do: If missing values, compute weighted mean on available obs, and center x using it. But don't insert additional missing values in x for missing weights ..
wbetween <- function(x, w, na.rm = FALSE) {
  if(na.rm) {
    xcc <- !is.na(x)
    cc <- xcc & !is.na(w)
    w <- w[cc]
    wm <- bsum(w * x[cc]) / bsum(w)
    x[xcc] <- rep(wm, bsum(xcc))
    return(x)
  } else {
    wm <- bsum(w * x) / bsum(w)
    return(rep(wm, length(x)))
  }
}
wwithin <- function(x, w, na.rm = FALSE, mean = 0) {
  if(na.rm) {
    cc <- complete.cases(x, w)
    w <- w[cc]
    wm <- bsum(w * x[cc]) / bsum(w)
  } else wm <- bsum(w * x) / bsum(w)
    return(x - wm + mean)
}


# fbetween

test_that("fbetween performs like between", {
  expect_equal(fbetween(NA), as.double(between(NA)))
  expect_equal(fbetween(NA, na.rm = FALSE), as.double(between(NA)))
  expect_equal(fbetween(1), between(1, na.rm = TRUE))
  expect_equal(fbetween(1:3), between(1:3, na.rm = TRUE))
  expect_equal(fbetween(-1:1), between(-1:1, na.rm = TRUE))
  expect_equal(fbetween(1, na.rm = FALSE), between(1))
  expect_equal(fbetween(1:3, na.rm = FALSE), between(1:3))
  expect_equal(fbetween(-1:1, na.rm = FALSE), between(-1:1))
  expect_equal(fbetween(x), between(x, na.rm = TRUE))
  expect_equal(fbetween(x, na.rm = FALSE), between(x))
  expect_equal(fbetween(xNA, na.rm = FALSE), between(xNA))
  expect_equal(fbetween(xNA), between(xNA, na.rm = TRUE))
  expect_equal(qM(fbetween(mtcars)), fbetween(m))
  expect_equal(fbetween(m), dapply(m, between, na.rm = TRUE))
  expect_equal(fbetween(m, na.rm = FALSE), dapply(m, between))
  expect_equal(fbetween(mNA, na.rm = FALSE), dapply(mNA, between))
  expect_equal(fbetween(mNA), dapply(mNA, between, na.rm = TRUE))
  expect_equal(fbetween(mtcars), dapply(mtcars, between, na.rm = TRUE))
  expect_equal(fbetween(mtcars, na.rm = FALSE), dapply(mtcars, between))
  expect_equal(fbetween(mtcNA, na.rm = FALSE), dapply(mtcNA, between))
  expect_equal(fbetween(mtcNA), dapply(mtcNA, between, na.rm = TRUE))
  expect_equal(fbetween(x, f), BY(x, f, between, na.rm = TRUE, use.g.names = FALSE))
  expect_equal(fbetween(x, f, na.rm = FALSE), BY(x, f, between, use.g.names = FALSE))
  expect_equal(fbetween(xNA, f, na.rm = FALSE), BY(xNA, f, between, use.g.names = FALSE))
  expect_equal(fbetween(xNA, f), BY(xNA, f, between, na.rm = TRUE, use.g.names = FALSE))
  expect_equal(fbetween(m, g), BY(m, g, between, na.rm = TRUE, use.g.names = FALSE))
  expect_equal(fbetween(m, g, na.rm = FALSE), BY(m, g, between, use.g.names = FALSE))
  expect_equal(fbetween(mNA, g, na.rm = FALSE), BY(mNA, g, between, use.g.names = FALSE))
  expect_equal(fbetween(mNA, g), BY(mNA, g, between, na.rm = TRUE, use.g.names = FALSE))
  expect_equal(fbetween(mtcars, g), BY(mtcars, g, between, na.rm = TRUE, use.g.names = FALSE))
  expect_equal(fbetween(mtcars, g, na.rm = FALSE), BY(mtcars, g, between, use.g.names = FALSE))
  expect_equal(fbetween(mtcNA, g, na.rm = FALSE), BY(mtcNA, g, between, use.g.names = FALSE))
  expect_equal(fbetween(mtcNA, g), BY(mtcNA, g, between, na.rm = TRUE, use.g.names = FALSE))
})

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

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

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

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

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

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

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

  expect_equal(fbetween(NA, w = NA), NA_real_)
  expect_equal(fbetween(NaN, w = NA), NaN)
  expect_equal(fbetween(Inf, w = NA), NA_real_)
  expect_equal(fbetween(c(Inf,Inf), w = c(NA,2)), c(Inf,Inf))
  expect_equal(fbetween(-Inf, w = NA), NA_real_)
  expect_equal(fbetween(c(-Inf,-Inf), w = c(NA,2)), c(-Inf,-Inf))
  expect_equal(fbetween(TRUE, w = NA), NA_real_)
  expect_equal(fbetween(FALSE, w = NA), NA_real_)
  expect_equal(fbetween(NA, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fbetween(NaN, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fbetween(Inf, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fbetween(c(Inf,Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fbetween(-Inf, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fbetween(c(-Inf,-Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fbetween(TRUE, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fbetween(FALSE, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fbetween(c(1,NA), w = c(NA,2)), c(NA_real_,NA_real_))
  expect_equal(fbetween(c(1,NaN), w = c(NA,2)), c(NA_real_,NA_real_))
  expect_equal(fbetween(c(1,Inf), w = c(NA,2)), c(Inf,Inf))
  expect_equal(fbetween(c(1,-Inf), w = c(NA,2)), c(-Inf,-Inf))
  expect_equal(fbetween(c(1,Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fbetween(c(1,-Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fbetween(c(NA,-Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fbetween(c(FALSE,FALSE), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fbetween(c(1,NA), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fbetween(1:3, w = c(NA,Inf,3)), c(NA_real_,NA_real_,NA_real_))
  expect_equal(fbetween(1:3, w = c(NA,-Inf,3)), c(NA_real_,NA_real_,NA_real_))
  expect_equal(fbetween(1:3, w = c(NA,Inf,3), na.rm = FALSE), c(NA_real_,NA_real_,NA_real_))
  expect_equal(fbetween(1:3, w = c(NA,-Inf,3), na.rm = FALSE), c(NA_real_,NA_real_,NA_real_))
})

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

# B

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

test_that("B.data.frame method is foolproof", {
  expect_visible(B(wlddev))
  expect_visible(B(wlddev, w = wlddev$year))
  expect_visible(B(wlddev, w = ~year))
  expect_visible(B(wlddev, wlddev$iso3c))
  expect_visible(B(wlddev, ~iso3c))
  expect_visible(B(wlddev, ~iso3c + region))
  expect_visible(B(wlddev, wlddev$iso3c, wlddev$year))
  expect_visible(B(wlddev, ~iso3c, ~year))
  expect_visible(B(wlddev, cols = 9:12))
  expect_visible(B(wlddev, w = wlddev$year, cols = 9:12))
  expect_visible(B(wlddev, w = ~year, cols = 9:12))
  expect_visible(B(wlddev, wlddev$iso3c, cols = 9:12))
  expect_visible(B(wlddev, ~iso3c, cols = 9:12))
  expect_visible(B(wlddev, wlddev$iso3c, wlddev$year, cols = 9:12))
  expect_visible(B(wlddev, ~iso3c, ~year, cols = 9:12))
  expect_visible(B(wlddev, cols = c("PCGDP","LIFEEX")))
  expect_visible(B(wlddev, w = wlddev$year, cols = c("PCGDP","LIFEEX")))
  expect_visible(B(wlddev, w = ~year, cols = c("PCGDP","LIFEEX")))
  expect_visible(B(wlddev, wlddev$iso3c, cols = c("PCGDP","LIFEEX")))
  expect_visible(B(wlddev, ~iso3c, cols = c("PCGDP","LIFEEX")))
  expect_visible(B(wlddev, wlddev$iso3c, wlddev$year, cols = c("PCGDP","LIFEEX")))
  expect_visible(B(wlddev, ~iso3c, ~year, cols = c("PCGDP","LIFEEX")))

  expect_error(B(wlddev, cols = NULL))
  expect_error(B(wlddev, w = wlddev$year, cols = NULL))
  expect_error(B(wlddev, w = ~year, cols = NULL))
  expect_error(B(wlddev, wlddev$iso3c, cols = NULL))
  expect_error(B(wlddev, ~iso3c, cols = NULL))
  expect_error(B(wlddev, wlddev$iso3c, wlddev$year, cols = NULL))
  expect_error(B(wlddev, ~iso3c, ~year, cols = NULL))
  expect_error(B(wlddev, cols = 9:14))
  expect_error(B(wlddev, w = wlddev$year, cols = 9:14))
  expect_error(B(wlddev, w = ~year, cols = 9:14))
  expect_error(B(wlddev, wlddev$iso3c, cols = 9:14))
  expect_error(B(wlddev, ~iso3c, cols = 9:14))
  expect_error(B(wlddev, wlddev$iso3c, wlddev$year, cols = 9:14))
  expect_error(B(wlddev, ~iso3c, ~year, cols = 9:14))
  expect_error(B(wlddev, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(B(wlddev, w = wlddev$year, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(B(wlddev, w = ~year, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(B(wlddev, wlddev$iso3c, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(B(wlddev, ~iso3c, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(B(wlddev, wlddev$iso3c, wlddev$year, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(B(wlddev, ~iso3c, ~year, cols = c("PCGDP","LIFEEX","bla")))

  expect_error(B(wlddev, w = mtcars))
  expect_error(B(wlddev, w = 4))
  expect_error(B(wlddev, w = "year"))
  expect_error(B(wlddev, w = ~year2))
  # suppressWarnings(expect_error(B(wlddev, w = ~year + region)))
  expect_error(B(wlddev, mtcars))
  expect_error(B(wlddev, 2))
  expect_error(B(wlddev, "iso3c"))
  expect_error(B(wlddev, ~iso3c2))
  expect_error(B(wlddev, ~iso3c + bla))
  expect_error(B(wlddev, mtcars$mpg, mtcars$cyl))
  expect_error(B(wlddev, 2, 4))
  expect_error(B(wlddev, ~iso3c2, ~year2))
  expect_error(B(wlddev, cols = ~bla))
  expect_error(B(wlddev, w = ~bla, cols = 9:12))
  expect_error(B(wlddev, w = 4, cols = 9:12))
  expect_error(B(wlddev, w = "year", cols = 9:12))
  expect_error(B(wlddev, w = ~yewar, cols = 9:12))
  expect_error(B(wlddev, mtcars$mpg, cols = 9:12))
  expect_error(B(wlddev, ~iso3c + ss, cols = 9:12))
  expect_error(B(wlddev, 2, cols = 9:12))
  expect_error(B(wlddev, "iso3c", cols = 9:12))
  expect_error(B(wlddev, wlddev$iso3c, ~year + bla, cols = 9:12))
  expect_error(B(wlddev, ~iso3c3, ~year, cols = 9:12))
  expect_error(B(wlddev, cols = c("PC3GDP","LIFEEX")))
})


# fwithin

test_that("fwithin performs like within", {
  expect_equal(fwithin(NA), as.double(within(NA)))
  expect_equal(fwithin(NA, na.rm = FALSE), as.double(within(NA)))
  expect_equal(fwithin(1), within(1, na.rm = TRUE))
  expect_equal(fwithin(1:3), within(1:3, na.rm = TRUE))
  expect_equal(fwithin(-1:1), within(-1:1, na.rm = TRUE))
  expect_equal(fwithin(1, na.rm = FALSE), within(1))
  expect_equal(fwithin(1:3, na.rm = FALSE), within(1:3))
  expect_equal(fwithin(-1:1, na.rm = FALSE), within(-1:1))
  expect_equal(fwithin(x), within(x, na.rm = TRUE))
  expect_equal(fwithin(x, na.rm = FALSE), within(x))
  expect_equal(fwithin(xNA, na.rm = FALSE), within(xNA))
  expect_equal(fwithin(xNA), within(xNA, na.rm = TRUE))
  expect_equal(qM(fwithin(mtcars)), fwithin(m))
  expect_equal(fwithin(m), dapply(m, within, na.rm = TRUE))
  expect_equal(fwithin(m, na.rm = FALSE), dapply(m, within))
  expect_equal(fwithin(mNA, na.rm = FALSE), dapply(mNA, within))
  expect_equal(fwithin(mNA), dapply(mNA, within, na.rm = TRUE))
  expect_equal(fwithin(mtcars), dapply(mtcars, within, na.rm = TRUE))
  expect_equal(fwithin(mtcars, na.rm = FALSE), dapply(mtcars, within))
  expect_equal(fwithin(mtcNA, na.rm = FALSE), dapply(mtcNA, within))
  expect_equal(fwithin(mtcNA), dapply(mtcNA, within, na.rm = TRUE))
  expect_equal(fwithin(x, f), BY(x, f, within, na.rm = TRUE))
  expect_equal(fwithin(x, f, na.rm = FALSE), BY(x, f, within))
  expect_equal(fwithin(xNA, f, na.rm = FALSE), BY(xNA, f, within))
  expect_equal(fwithin(xNA, f), BY(xNA, f, within, na.rm = TRUE))
  expect_equal(fwithin(m, g), BY(m, g, within, na.rm = TRUE))
  expect_equal(fwithin(m, g, na.rm = FALSE), BY(m, g, within))
  expect_equal(fwithin(mNA, g, na.rm = FALSE), BY(mNA, g, within))
  expect_equal(fwithin(mNA, g), BY(mNA, g, within, na.rm = TRUE))
  expect_equal(fwithin(mtcars, g), BY(mtcars, g, within, na.rm = TRUE))
  expect_equal(fwithin(mtcars, g, na.rm = FALSE), BY(mtcars, g, within))
  expect_equal(fwithin(mtcNA, g, na.rm = FALSE), BY(mtcNA, g, within))
  expect_equal(fwithin(mtcNA, g), BY(mtcNA, g, within, na.rm = TRUE))
})

test_that("fwithin with custom mean performs like within (defined above)", {
  expect_equal(fwithin(x, mean = 4.8456), within(x, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(x, na.rm = FALSE, mean = 4.8456), within(x, mean = 4.8456))
  expect_equal(fwithin(xNA, na.rm = FALSE, mean = 4.8456), within(xNA, mean = 4.8456))
  expect_equal(fwithin(xNA, mean = 4.8456), within(xNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(qM(fwithin(mtcars, mean = 4.8456)), fwithin(m, mean = 4.8456))
  expect_equal(fwithin(m, mean = 4.8456), dapply(m, within, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(m, na.rm = FALSE, mean = 4.8456), dapply(m, within, mean = 4.8456))
  expect_equal(fwithin(mNA, na.rm = FALSE, mean = 4.8456), dapply(mNA, within, mean = 4.8456))
  expect_equal(fwithin(mNA, mean = 4.8456), dapply(mNA, within, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(x, f, mean = 4.8456), BY(x, f, within, na.rm = TRUE, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(x, f, na.rm = FALSE, mean = 4.8456), BY(x, f, within, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(xNA, f, na.rm = FALSE, mean = 4.8456), BY(xNA, f, within, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(xNA, f, mean = 4.8456), BY(xNA, f, within, na.rm = TRUE, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(m, g, mean = 4.8456), BY(m, g, within, na.rm = TRUE, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(m, g, na.rm = FALSE, mean = 4.8456), BY(m, g, within, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(mNA, g, na.rm = FALSE, mean = 4.8456), BY(mNA, g, within, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(mNA, g, mean = 4.8456), BY(mNA, g, within, na.rm = TRUE, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(mtcars, g, mean = 4.8456), BY(mtcars, g, within, na.rm = TRUE, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(mtcars, g, na.rm = FALSE, mean = 4.8456), BY(mtcars, g, within, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(mtcNA, g, na.rm = FALSE, mean = 4.8456), BY(mtcNA, g, within, use.g.names = FALSE, mean = 4.8456))
  expect_equal(fwithin(mtcNA, g, mean = 4.8456), BY(mtcNA, g, within, na.rm = TRUE, use.g.names = FALSE, mean = 4.8456))
})

test_that("Centering on overall mean performs as indended", {
  expect_equal(fwithin(x, f, mean = "overall.mean"), BY(x, f, within, na.rm = TRUE, use.g.names = FALSE) + ave(x))
  expect_equal(fwithin(x, f, na.rm = FALSE, mean = "overall.mean"), BY(x, f, within, use.g.names = FALSE) + ave(x))
  # expect_equal(fwithin(xNA, f, na.rm = FALSE, mean = "overall.mean"), BY(xNA, f, within, use.g.names = FALSE) + B(xNA)) # Not the same !!
  expect_equal(fwithin(xNA, f, mean = "overall.mean"), BY(xNA, f, within, na.rm = TRUE, use.g.names = FALSE) + B(xNA))
  expect_equal(fwithin(m, g, mean = "overall.mean"), BY(m, g, within, na.rm = TRUE, use.g.names = FALSE) + B(m))
  expect_equal(fwithin(m, g, na.rm = FALSE, mean = "overall.mean"), BY(m, g, within, use.g.names = FALSE) + B(m))
  # expect_equal(fwithin(mNA, g, na.rm = FALSE, mean = "overall.mean"), BY(mNA, g, within, use.g.names = FALSE) + B(mNA))
  expect_equal(fwithin(mNA, g, mean = "overall.mean"), BY(mNA, g, within, na.rm = TRUE, use.g.names = FALSE) + B(mNA))
  expect_equal(fwithin(mtcars, g, mean = "overall.mean"), BY(mtcars, g, within, na.rm = TRUE, use.g.names = FALSE) + B(mtcars))
  expect_equal(fwithin(mtcars, g, na.rm = FALSE, mean = "overall.mean"), BY(mtcars, g, within, use.g.names = FALSE) + B(mtcars))
  # expect_equal(fwithin(mtcNA, g, na.rm = FALSE, mean = "overall.mean"), BY(mtcNA, g, within, use.g.names = FALSE) + B(mtcNA))
  expect_equal(fwithin(mtcNA, g, mean = "overall.mean"), BY(mtcNA, g, within, na.rm = TRUE, use.g.names = FALSE)+ B(mtcNA))
})

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

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

test_that("fwithin with custom mean and weights performs like wwithin (defined above)", {
  # complete weights
  expect_equal(fwithin(x, w = w, mean = 4.8456), wwithin(x, w, mean = 4.8456))
  expect_equal(fwithin(x, w = w, na.rm = FALSE, mean = 4.8456), wwithin(x, w, mean = 4.8456))
  expect_equal(fwithin(xNA, w = w, na.rm = FALSE, mean = 4.8456), wwithin(xNA, w, mean = 4.8456))
  expect_equal(fwithin(xNA, w = w, mean = 4.8456), wwithin(xNA, w, na.rm = TRUE, mean = 4.8456))
  expect_equal(qM(fwithin(mtcars, w = wdat, mean = 4.8456)), fwithin(m, w = wdat, mean = 4.8456))
  expect_equal(fwithin(m, w = wdat, mean = 4.8456), dapply(m, wwithin, wdat, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(m, w = wdat, na.rm = FALSE, mean = 4.8456), dapply(m, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mNA, w = wdat, na.rm = FALSE, mean = 4.8456), dapply(mNA, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mNA, w = wdat, mean = 4.8456), dapply(mNA, wwithin, wdat, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(mtcars, w = wdat, mean = 4.8456), dapply(mtcars, wwithin, wdat, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(mtcars, w = wdat, na.rm = FALSE, mean = 4.8456), dapply(mtcars, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mtcNA, w = wdat, na.rm = FALSE, mean = 4.8456), dapply(mtcNA, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mtcNA, w = wdat, mean = 4.8456), dapply(mtcNA, wwithin, wdat, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(x, f, w, mean = 4.8456), BY(x, f, wwithin, w, mean = 4.8456))
  expect_equal(fwithin(x, f, w, na.rm = FALSE, mean = 4.8456), BY(x, f, wwithin, w, mean = 4.8456))
  expect_equal(fwithin(xNA, f, w, na.rm = FALSE, mean = 4.8456), BY(xNA, f, wwithin, w, mean = 4.8456))
  expect_equal(fwithin(xNA, f, w, mean = 4.8456), BY(xNA, f, wwithin, w, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(m, g, wdat, mean = 4.8456), BY(m, g, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(m, g, wdat, na.rm = FALSE, mean = 4.8456), BY(m, g, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mNA, g, wdat, na.rm = FALSE, mean = 4.8456),  BY(mNA, g, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mNA, g, wdat, mean = 4.8456), BY(mNA, g, wwithin, wdat, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(mtcars, g, wdat, mean = 4.8456), BY(mtcars, g, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mtcars, g, wdat, na.rm = FALSE, mean = 4.8456), BY(mtcars, g, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mtcNA, g, wdat, na.rm = FALSE, mean = 4.8456), BY(mtcNA, g, wwithin, wdat, mean = 4.8456))
  expect_equal(fwithin(mtcNA, g, wdat, mean = 4.8456), BY(mtcNA, g, wwithin, wdat, na.rm = TRUE, mean = 4.8456))
  # missing weights
  expect_equal(fwithin(x, w = wNA, mean = 4.8456), wwithin(x, wNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(x, w = wNA, na.rm = FALSE, mean = 4.8456), wwithin(x, wNA, mean = 4.8456))
  expect_equal(fwithin(xNA, w = wNA, na.rm = FALSE, mean = 4.8456), wwithin(xNA, wNA, mean = 4.8456))
  expect_equal(fwithin(xNA, w = wNA, mean = 4.8456), wwithin(xNA, wNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(qM(fwithin(mtcars, w = wdatNA, mean = 4.8456)), fwithin(m, w = wdatNA, mean = 4.8456))
  expect_equal(fwithin(m, w = wdatNA, mean = 4.8456), dapply(m, wwithin, wdatNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(m, w = wdatNA, na.rm = FALSE, mean = 4.8456), dapply(m, wwithin, wdatNA, mean = 4.8456))
  expect_equal(fwithin(mNA, w = wdatNA, na.rm = FALSE, mean = 4.8456), dapply(mNA, wwithin, wdatNA, mean = 4.8456))
  expect_equal(fwithin(mNA, w = wdatNA, mean = 4.8456), dapply(mNA, wwithin, wdatNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(mtcars, w = wdatNA, mean = 4.8456), dapply(mtcars, wwithin, wdatNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(mtcars, w = wdatNA, na.rm = FALSE, mean = 4.8456), dapply(mtcars, wwithin, wdatNA, mean = 4.8456))
  expect_equal(fwithin(mtcNA, w = wdatNA, na.rm = FALSE, mean = 4.8456), dapply(mtcNA, wwithin, wdatNA, mean = 4.8456))
  expect_equal(fwithin(mtcNA, w = wdatNA, mean = 4.8456), dapply(mtcNA, wwithin, wdatNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(x, f, wNA, mean = 4.8456), BY(x, f, wwithin, wNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(x, f, wNA, na.rm = FALSE, mean = 4.8456), BY(x, f, wwithin, wNA, mean = 4.8456))
  expect_equal(fwithin(xNA, f, wNA, na.rm = FALSE, mean = 4.8456), BY(xNA, f, wwithin, wNA, mean = 4.8456))
  expect_equal(fwithin(xNA, f, wNA, mean = 4.8456), BY(xNA, f, wwithin, wNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(m, g, wdatNA, mean = 4.8456), BY(m, g, wwithin, wdatNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(m, g, wdatNA, na.rm = FALSE, mean = 4.8456), BY(m, g, wwithin, wdatNA, mean = 4.8456))
  expect_equal(fwithin(mNA, g, wdatNA, na.rm = FALSE, mean = 4.8456),  BY(mNA, g, wwithin, wdatNA, mean = 4.8456))
  expect_equal(fwithin(mNA, g, wdatNA, mean = 4.8456), BY(mNA, g, wwithin, wdatNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(mtcars, g, wdatNA, mean = 4.8456), BY(mtcars, g, wwithin, wdatNA, na.rm = TRUE, mean = 4.8456))
  expect_equal(fwithin(mtcars, g, wdatNA, na.rm = FALSE, mean = 4.8456), BY(mtcars, g, wwithin, wdatNA, mean = 4.8456))
  expect_equal(fwithin(mtcNA, g, wdatNA, na.rm = FALSE, mean = 4.8456), BY(mtcNA, g, wwithin, wdatNA, mean = 4.8456))
  expect_equal(fwithin(mtcNA, g, wdatNA, mean = 4.8456), BY(mtcNA, g, wwithin, wdatNA, na.rm = TRUE, mean = 4.8456))
})

test_that("Weighted centering on overall weighted mean performs as indended", {
  # complete weights
  expect_equal(fwithin(x, f, w, mean = "overall.mean"), as.numeric(mapply(wwithin, split(x, f), split(w, f), na.rm = TRUE)) + B(x, w = w))
  expect_equal(fwithin(x, f, w, na.rm = FALSE, mean = "overall.mean"), as.numeric(mapply(wwithin, split(x, f), split(w, f))) + B(x, w = w))
  # expect_equal(fwithin(xNA, f, w, na.rm = FALSE, mean = "overall.mean"), as.numeric(mapply(wwithin, split(xNA, f), split(w, f))) + B(xNA, w = w)) # Not the same !!
  expect_equal(fwithin(xNA, f, w, mean = "overall.mean"), as.numeric(mapply(wwithin, split(xNA, f), split(w, f), na.rm = TRUE)) + B(xNA, w = w))
})
# Do more than this to test the rest ...

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

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

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

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

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

  expect_equal(fwithin(NA, w = NA), NA_real_)
  expect_equal(fwithin(NaN, w = NA), NaN)
  expect_equal(fwithin(Inf, w = NA), NaN)
  expect_equal(fwithin(c(Inf,Inf), w = c(NA,2)), c(NaN,NaN))
  expect_equal(fwithin(-Inf, w = NA), NA_real_)
  expect_equal(fwithin(c(-Inf,-Inf), w = c(NA,2)), c(NaN,NaN))
  expect_equal(fwithin(TRUE, w = NA), NA_real_)
  expect_equal(fwithin(FALSE, w = NA), NA_real_)
  expect_equal(fwithin(NA, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fwithin(NaN, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fwithin(Inf, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fwithin(c(Inf,Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fwithin(-Inf, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fwithin(c(-Inf,-Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fwithin(TRUE, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fwithin(FALSE, w = NA, na.rm = FALSE), NA_real_)
  expect_equal(fwithin(c(1,NA), w = c(NA,2)), c(NA_real_,NA_real_))
  expect_equal(fwithin(c(1,NaN), w = c(NA,2)), c(NA_real_,NA_real_))
  expect_equal(fwithin(c(1,Inf), w = c(NA,2)), c(-Inf,NaN))
  expect_equal(fwithin(c(1,-Inf), w = c(NA,2)), c(Inf,NaN))
  expect_equal(fwithin(c(1,Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fwithin(c(1,-Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fwithin(c(NA,-Inf), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fwithin(c(FALSE,FALSE), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fwithin(c(1,NA), w = c(NA,2), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fwithin(1:3, w = c(NA,Inf,3)), c(NaN,NaN,NaN))
  expect_equal(fwithin(1:3, w = c(NA,-Inf,3)), c(NaN,NaN,NaN))
  expect_equal(fwithin(1:3, w = c(NA,Inf,3), na.rm = FALSE), c(NA_real_,NA_real_,NA_real_))
  expect_equal(fwithin(1:3, w = c(NA,-Inf,3), na.rm = FALSE), c(NA_real_,NA_real_,NA_real_))
})

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

test_that("fwithin shoots errors for wrong input to mean", {
  expect_error(fwithin(x, mean = FALSE))
  expect_error(fwithin(m, mean = FALSE))
  expect_error(fwithin(mtcars, mean = FALSE))
  expect_error(fwithin(x, mean = "overall.mean"))
  expect_error(fwithin(m, mean = "overall.mean"))
  expect_error(fwithin(mtcars, mean = "overall.mean"))
  expect_error(fwithin(m, mean = fmean(m)))
  expect_error(fwithin(mtcars, mean = fmean(mtcars)))
})


# W

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

test_that("W.data.frame method is foolproof", {
  expect_visible(W(wlddev))
  expect_visible(W(wlddev, w = wlddev$year))
  expect_visible(W(wlddev, w = ~year))
  expect_visible(W(wlddev, wlddev$iso3c))
  expect_visible(W(wlddev, ~iso3c))
  expect_visible(W(wlddev, ~iso3c + region))
  expect_visible(W(wlddev, wlddev$iso3c, wlddev$year))
  expect_visible(W(wlddev, ~iso3c, ~year))
  expect_visible(W(wlddev, cols = 9:12))
  expect_visible(W(wlddev, w = wlddev$year, cols = 9:12))
  expect_visible(W(wlddev, w = ~year, cols = 9:12))
  expect_visible(W(wlddev, wlddev$iso3c, cols = 9:12))
  expect_visible(W(wlddev, ~iso3c, cols = 9:12))
  expect_visible(W(wlddev, wlddev$iso3c, wlddev$year, cols = 9:12))
  expect_visible(W(wlddev, ~iso3c, ~year, cols = 9:12))
  expect_visible(W(wlddev, cols = c("PCGDP","LIFEEX")))
  expect_visible(W(wlddev, w = wlddev$year, cols = c("PCGDP","LIFEEX")))
  expect_visible(W(wlddev, w = ~year, cols = c("PCGDP","LIFEEX")))
  expect_visible(W(wlddev, wlddev$iso3c, cols = c("PCGDP","LIFEEX")))
  expect_visible(W(wlddev, ~iso3c, cols = c("PCGDP","LIFEEX")))
  expect_visible(W(wlddev, wlddev$iso3c, wlddev$year, cols = c("PCGDP","LIFEEX")))
  expect_visible(W(wlddev, ~iso3c, ~year, cols = c("PCGDP","LIFEEX")))

  expect_error(W(wlddev, cols = NULL))
  expect_error(W(wlddev, w = wlddev$year, cols = NULL))
  expect_error(W(wlddev, w = ~year, cols = NULL))
  expect_error(W(wlddev, wlddev$iso3c, cols = NULL))
  expect_error(W(wlddev, ~iso3c, cols = NULL))
  expect_error(W(wlddev, wlddev$iso3c, wlddev$year, cols = NULL))
  expect_error(W(wlddev, ~iso3c, ~year, cols = NULL))
  expect_error(W(wlddev, cols = 9:14))
  expect_error(W(wlddev, w = wlddev$year, cols = 9:14))
  expect_error(W(wlddev, w = ~year, cols = 9:14))
  expect_error(W(wlddev, wlddev$iso3c, cols = 9:14))
  expect_error(W(wlddev, ~iso3c, cols = 9:14))
  expect_error(W(wlddev, wlddev$iso3c, wlddev$year, cols = 9:14))
  expect_error(W(wlddev, ~iso3c, ~year, cols = 9:14))
  expect_error(W(wlddev, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(W(wlddev, w = wlddev$year, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(W(wlddev, w = ~year, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(W(wlddev, wlddev$iso3c, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(W(wlddev, ~iso3c, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(W(wlddev, wlddev$iso3c, wlddev$year, cols = c("PCGDP","LIFEEX","bla")))
  expect_error(W(wlddev, ~iso3c, ~year, cols = c("PCGDP","LIFEEX","bla")))

  expect_error(W(wlddev, w = mtcars))
  expect_error(W(wlddev, w = 4))
  expect_error(W(wlddev, w = "year"))
  expect_error(W(wlddev, w = ~year2))
  # suppressWarnings(expect_error(W(wlddev, w = ~year + region)))
  expect_error(W(wlddev, mtcars))
  expect_error(W(wlddev, 2))
  expect_error(W(wlddev, "iso3c"))
  expect_error(W(wlddev, ~iso3c2))
  expect_error(W(wlddev, ~iso3c + bla))
  expect_error(W(wlddev, mtcars$mpg, mtcars$cyl))
  expect_error(W(wlddev, 2, 4))
  expect_error(W(wlddev, ~iso3c2, ~year2))
  expect_error(W(wlddev, cols = ~bla))
  expect_error(W(wlddev, w = ~bla, cols = 9:12))
  expect_error(W(wlddev, w = 4, cols = 9:12))
  expect_error(W(wlddev, w = "year", cols = 9:12))
  expect_error(W(wlddev, w = ~yewar, cols = 9:12))
  expect_error(W(wlddev, mtcars$mpg, cols = 9:12))
  expect_error(W(wlddev, ~iso3c + ss, cols = 9:12))
  expect_error(W(wlddev, 2, cols = 9:12))
  expect_error(W(wlddev, "iso3c", cols = 9:12))
  expect_error(W(wlddev, wlddev$iso3c, ~year + bla, cols = 9:12))
  expect_error(W(wlddev, ~iso3c3, ~year, cols = 9:12))
  expect_error(W(wlddev, cols = c("PC3GDP","LIFEEX")))
})