context("fscale / STD")



bsum <- base::sum

# TODO: Still a few uneccessary infinity values generated with weights when the sd is null. search replace_Inf to find them.

# 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"

bscale <- function(x, na.rm = FALSE, mean = 0, sd = 1) {
  if(na.rm || !anyNA(x)) `attributes<-`(drop(base::scale(x)), NULL) * sd + mean else
    rep(NA_real_, length(x))
}
# NOTE: This is what fscale currently does: If missing values, compute weighted mean and sd on available obs, and scale x using it. but don't insert additional missing values in x for missing weights ..
wbscale <- function(x, w, na.rm = FALSE, mean = 0, sd = 1) {
  if(na.rm) {
    x2 <- x
    cc <- complete.cases(x, w)
    x <- x[cc]
    # if(length(x) < 2L) return(rep(NA_real_, length(x2))) # wbscale(NA, 1, na.rm = TRUE) gives length 0
    if(length(x) < 2L || all(x[1L] == x[-1L])) return(rep(NA_real_, length(x2)))
    w <- w[cc]
  } else {
    if(length(x) < 2L) return(NA_real_)
    ck <- all(x[1L] == x[-1L])
    if(is.na(ck) || all(ck)) return(rep(NA_real_, length(x)))
  }
  sw <- bsum(w)
  wm <- bsum(w * x) / sw
  xdm <- x - wm
  wsd <- sqrt(bsum(w * xdm^2) / (sw - 1)) / sd
  if(!na.rm) return(xdm / wsd + mean)
  return((x2 - wm) / wsd + mean)
}


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

su <- function(x) if(is.null(dim(x))) `attributes<-`(qsu.default(x)[2:3], NULL) else `attributes<-`(qsu(x)[,2:3], NULL)
suby <- function(x, f) if(is.null(dim(x))) `attributes<-`(qsu.default(x, f)[, 2:3], NULL) else `attributes<-`(qsu(x, f)[,2:3,], NULL)

miss <- unname(rep(ifelse(dapply(mNA, anyNA), NA_real_, 0), 2))

test_that("Unweighted customized scaling works as intended", {
  expect_equal(su(fscale(x, mean = 5.1, sd = 3.9)), c(5.1, 3.9))
  expect_equal(su(fscale(x, mean = 5.1, sd = 3.9, na.rm = FALSE)), c(5.1, 3.9))
  expect_equal(su(fscale(xNA, mean = 5.1, sd = 3.9, na.rm = FALSE)), c(NaN, NA))
  expect_equal(su(fscale(xNA, mean = 5.1, sd = 3.9)), c(5.1, 3.9))
  expect_equal(qM(fscale(mtcars, mean = 5.1, sd = 3.9)), fscale(m, mean = 5.1, sd = 3.9))
  expect_equal(su(fscale(m, mean = 5.1, sd = 3.9)), rep(c(5.1, 3.9), c(11, 11)))
  expect_equal(su(fscale(m, mean = 5.1, sd = 3.9, na.rm = FALSE)), rep(c(5.1, 3.9), c(11, 11)))
  expect_equal(su(fscale(mNA, mean = 5.1, sd = 3.9, na.rm = FALSE)), rep(c(5.1, 3.9), c(11, 11))+miss)
  expect_equal(su(fscale(mNA, mean = 5.1, sd = 3.9)), rep(c(5.1, 3.9), c(11, 11)))
  expect_equal(suby(fscale(x, f, mean = 5.1, sd = 3.9), f), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))
  expect_equal(suby(fscale(x, f, mean = 5.1, sd = 3.9, na.rm = FALSE), f), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))
  expect_equal(suby(fscale(xNA, f, mean = 5.1, sd = 3.9), f), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))
})

