context("fvar and fsd")

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

bvar <- stats::var
bsd <- stats::sd
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"

na20 <- function(x) {
  x[is.na(x)] <- 0
  x
}

# This is correct, including Bessels correction.
wvar <- function(x, w, na.rm = FALSE) {
  if(na.rm) {
    cc <- complete.cases(x, w)
    x <- x[cc]
    # if(length(x) < 2L) return(NA_real_)
    w <- w[cc]
  } # else if(length(x) < 2L) return(if(is.na(x)) NA_real_ else 0)
  bsum(w*(x-weighted.mean(x,w))^2)/(bsum(w)-1)
}


# fvar using Welford's Algoritm (default)

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

test_that("fvar with weights performs as intended (unbiased)", {
  expect_equal(fvar(c(2,2,4,5,5,5)), fvar(c(2,4,5), w = c(2,1,3)))
  expect_equal(fvar(c(2,2,4,5,5,5), na.rm = FALSE), fvar(c(2,4,5), w = c(2,1,3), na.rm = FALSE))
  expect_equal(fvar(c(2.456,2.456,4.123,5.009,5.009,5.009)), fvar(c(2.456,4.123,5.009), w = c(2,1,3)))
  expect_equal(fvar(c(2.456,2.456,4.123,5.009,5.009,5.009), na.rm = FALSE), fvar(c(2.456,4.123,5.009), w = c(2,1,3), na.rm = FALSE))
  expect_equal(fvar(c(2,2,NA,5,5,5)), fvar(c(2,NA,5), w = c(2,1,3)))
  expect_equal(fvar(c(2,2,NA,5,5,5), na.rm = FALSE), fvar(c(2,NA,5), w = c(2,1,3), na.rm = FALSE))
  expect_equal(fvar(c(2,2,NA,5,5,5)), fvar(c(2,4,5), w = c(2,NA,3)))
  expect_equal(fvar(c(2,2,NA,5,5,5), na.rm = FALSE), fvar(c(2,4,5), w = c(2,NA,3), na.rm = FALSE))
  expect_equal(fvar(c(NA,NA,4.123,5.009,5.009,5.009)), fvar(c(NA,4.123,5.009), w = c(2,1,3)))
  expect_equal(fvar(c(NA,NA,4.123,5.009,5.009,5.009), na.rm = FALSE), fvar(c(NA,4.123,5.009), w = c(2,1,3), na.rm = FALSE))
  expect_equal(fvar(c(NA,NA,4.123,5.009,5.009,5.009)), fvar(c(2.456,4.123,5.009), w = c(NA,1,3)))
  expect_equal(fvar(c(NA,NA,4.123,5.009,5.009,5.009), na.rm = FALSE), fvar(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(fvar(v, f), fvar(vs, fs, w))
  expect_equal(fvar(v, f, na.rm = FALSE), fvar(vs, fs, w, na.rm = FALSE))
  expect_equal(fvar(v2, f), fvar(v2s, fs, w))
  expect_equal(fvar(v2, f, na.rm = FALSE), fvar(v2s, fs, w, na.rm = FALSE))
  v[c(3,9)] <- NA; vs[c(2,5)] <- NA
  expect_equal(fvar(v, f), fvar(vs, fs, w))
  expect_equal(fvar(v, f, na.rm = FALSE), fvar(vs, fs, w, na.rm = FALSE))
  vs[c(2,5)] <- 4; w[c(2,5)] <- NA
  expect_equal(fvar(v, f), fvar(vs, fs, w))
  expect_equal(fvar(v, f, na.rm = FALSE), fvar(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(fvar(v2, f), fvar(v2s, fs, w))
  expect_equal(fvar(v2, f, na.rm = FALSE), fvar(v2s, fs, w, na.rm = FALSE))
  v2s[c(1,4)] <- 2.456; w[c(1,4)] <- NA
  expect_equal(fvar(v2, f), fvar(v2s, fs, w))
  expect_equal(fvar(v2, f, na.rm = FALSE), fvar(v2s, fs, w, na.rm = FALSE))
})

test_that("fvar performs like fvar with unit weights", {
  expect_equal(fvar(NA), fvar(NA, w = 1))
  expect_equal(fvar(NA, na.rm = FALSE), fvar(NA, w = 1, na.rm = FALSE))
  expect_equal(fvar(1), fvar(1, w = 1))
  expect_equal(fvar(1:3), fvar(1:3, w = rep(1,3)))
  expect_equal(fvar(-1:1), fvar(-1:1, w = rep(1,3)))
  expect_equal(fvar(1, na.rm = FALSE), fvar(1, w = 1, na.rm = FALSE))
  expect_equal(fvar(1:3, na.rm = FALSE), fvar(1:3, w = rep(1, 3), na.rm = FALSE))
  expect_equal(fvar(-1:1, na.rm = FALSE), fvar(-1:1, w = rep(1, 3), na.rm = FALSE))
  expect_equal(fvar(x), fvar(x, w = rep(1,100)))
  expect_equal(fvar(x, na.rm = FALSE), fvar(x, w = rep(1, 100), na.rm = FALSE))
  expect_equal(fvar(xNA, na.rm = FALSE), fvar(xNA, w = rep(1, 100), na.rm = FALSE))
  expect_equal(fvar(xNA), fvar(xNA, w = rep(1, 100)))
  expect_equal(fvar(m), fvar(m, w = rep(1, 32)))
  expect_equal(fvar(m, na.rm = FALSE), fvar(m, w = rep(1, 32), na.rm = FALSE))
  expect_equal(fvar(mNA, na.rm = FALSE), fvar(mNA, w = rep(1, 32), na.rm = FALSE))
  expect_equal(fvar(mNA), fvar(mNA, w = rep(1, 32)))
  expect_equal(fvar(mtcars), fvar(mtcars, w = rep(1, 32)))
  expect_equal(fvar(mtcars, na.rm = FALSE), fvar(mtcars, w = rep(1, 32), na.rm = FALSE))
  expect_equal(fvar(mtcNA, na.rm = FALSE), fvar(mtcNA, w = rep(1, 32), na.rm = FALSE))
  expect_equal(fvar(mtcNA), fvar(mtcNA, w = rep(1, 32)))
  expect_equal(fvar(x, f), fvar(x, f, rep(1,100)))
  expect_equal(fvar(x, f, na.rm = FALSE), fvar(x, f, rep(1,100), na.rm = FALSE))
  expect_equal(fvar(xNA, f, na.rm = FALSE), fvar(xNA, f, rep(1,100), na.rm = FALSE))
  expect_equal(fvar(xNA, f), fvar(xNA, f, rep(1,100)))
  expect_equal(fvar(m, g), fvar(m, g, rep(1,32)))
  expect_equal(fvar(m, g, na.rm = FALSE), fvar(m, g, rep(1,32), na.rm = FALSE))
  expect_equal(fvar(mNA, g, na.rm = FALSE), fvar(mNA, g, rep(1,32), na.rm = FALSE))
  expect_equal(fvar(mNA, g), fvar(mNA, g, rep(1,32)))
  expect_equal(fvar(mtcars, g), fvar(mtcars, g, rep(1,32)))
  expect_equal(fvar(mtcars, g, na.rm = FALSE), fvar(mtcars, g, rep(1,32), na.rm = FALSE))
  expect_equal(fvar(mtcNA, g, na.rm = FALSE), fvar(mtcNA, g, rep(1,32), na.rm = FALSE))
  expect_equal(fvar(mtcNA, g), fvar(mtcNA, g, rep(1,32)))
})

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

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

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

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

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

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

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


# Repeating all tests for the other algorithm


test_that("fvar with direct algorithm performs like base::var", {
  expect_equal(fvar(NA, stable.algo = FALSE), bvar(NA))
  expect_equal(fvar(NA, na.rm = FALSE, stable.algo = FALSE), bvar(NA))
  expect_equal(fvar(1, stable.algo = FALSE), bvar(1, na.rm = TRUE))
  expect_equal(fvar(1:3, stable.algo = FALSE), bvar(1:3, na.rm = TRUE))
  expect_equal(fvar(-1:1, stable.algo = FALSE), bvar(-1:1, na.rm = TRUE))
  expect_equal(fvar(1, na.rm = FALSE, stable.algo = FALSE), bvar(1))
  expect_equal(fvar(1:3, na.rm = FALSE, stable.algo = FALSE), bvar(1:3))
  expect_equal(fvar(-1:1, na.rm = FALSE, stable.algo = FALSE), bvar(-1:1))
  expect_equal(fvar(x, stable.algo = FALSE), bvar(x, na.rm = TRUE))
  expect_equal(fvar(x, na.rm = FALSE, stable.algo = FALSE), bvar(x))
  expect_equal(fvar(xNA, na.rm = FALSE, stable.algo = FALSE), bvar(xNA))
  expect_equal(fvar(xNA, stable.algo = FALSE), bvar(xNA, na.rm = TRUE))
  expect_equal(fvar(mtcars, stable.algo = FALSE), fvar(m))
  expect_equal(fvar(m, stable.algo = FALSE), dapply(m, bvar, na.rm = TRUE))
  expect_equal(fvar(m, na.rm = FALSE, stable.algo = FALSE), dapply(m, bvar))
  expect_equal(fvar(mNA, na.rm = FALSE, stable.algo = FALSE), dapply(mNA, bvar))
  expect_equal(fvar(mNA, stable.algo = FALSE), dapply(mNA, bvar, na.rm = TRUE))
  expect_equal(fvar(mtcars, stable.algo = FALSE), dapply(mtcars, bvar, na.rm = TRUE))
  expect_equal(fvar(mtcars, na.rm = FALSE, stable.algo = FALSE), dapply(mtcars, bvar))
  expect_equal(fvar(mtcNA, na.rm = FALSE, stable.algo = FALSE), dapply(mtcNA, bvar))
  expect_equal(fvar(mtcNA, stable.algo = FALSE), dapply(mtcNA, bvar, na.rm = TRUE))
  expect_equal(fvar(x, f, stable.algo = FALSE), BY(x, f, bvar, na.rm = TRUE))
  expect_equal(fvar(x, f, na.rm = FALSE, stable.algo = FALSE), BY(x, f, bvar))
  expect_equal(fvar(xNA, f, na.rm = FALSE, stable.algo = FALSE), BY(xNA, f, bvar))
  expect_equal(fvar(xNA, f, stable.algo = FALSE), BY(xNA, f, bvar, na.rm = TRUE)) # failed?
  # expect_equal(fvar(m, g, stable.algo = FALSE), BY(m, g, bvar, na.rm = TRUE)) # failed on arch i386 CMD check and patched-solaris-x86
  # expect_equal(fvar(m, g, na.rm = FALSE, stable.algo = FALSE), BY(m, g, bvar)) # failed on arch i386 CMD check and patched-solaris-x86
  # expect_equal(fvar(mNA, g, na.rm = FALSE, stable.algo = FALSE), BY(mNA, g, bvar)) # failed on arch i386 CMD check and patched-solaris-x86
  # expect_equal(fvar(mNA, g, stable.algo = FALSE), BY(mNA, g, bvar, na.rm = TRUE)) # failed on arch i386 CMD check and patched-solaris-x86
  # expect_equal(fvar(mtcars, g, stable.algo = FALSE), BY(mtcars, g, bvar, na.rm = TRUE)) # failed on arch i386 CMD check and patched-solaris-x86
  # expect_equal(fvar(mtcars, g, na.rm = FALSE, stable.algo = FALSE), BY(mtcars, g, bvar)) # failed on arch i386 CMD check and patched-solaris-x86
  # expect_equal(fvar(mtcNA, g, na.rm = FALSE, stable.algo = FALSE), BY(mtcNA, g, bvar)) # failed on arch i386 CMD check and patched-solaris-x86
  # expect_equal(fvar(mtcNA, g, stable.algo = FALSE), BY(mtcNA, g, bvar, na.rm = TRUE)) # failed on arch i386 CMD check and patched-solaris-x86
})

test_that("fvar with with direct algorithm and weights performs as intended (unbiased)", {
  expect_equal(fvar(c(2,2,4,5,5,5), stable.algo = FALSE), fvar(c(2,4,5), w = c(2,1,3), stable.algo = FALSE))
  expect_equal(fvar(c(2,2,4,5,5,5), na.rm = FALSE, stable.algo = FALSE), fvar(c(2,4,5), w = c(2,1,3), na.rm = FALSE), stable.algo = FALSE)
  expect_equal(fvar(c(2.456,2.456,4.123,5.009,5.009,5.009), stable.algo = FALSE), fvar(c(2.456,4.123,5.009), w = c(2,1,3), stable.algo = FALSE))
  expect_equal(fvar(c(2.456,2.456,4.123,5.009,5.009,5.009), na.rm = FALSE, stable.algo = FALSE), fvar(c(2.456,4.123,5.009), w = c(2,1,3), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(c(2,2,NA,5,5,5), stable.algo = FALSE), fvar(c(2,NA,5), w = c(2,1,3), stable.algo = FALSE))
  expect_equal(fvar(c(2,2,NA,5,5,5), na.rm = FALSE, stable.algo = FALSE), fvar(c(2,NA,5), w = c(2,1,3), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(c(2,2,NA,5,5,5), stable.algo = FALSE), fvar(c(2,4,5), w = c(2,NA,3), stable.algo = FALSE))
  expect_equal(fvar(c(2,2,NA,5,5,5), na.rm = FALSE, stable.algo = FALSE), fvar(c(2,4,5), w = c(2,NA,3), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(c(NA,NA,4.123,5.009,5.009,5.009), stable.algo = FALSE), fvar(c(NA,4.123,5.009), w = c(2,1,3), stable.algo = FALSE))
  expect_equal(fvar(c(NA,NA,4.123,5.009,5.009,5.009), na.rm = FALSE, stable.algo = FALSE), fvar(c(NA,4.123,5.009), w = c(2,1,3), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(c(NA,NA,4.123,5.009,5.009,5.009), stable.algo = FALSE), fvar(c(2.456,4.123,5.009), w = c(NA,1,3), stable.algo = FALSE))
  expect_equal(fvar(c(NA,NA,4.123,5.009,5.009,5.009), na.rm = FALSE, stable.algo = FALSE), fvar(c(2.456,4.123,5.009), w = c(NA,1,3), na.rm = FALSE, stable.algo = 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(fvar(v, f, stable.algo = FALSE), fvar(vs, fs, w, stable.algo = FALSE))
  expect_equal(fvar(v, f, na.rm = FALSE, stable.algo = FALSE), fvar(vs, fs, w, na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(v2, f, stable.algo = FALSE), fvar(v2s, fs, w, stable.algo = FALSE))
  expect_equal(fvar(v2, f, na.rm = FALSE, stable.algo = FALSE), fvar(v2s, fs, w, na.rm = FALSE, stable.algo = FALSE))
  v[c(3,9)] <- NA; vs[c(2,5)] <- NA
  expect_equal(fvar(v, f, stable.algo = FALSE), fvar(vs, fs, w, stable.algo = FALSE))
  expect_equal(fvar(v, f, na.rm = FALSE, stable.algo = FALSE), fvar(vs, fs, w, na.rm = FALSE, stable.algo = FALSE))
  vs[c(2,5)] <- 4; w[c(2,5)] <- NA
  expect_equal(fvar(v, f, stable.algo = FALSE), fvar(vs, fs, w, stable.algo = FALSE))
  expect_equal(fvar(v, f, na.rm = FALSE, stable.algo = FALSE), fvar(vs, fs, w, na.rm = FALSE, stable.algo = FALSE))
  w[c(2,5)] <- 1; v2[c(1:2,7:8)] <- NA; v2s[c(1,4)] <- NA
  expect_equal(fvar(v2, f, stable.algo = FALSE), fvar(v2s, fs, w, stable.algo = FALSE))
  expect_equal(fvar(v2, f, na.rm = FALSE, stable.algo = FALSE), fvar(v2s, fs, w, na.rm = FALSE, stable.algo = FALSE))
  v2s[c(1,4)] <- 2.456; w[c(1,4)] <- NA
  expect_equal(fvar(v2, f, stable.algo = FALSE), fvar(v2s, fs, w, stable.algo = FALSE))
  expect_equal(fvar(v2, f, na.rm = FALSE, stable.algo = FALSE), fvar(v2s, fs, w, na.rm = FALSE, stable.algo = FALSE))
})

test_that("fvar with direct algorithm performs like fvar with unit weights", {
  expect_equal(fvar(NA, stable.algo = FALSE), fvar(NA, w = 1, stable.algo = FALSE))
  expect_equal(fvar(NA, na.rm = FALSE, stable.algo = FALSE), fvar(NA, w = 1, na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(1, stable.algo = FALSE), fvar(1, w = 1, stable.algo = FALSE))
  expect_equal(fvar(1:3, stable.algo = FALSE), fvar(1:3, w = rep(1,3), stable.algo = FALSE))
  expect_equal(fvar(-1:1, stable.algo = FALSE), fvar(-1:1, w = rep(1,3), stable.algo = FALSE))
  expect_equal(fvar(1, na.rm = FALSE, stable.algo = FALSE), fvar(1, w = 1, na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(1:3, na.rm = FALSE, stable.algo = FALSE), fvar(1:3, w = rep(1, 3), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(-1:1, na.rm = FALSE, stable.algo = FALSE), fvar(-1:1, w = rep(1, 3), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(x, stable.algo = FALSE), fvar(x, w = rep(1,100), stable.algo = FALSE))
  expect_equal(fvar(x, na.rm = FALSE, stable.algo = FALSE), fvar(x, w = rep(1, 100), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(xNA, na.rm = FALSE, stable.algo = FALSE), fvar(xNA, w = rep(1, 100), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(xNA, stable.algo = FALSE), fvar(xNA, w = rep(1, 100), stable.algo = FALSE))
  expect_equal(fvar(m, stable.algo = FALSE), fvar(m, w = rep(1, 32), stable.algo = FALSE))
  expect_equal(fvar(m, na.rm = FALSE, stable.algo = FALSE), fvar(m, w = rep(1, 32), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(mNA, na.rm = FALSE, stable.algo = FALSE), fvar(mNA, w = rep(1, 32), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(mNA, stable.algo = FALSE), fvar(mNA, w = rep(1, 32), stable.algo = FALSE))
  expect_equal(fvar(mtcars, stable.algo = FALSE), fvar(mtcars, w = rep(1, 32), stable.algo = FALSE))
  expect_equal(fvar(mtcars, na.rm = FALSE, stable.algo = FALSE), fvar(mtcars, w = rep(1, 32), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(mtcNA, na.rm = FALSE, stable.algo = FALSE), fvar(mtcNA, w = rep(1, 32), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(mtcNA, stable.algo = FALSE), fvar(mtcNA, w = rep(1, 32), stable.algo = FALSE))
  expect_equal(fvar(x, f, stable.algo = FALSE), fvar(x, f, rep(1,100), stable.algo = FALSE))
  expect_equal(fvar(x, f, na.rm = FALSE, stable.algo = FALSE), fvar(x, f, rep(1,100), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(xNA, f, na.rm = FALSE, stable.algo = FALSE), fvar(xNA, f, rep(1,100), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(xNA, f, stable.algo = FALSE), fvar(xNA, f, rep(1,100), stable.algo = FALSE))
  expect_equal(fvar(m, g, stable.algo = FALSE), fvar(m, g, rep(1,32), stable.algo = FALSE))
  expect_equal(fvar(m, g, na.rm = FALSE, stable.algo = FALSE), fvar(m, g, rep(1,32), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(mNA, g, na.rm = FALSE, stable.algo = FALSE), fvar(mNA, g, rep(1,32), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(mNA, g, stable.algo = FALSE), fvar(mNA, g, rep(1,32), stable.algo = FALSE))
  expect_equal(fvar(mtcars, g, stable.algo = FALSE), fvar(mtcars, g, rep(1,32), stable.algo = FALSE))
  expect_equal(fvar(mtcars, g, na.rm = FALSE, stable.algo = FALSE), fvar(mtcars, g, rep(1,32), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(mtcNA, g, na.rm = FALSE, stable.algo = FALSE), fvar(mtcNA, g, rep(1,32), na.rm = FALSE, stable.algo = FALSE))
  expect_equal(fvar(mtcNA, g, stable.algo = FALSE), fvar(mtcNA, g, rep(1,32), stable.algo = FALSE))
})

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

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

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

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

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

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

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


# fsd (not necessary to test in the same way because it's just sqrt(fvar()))

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

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

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

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