test_that("Unweighted customized scaling works like bscale (defined above)", {
  expect_equal(fscale(x, mean = 5.1, sd = 3.9), bscale(x, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, mean = 5.1, sd = 3.9, na.rm = FALSE), bscale(x, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, mean = 5.1, sd = 3.9, na.rm = FALSE), bscale(xNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, mean = 5.1, sd = 3.9), bscale(xNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(qM(fscale(mtcars, mean = 5.1, sd = 3.9)), fscale(m, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, mean = 5.1, sd = 3.9), dapply(m, bscale, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, mean = 5.1, sd = 3.9, na.rm = FALSE), dapply(m, bscale, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, mean = 5.1, sd = 3.9, na.rm = FALSE), dapply(mNA, bscale, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, mean = 5.1, sd = 3.9), dapply(mNA, bscale, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, f, mean = 5.1, sd = 3.9), BY(x, f, bscale, na.rm = TRUE, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, f, mean = 5.1, sd = 3.9, na.rm = FALSE), BY(x, f, bscale, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, f, mean = 5.1, sd = 3.9, na.rm = FALSE), BY(xNA, f, bscale, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, f, mean = 5.1, sd = 3.9), BY(xNA, f, bscale, na.rm = TRUE, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, g, mean = 5.1, sd = 3.9), BY(m, g, bscale, na.rm = TRUE, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, g, mean = 5.1, sd = 3.9, na.rm = FALSE), BY(m, g, bscale, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, g, mean = 5.1, sd = 3.9, na.rm = FALSE), BY(mNA, g, bscale, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, g, mean = 5.1, sd = 3.9), BY(mNA, g, bscale, na.rm = TRUE, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, g, mean = 5.1, sd = 3.9), BY(mtcars, g, bscale, na.rm = TRUE, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, g, mean = 5.1, sd = 3.9, na.rm = FALSE), BY(mtcars, g, bscale, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, g, mean = 5.1, sd = 3.9, na.rm = FALSE), BY(mtcNA, g, bscale, use.g.names = FALSE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, g, mean = 5.1, sd = 3.9), BY(mtcNA, g, bscale, na.rm = TRUE, use.g.names = FALSE, mean = 5.1, sd = 3.9))
})

test_that("Unweighted customized scaling special cases perform as intended ", {
  # No mean / centering
  expect_equal(fscale(x, mean = FALSE, sd = 3.9), bscale(x, na.rm = TRUE, mean = fmean(x), sd = 3.9))
  expect_equal(fscale(x, mean = FALSE, sd = 3.9, na.rm = FALSE), bscale(x, mean = fmean(x), sd = 3.9))
  expect_equal(fscale(xNA, mean = FALSE, sd = 3.9, na.rm = FALSE), bscale(xNA, mean = fmean(xNA), sd = 3.9))
  expect_equal(fscale(xNA, mean = FALSE, sd = 3.9), bscale(xNA, na.rm = TRUE, mean = fmean(xNA), sd = 3.9))
  expect_equal(qM(fscale(mtcars, mean = FALSE, sd = 3.9)), fscale(m, mean = FALSE, sd = 3.9))
  expect_equal(fmean(fscale(mtcars, mean = FALSE, sd = 3.9)), fmean(mtcars))
  expect_equal(unname(fsd(fscale(mtcars, mean = FALSE, sd = 3.9))), rep(3.9, length(mtcars)))

  expect_equal(fscale(x, f, mean = FALSE), BY(x, f, bscale, na.rm = TRUE, use.g.names = FALSE) + B(x, f))
  expect_equal(fscale(x, f, na.rm = FALSE, mean = FALSE), BY(x, f, bscale, use.g.names = FALSE) + B(x, f))
  expect_equal(fscale(xNA, f, na.rm = FALSE, mean = FALSE), BY(xNA, f, bscale, use.g.names = FALSE) + B(xNA, f))
  expect_equal(fscale(xNA, f, mean = FALSE), BY(xNA, f, bscale, na.rm = TRUE, use.g.names = FALSE) + B(xNA, f))
  expect_equal(fscale(m, g, mean = FALSE), BY(m, g, bscale, na.rm = TRUE, use.g.names = FALSE) + B(m, g))
  expect_equal(fscale(m, g, na.rm = FALSE, mean = FALSE), BY(m, g, bscale, use.g.names = FALSE) + B(m, g))
  expect_equal(fscale(mNA, g, na.rm = FALSE, mean = FALSE), BY(mNA, g, bscale, use.g.names = FALSE) + B(mNA, g))
  expect_equal(fscale(mNA, g, mean = FALSE), BY(mNA, g, bscale, na.rm = TRUE, use.g.names = FALSE) + B(mNA, g))
  expect_equal(fscale(mtcars, g, mean = FALSE), BY(mtcars, g, bscale, na.rm = TRUE, use.g.names = FALSE) + B(mtcars, g))
  expect_equal(fscale(mtcars, g, na.rm = FALSE, mean = FALSE), BY(mtcars, g, bscale, use.g.names = FALSE) + B(mtcars, g))
  expect_equal(fscale(mtcNA, g, na.rm = FALSE, mean = FALSE), BY(mtcNA, g, bscale, use.g.names = FALSE) + B(mtcNA, g))
  expect_equal(fscale(mtcNA, g, mean = FALSE), BY(mtcNA, g, bscale, na.rm = TRUE, use.g.names = FALSE)+ B(mtcNA, g))

  # Centering on overall mean
  expect_equal(fscale(x, f, mean = "overall.mean"), BY(x, f, bscale, na.rm = TRUE, use.g.names = FALSE) + ave(x))
  expect_equal(fscale(x, f, na.rm = FALSE, mean = "overall.mean"), BY(x, f, bscale, use.g.names = FALSE) + ave(x))
  # expect_equal(fscale(xNA, f, na.rm = FALSE, mean = "overall.mean"), BY(xNA, f, bscale, use.g.names = FALSE) + B(xNA)) # Not the same !!
  expect_equal(fscale(xNA, f, mean = "overall.mean"), BY(xNA, f, bscale, na.rm = TRUE, use.g.names = FALSE) + B(xNA))
  expect_equal(fscale(m, g, mean = "overall.mean"), BY(m, g, bscale, na.rm = TRUE, use.g.names = FALSE) + B(m))
  expect_equal(fscale(m, g, na.rm = FALSE, mean = "overall.mean"), BY(m, g, bscale, use.g.names = FALSE) + B(m))
  # expect_equal(fscale(mNA, g, na.rm = FALSE, mean = "overall.mean"), BY(mNA, g, bscale, use.g.names = FALSE) + B(mNA))
  expect_equal(fscale(mNA, g, mean = "overall.mean"), BY(mNA, g, bscale, na.rm = TRUE, use.g.names = FALSE) + B(mNA))
  expect_equal(fscale(mtcars, g, mean = "overall.mean"), BY(mtcars, g, bscale, na.rm = TRUE, use.g.names = FALSE) + B(mtcars))
  expect_equal(fscale(mtcars, g, na.rm = FALSE, mean = "overall.mean"), BY(mtcars, g, bscale, use.g.names = FALSE) + B(mtcars))
  # expect_equal(fscale(mtcNA, g, na.rm = FALSE, mean = "overall.mean"), BY(mtcNA, g, bscale, use.g.names = FALSE) + B(mtcNA))
  expect_equal(fscale(mtcNA, g, mean = "overall.mean"), BY(mtcNA, g, bscale, na.rm = TRUE, use.g.names = FALSE)+ B(mtcNA))

  # Scaling by within-sd
  expect_equal(fscale(x, f, sd = "within.sd"), BY(x, f, bscale, na.rm = TRUE, use.g.names = FALSE) * fsd(W(x, f)))
  expect_equal(fscale(x, f, na.rm = FALSE, sd = "within.sd"), BY(x, f, bscale, use.g.names = FALSE) * fsd(W(x, f)))
  # expect_equal(fscale(xNA, f, na.rm = FALSE, sd = "within.sd"), BY(xNA, f, bscale, use.g.names = FALSE) * fsd(W(xNA, f))) # Not the same !!
  expect_equal(fscale(xNA, f, sd = "within.sd"), BY(xNA, f, bscale, na.rm = TRUE, use.g.names = FALSE) * fsd(W(xNA, f)))
  expect_equal(fscale(m, g, sd = "within.sd"), BY(m, g, bscale, na.rm = TRUE, use.g.names = FALSE) * TRA(m,fsd(W(m, g)),1L))
  expect_equal(fscale(m, g, na.rm = FALSE, sd = "within.sd"), BY(m, g, bscale, use.g.names = FALSE) * TRA(m,fsd(W(m, g)),1L))
  # expect_equal(fscale(mNA, g, na.rm = FALSE, sd = "within.sd"), BY(mNA, g, bscale, use.g.names = FALSE) * TRA(mNA,fsd(W(mNA, g)),1L))
  expect_equal(fscale(mNA, g, sd = "within.sd"), BY(mNA, g, bscale, na.rm = TRUE, use.g.names = FALSE) * TRA(mNA,fsd(W(mNA, g)),1L))
  expect_equal(fscale(mtcars, g, sd = "within.sd"), BY(mtcars, g, bscale, na.rm = TRUE, use.g.names = FALSE) * TRA(mtcars,fsd(W(mtcars, g)),1L))
  expect_equal(fscale(mtcars, g, na.rm = FALSE, sd = "within.sd"), BY(mtcars, g, bscale, use.g.names = FALSE) * TRA(mtcars,fsd(W(mtcars, g)),1L))
  # expect_equal(fscale(mtcNA, g, na.rm = FALSE, sd = "within.sd"), BY(mtcNA, g, bscale, use.g.names = FALSE) * TRA(mtcNA,fsd(W(mtcNA, g)),1L))
  expect_equal(fscale(mtcNA, g, sd = "within.sd"), BY(mtcNA, g, bscale, na.rm = TRUE, use.g.names = FALSE) * TRA(mtcNA,fsd(W(mtcNA, g)),1L))

  # Centering on overall mean and scaling by within-sd
  expect_equal(fscale(x, f, mean = "overall.mean", sd = "within.sd"), BY(x, f, bscale, na.rm = TRUE, use.g.names = FALSE) * fsd(W(x, f)) + ave(x))
  expect_equal(fscale(x, f, na.rm = FALSE, mean = "overall.mean", sd = "within.sd"), BY(x, f, bscale, use.g.names = FALSE) * fsd(W(x, f)) + ave(x))
  # expect_equal(fscale(xNA, f, na.rm = FALSE, mean = "overall.mean", sd = "within.sd"), BY(xNA, f, bscale, use.g.names = FALSE) * fsd(W(xNA, f)) + B(xNA)) # Not the same !!
  expect_equal(fscale(xNA, f, mean = "overall.mean", sd = "within.sd"), BY(xNA, f, bscale, na.rm = TRUE, use.g.names = FALSE) * fsd(W(xNA, f)) + B(xNA))
  expect_equal(fscale(m, g, mean = "overall.mean", sd = "within.sd"), BY(m, g, bscale, na.rm = TRUE, use.g.names = FALSE) * TRA(m,fsd(W(m, g)),1L) + B(m))
  expect_equal(fscale(m, g, na.rm = FALSE, mean = "overall.mean", sd = "within.sd"), BY(m, g, bscale, use.g.names = FALSE) * TRA(m,fsd(W(m, g)),1L) + B(m))
  # expect_equal(fscale(mNA, g, na.rm = FALSE, mean = "overall.mean", sd = "within.sd"), BY(mNA, g, bscale, use.g.names = FALSE) * TRA(mNA,fsd(W(mNA, g)),1L) + B(mNA))
  expect_equal(fscale(mNA, g, mean = "overall.mean", sd = "within.sd"), BY(mNA, g, bscale, na.rm = TRUE, use.g.names = FALSE) * TRA(mNA,fsd(W(mNA, g)),1L) + B(mNA))
  expect_equal(fscale(mtcars, g, mean = "overall.mean", sd = "within.sd"), BY(mtcars, g, bscale, na.rm = TRUE, use.g.names = FALSE) * TRA(mtcars,fsd(W(mtcars, g)),1L) + B(mtcars))
  expect_equal(fscale(mtcars, g, na.rm = FALSE, mean = "overall.mean", sd = "within.sd"), BY(mtcars, g, bscale, use.g.names = FALSE) * TRA(mtcars,fsd(W(mtcars, g)),1L) + B(mtcars))
  # expect_equal(fscale(mtcNA, g, na.rm = FALSE, mean = "overall.mean", sd = "within.sd"), BY(mtcNA, g, bscale, use.g.names = FALSE) * TRA(mtcNA,fsd(W(mtcNA, g)),1L) + B(mtcNA))
  expect_equal(fscale(mtcNA, g, mean = "overall.mean", sd = "within.sd"), BY(mtcNA, g, bscale, na.rm = TRUE, use.g.names = FALSE) * TRA(mtcNA,fsd(W(mtcNA, g)),1L) + B(mtcNA))

})

# Still test weighted special cases ...

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

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

wsu <- function(x, w) if(is.null(dim(x))) `attributes<-`(qsu.default(x, w = w)[2:3], NULL) else `attributes<-`(qsu(x, w = w)[,2:3], NULL)
wsuby <- function(x, f, w) if(is.null(dim(x))) `attributes<-`(qsu.default(x, f, w = w)[, 2:3], NULL) else `attributes<-`(qsu(x, f, w = w)[,2:3,], NULL)

test_that("Weighted customized scaling works as intended", {
  expect_equal(wsu(fscale(x, w = w, mean = 5.1, sd = 3.9), w = w), c(5.1, 3.9))
  expect_equal(wsu(fscale(x, w = w, mean = 5.1, sd = 3.9, na.rm = FALSE), w = w), c(5.1, 3.9))
  expect_equal(wsu(fscale(xNA, w = w, mean = 5.1, sd = 3.9, na.rm = FALSE), w = w), c(NaN, NA))
  expect_equal(wsu(fscale(xNA, w = w, mean = 5.1, sd = 3.9), w = w), c(5.1, 3.9), w = w)
  expect_equal(qM(fscale(mtcars, w = wdat, mean = 5.1, sd = 3.9)), fscale(m, w = wdat, mean = 5.1, sd = 3.9))
  expect_equal(wsu(fscale(m, w = wdat, mean = 5.1, sd = 3.9), w = wdat), rep(c(5.1, 3.9), c(11, 11)))
  expect_equal(wsu(fscale(m, w = wdat, mean = 5.1, sd = 3.9, na.rm = FALSE), w = wdat), rep(c(5.1, 3.9), c(11, 11)))
  expect_equal(wsu(fscale(mNA, w = wdat, mean = 5.1, sd = 3.9, na.rm = FALSE), w = wdat), rep(c(5.1, 3.9), c(11, 11))+miss)
  expect_equal(wsu(fscale(mNA, w = wdat, mean = 5.1, sd = 3.9), w = wdat), rep(c(5.1, 3.9), c(11, 11)))
  expect_equal(wsuby(fscale(x, f, w = w, mean = 5.1, sd = 3.9), f, w = w), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))
  expect_equal(wsuby(fscale(x, f, w = w, mean = 5.1, sd = 3.9, na.rm = FALSE), f, w = w), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))
  expect_equal(wsuby(fscale(xNA, f, w = w, mean = 5.1, sd = 3.9), f, w = w), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))

  expect_equal(wsu(fscale(x, w = w, mean = FALSE, sd = 3.9), w = w), c(fmean(x, w = w), 3.9))
  expect_equal(wsu(fscale(x, w = w, mean = FALSE, sd = 3.9, na.rm = FALSE), w = w), c(fmean(x, w = w), 3.9))
  expect_equal(wsu(fscale(xNA, w = w, mean = FALSE, sd = 3.9, na.rm = FALSE), w = w), c(NaN, NA))
  expect_equal(wsu(fscale(xNA, w = w, mean = FALSE, sd = 3.9), w = w), c(fmean(xNA, w = w), 3.9))
  expect_equal(qM(fscale(mtcars, w = wdat, mean = FALSE, sd = 3.9)), fscale(m, w = wdat, mean = FALSE, sd = 3.9))

  # ...
  # expect_equal(wsuby(fscale(x, f, w = w, mean = "overall.mean", sd = 3.9), f, w = w), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))
  # expect_equal(wsuby(fscale(x, f, w = w, mean = FALSE, sd = 3.9, na.rm = FALSE), f, w = w), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))
  # expect_equal(wsuby(fscale(xNA, f, w = w, mean = FALSE, sd = 3.9), f, w = w), rep(c(5.1, 3.9), rep(fnlevels(f), 2)))

})

test_that("Weighted customized scaling performs like wbscale (defined above)", {
  # complete weights
  expect_equal(fscale(NA, w = 1, mean = 5.1, sd = 3.9), wbscale(NA, 1, mean = 5.1, sd = 3.9))
  expect_equal(fscale(NA, w = 1, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(NA, 1, mean = 5.1, sd = 3.9))
  expect_equal(fscale(1, w = 1, mean = 5.1, sd = 3.9), wbscale(1, w = 1, mean = 5.1, sd = 3.9))
  expect_equal(fscale(1:3, w = 1:3, mean = 5.1, sd = 3.9), wbscale(1:3, 1:3, mean = 5.1, sd = 3.9))
  expect_equal(fscale(-1:1, w = 1:3, mean = 5.1, sd = 3.9), wbscale(-1:1, 1:3, mean = 5.1, sd = 3.9))
  expect_equal(fscale(1, w = 1, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(1, 1, mean = 5.1, sd = 3.9))
  expect_equal(fscale(1:3, w = c(0.99,3454,1.111), na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(1:3, c(0.99,3454,1.111), mean = 5.1, sd = 3.9))
  expect_equal(fscale(-1:1, w = 1:3, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(-1:1, 1:3, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, w = w, mean = 5.1, sd = 3.9), wbscale(x, w, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, w = w, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(x, w, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, w = w, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(xNA, w, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, w = w, mean = 5.1, sd = 3.9), wbscale(xNA, w, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(qM(fscale(mtcars, w = wdat, mean = 5.1, sd = 3.9)), fscale(m, w = wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, w = wdat, mean = 5.1, sd = 3.9), dapply(m, wbscale, wdat, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, w = wdat, na.rm = FALSE, mean = 5.1, sd = 3.9), dapply(m, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, w = wdat, na.rm = FALSE, mean = 5.1, sd = 3.9), dapply(mNA, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, w = wdat, mean = 5.1, sd = 3.9), dapply(mNA, wbscale, wdat, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, w = wdat, mean = 5.1, sd = 3.9), dapply(mtcars, wbscale, wdat, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, w = wdat, na.rm = FALSE, mean = 5.1, sd = 3.9), dapply(mtcars, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, w = wdat, na.rm = FALSE, mean = 5.1, sd = 3.9), dapply(mtcNA, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, w = wdat, mean = 5.1, sd = 3.9), dapply(mtcNA, wbscale, wdat, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, f, w, mean = 5.1, sd = 3.9), BY(x, f, wbscale, w, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, f, w, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(x, f, wbscale, w, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, f, w, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(xNA, f, wbscale, w, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, f, w, mean = 5.1, sd = 3.9), BY(xNA, f, wbscale, w, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, g, wdat, mean = 5.1, sd = 3.9), BY(m, g, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, g, wdat, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(m, g, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, g, wdat, na.rm = FALSE, mean = 5.1, sd = 3.9),  BY(mNA, g, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, g, wdat, mean = 5.1, sd = 3.9), BY(mNA, g, wbscale, wdat, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, g, wdat, mean = 5.1, sd = 3.9), BY(mtcars, g, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, g, wdat, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(mtcars, g, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, g, wdat, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(mtcNA, g, wbscale, wdat, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, g, wdat, mean = 5.1, sd = 3.9), BY(mtcNA, g, wbscale, wdat, na.rm = TRUE, mean = 5.1, sd = 3.9))
  # missing weights
  expect_equal(fscale(NA, w = NA, mean = 5.1, sd = 3.9), wbscale(NA, NA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(NA, w = NA, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(NA, NA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(1, w = NA, mean = 5.1, sd = 3.9), wbscale(1, w = NA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(1:3, w = c(NA,1:2), mean = 5.1, sd = 3.9), wbscale(1:3, c(NA,1:2), na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(-1:1, w = c(NA,1:2), mean = 5.1, sd = 3.9), wbscale(-1:1, c(NA,1:2), na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(1, w = NA, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(1, NA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(1:3, w = c(NA,1:2), na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(1:3, c(NA,1:2), mean = 5.1, sd = 3.9))
  expect_equal(fscale(-1:1, w = c(NA,1:2), na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(-1:1, c(NA,1:2), mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, w = wNA, mean = 5.1, sd = 3.9), wbscale(x, wNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, w = wNA, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(x, wNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, w = wNA, na.rm = FALSE, mean = 5.1, sd = 3.9), wbscale(xNA, wNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, w = wNA, mean = 5.1, sd = 3.9), wbscale(xNA, wNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(qM(fscale(mtcars, w = wdatNA, mean = 5.1, sd = 3.9)), fscale(m, w = wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, w = wdatNA, mean = 5.1, sd = 3.9), dapply(m, wbscale, wdatNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, w = wdatNA, na.rm = FALSE, mean = 5.1, sd = 3.9), dapply(m, wbscale, wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, w = wdatNA, na.rm = FALSE, mean = 5.1, sd = 3.9), dapply(mNA, wbscale, wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, w = wdatNA, mean = 5.1, sd = 3.9), dapply(mNA, wbscale, wdatNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, w = wdatNA, mean = 5.1, sd = 3.9), dapply(mtcars, wbscale, wdatNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, w = wdatNA, na.rm = FALSE, mean = 5.1, sd = 3.9), dapply(mtcars, wbscale, wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, w = wdatNA, na.rm = FALSE, mean = 5.1, sd = 3.9), dapply(mtcNA, wbscale, wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, w = wdatNA, mean = 5.1, sd = 3.9), dapply(mtcNA, wbscale, wdatNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, f, wNA, mean = 5.1, sd = 3.9), BY(x, f, wbscale, wNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(x, f, wNA, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(x, f, wbscale, wNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(xNA, f, wNA, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(xNA, f, wbscale, wNA, mean = 5.1, sd = 3.9))
  # expect_equal(fscale(xNA, f, wNA, mean = 5.1, sd = 3.9), BY(xNA, f, wbscale, wNA, na.rm = TRUE, mean = 5.1, sd = 3.9)) # failed on release-windows-ix86+x86_64
  expect_equal(replace_Inf(fscale(m, g, wdatNA, mean = 5.1, sd = 3.9), NA), BY(m, g, wbscale, wdatNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(m, g, wdatNA, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(m, g, wbscale, wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mNA, g, wdatNA, na.rm = FALSE, mean = 5.1, sd = 3.9),  BY(mNA, g, wbscale, wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(replace_Inf(fscale(mNA, g, wdatNA, mean = 5.1, sd = 3.9), NA), BY(mNA, g, wbscale, wdatNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(replace_Inf(fscale(mtcars, g, wdatNA, mean = 5.1, sd = 3.9), NA), BY(mtcars, g, wbscale, wdatNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcars, g, wdatNA, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(mtcars, g, wbscale, wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(fscale(mtcNA, g, wdatNA, na.rm = FALSE, mean = 5.1, sd = 3.9), BY(mtcNA, g, wbscale, wdatNA, mean = 5.1, sd = 3.9))
  expect_equal(replace_Inf(fscale(mtcNA, g, wdatNA, mean = 5.1, sd = 3.9), NA), BY(mtcNA, g, wbscale, wdatNA, na.rm = TRUE, mean = 5.1, sd = 3.9))
})

test_that("Weighted customized scaling special cases perform as intended ", {
  # NOTE: These tests are currently only run with complete weights. STill implement them for missing weights ...

  # No mean / centering
  expect_equal(fscale(x, w = w, mean = FALSE, sd = 3.9), wbscale(x, na.rm = TRUE, w = w, mean = fmean(x, w = w), sd = 3.9))
  expect_equal(fscale(x, w = w, mean = FALSE, sd = 3.9, na.rm = FALSE), wbscale(x, w = w, mean = fmean(x, w = w), sd = 3.9))
  expect_equal(fscale(xNA, w = w, mean = FALSE, sd = 3.9, na.rm = FALSE), wbscale(xNA, w = w, mean = fmean(xNA, w = w), sd = 3.9))
  expect_equal(fscale(xNA, w = w, mean = FALSE, sd = 3.9), wbscale(xNA, na.rm = TRUE, w = w, mean = fmean(xNA, w = w), sd = 3.9))
  expect_equal(qM(fscale(mtcars, w = wdat, mean = FALSE, sd = 3.9)), fscale(m, w = wdat, mean = FALSE, sd = 3.9))
  expect_equal(fmean(fscale(mtcars, w = wdat, mean = FALSE, sd = 3.9), w = wdat), fmean(mtcars, w = wdat))
  expect_equal(unname(fsd(fscale(mtcars, w = wdat, mean = FALSE, sd = 3.9), w = wdat)), rep(3.9, length(mtcars)))

  expect_equal(fscale(x, f, w, mean = FALSE), unlist(Map(wbscale, split(x, f), split(w, f), na.rm = TRUE), use.names = FALSE) + B(x, f, w))
  expect_equal(fscale(x, f, w, mean = FALSE, na.rm = FALSE), unlist(Map(wbscale, split(x, f), split(w, f)), use.names = FALSE) + B(x, f, w))
  expect_equal(fscale(xNA, f, w, mean = FALSE, na.rm = FALSE), unlist(Map(wbscale, split(xNA, f), split(w, f)), use.names = FALSE) + B(xNA, f, w))
  expect_equal(fscale(xNA, f, w, mean = FALSE), unlist(Map(wbscale, split(xNA, f), split(w, f), na.rm = TRUE), use.names = FALSE) + B(xNA, f, w))

  # Centering on overall mean
  expect_equal(fscale(x, f, w, mean = "overall.mean"), unlist(Map(wbscale, split(x, f), split(w, f), na.rm = TRUE), use.names = FALSE) + B(x, w = w))
  expect_equal(fscale(x, f, w, mean = "overall.mean", na.rm = FALSE), unlist(Map(wbscale, split(x, f), split(w, f)), use.names = FALSE) + B(x, w = w))
  # expect_equal(fscale(xNA, f, w, mean = "overall.mean", na.rm = FALSE), unlist(Map(wbscale, split(xNA, f), split(w, f)), use.names = FALSE) + B(xNA, w = w)) # Not the same !!
  expect_equal(fscale(xNA, f, w, mean = "overall.mean"), unlist(Map(wbscale, split(xNA, f), split(w, f), na.rm = TRUE), use.names = FALSE) + B(xNA, w = w))

  # Scaling by within-sd
  expect_equal(fscale(x, f, w, sd = "within.sd"), unlist(Map(wbscale, split(x, f), split(w, f), na.rm = TRUE), use.names = FALSE) * fsd(W(x, f, w), w = w))
  expect_equal(fscale(x, f, w, sd = "within.sd", na.rm = FALSE), unlist(Map(wbscale, split(x, f), split(w, f)), use.names = FALSE) * fsd(W(x, f, w), w = w))
  # expect_equal(fscale(xNA, f, w, sd = "within.sd", na.rm = FALSE), unlist(Map(wbscale, split(xNA, f), split(w, f)), use.names = FALSE) * fsd(W(xNA, f, w), w = w)) # Not the same !!
  expect_equal(fscale(xNA, f, w, sd = "within.sd"), unlist(Map(wbscale, split(xNA, f), split(w, f), na.rm = TRUE), use.names = FALSE) * fsd(W(xNA, f, w), w = w))

  # Centering on overall mean and scaling by within-sd
  expect_equal(fscale(x, f, w, mean = "overall.mean", sd = "within.sd"), unlist(Map(wbscale, split(x, f), split(w, f), na.rm = TRUE), use.names = FALSE) * fsd(W(x, f, w), w = w) + B(x, w = w))
  expect_equal(fscale(x, f, w, mean = "overall.mean", sd = "within.sd", na.rm = FALSE), unlist(Map(wbscale, split(x, f), split(w, f)), use.names = FALSE) * fsd(W(x, f, w), w = w) + B(x, w = w))
  # expect_equal(fscale(xNA, f, w, mean = "overall.mean", sd = "within.sd", na.rm = FALSE), unlist(Map(wbscale, split(xNA, f), split(w, f)), use.names = FALSE) * fsd(W(xNA, f, w), w = w) + B(xNA, w = w)) # Not the same !!
  expect_equal(fscale(xNA, f, w, mean = "overall.mean", sd = "within.sd"), unlist(Map(wbscale, split(xNA, f), split(w, f), na.rm = TRUE), use.names = FALSE) * fsd(W(xNA, f, w), w = w) + B(xNA, w = w))


})

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

test_that("fscale customized scaling performs numerically stable", {
  expect_true(all_obj_equal(replicate(50, fscale(1, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(NA, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(NA, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(x, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(x, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(xNA, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(xNA, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(m, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(m, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mNA, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mNA, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mtcars, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mtcars, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mtcNA, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mtcNA, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(x, f, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(x, f, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(xNA, f, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(xNA, f, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(m, g, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(m, g, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mNA, g, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mNA, g, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mtcars, g, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mtcars, g, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mtcNA, g, na.rm = FALSE, mean = 5.1, sd = 3.9), simplify = FALSE)))
  expect_true(all_obj_equal(replicate(50, fscale(mtcNA, g, mean = 5.1, sd = 3.9), simplify = FALSE)))
})

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

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

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

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

# NOTE: fscale(c(a, a)) gives c(NaN, NaN) (sd is 0) !!!
test_that("fscale handles special values in the right way", {
  expect_equal(fscale(NA), NA_real_)
  expect_equal(fscale(NaN), NA_real_)
  expect_equal(fscale(Inf), NA_real_)
  expect_equal(fscale(-Inf), NA_real_)
  expect_equal(fscale(TRUE), NA_real_)
  expect_equal(fscale(FALSE), NA_real_)
  expect_equal(fscale(NA, na.rm = FALSE), NA_real_)
  expect_equal(fscale(NaN, na.rm = FALSE), NA_real_)
  expect_equal(fscale(Inf, na.rm = FALSE), NA_real_)
  expect_equal(fscale(-Inf, na.rm = FALSE), NA_real_)
  expect_equal(fscale(TRUE, na.rm = FALSE), NA_real_)
  expect_equal(fscale(FALSE, na.rm = FALSE), NA_real_)
  expect_equal(fscale(c(1,NA)), c(NA_real_,NA_real_))
  expect_equal(fscale(c(1,NaN)), c(NA_real_,NA_real_))
  expect_equal(fscale(c(1,Inf)), c(NA_real_,NA_real_))
  expect_equal(fscale(c(1,-Inf)), c(NA_real_,NA_real_))
  expect_equal(fscale(c(1,Inf), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fscale(c(1,-Inf), na.rm = FALSE), c(NA_real_,NA_real_))
  expect_equal(fscale(c(FALSE,FALSE), na.rm = FALSE), c(NaN,NaN))
  expect_equal(fscale(c(1,1), na.rm = FALSE), c(NaN,NaN))
})

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

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

test_that("fscale shoots errors for wrong input to mean and sd", {
  expect_error(fscale(x, sd = FALSE))
  expect_error(fscale(m, sd = FALSE))
  expect_error(fscale(mtcars, sd = FALSE))
  expect_error(fscale(x, sd = "bla"))
  expect_error(fscale(x, mean = "bla"))
  expect_error(fscale(x, sd = "within.sd"))
  expect_error(fscale(m, sd = "within.sd"))
  expect_error(fscale(mtcars, sd = "within.sd"))
  expect_error(fscale(x, mean = "overall.mean"))
  expect_error(fscale(m, mean = "overall.mean"))
  expect_error(fscale(mtcars, mean = "overall.mean"))
  expect_error(fscale(m, mean = fmean(m)))
  expect_error(fscale(mtcars, mean = fmean(mtcars)))
  expect_error(fscale(m, sd = fsd(m)))
  expect_error(fscale(mtcars, sd = fsd(mtcars)))
})




# Testing STD: Only testing wrong inputs, especially for data.frame method. Otherwise it is identical to fscale

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

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

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

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