test_that("the back-transform preserves the linear predictor exactly", {
  set.seed(31)
  n <- 50; p <- 4
  Xraw <- matrix(stats::rnorm(n * p, mean = c(0, 10, -5, 100),
                              sd = c(1, 0.01, 50, 3)), n, p, byrow = TRUE)
  x_mu <- colMeans(Xraw)
  x_sd <- apply(Xraw, 2L, stats::sd)
  Xstd <- sweep(sweep(Xraw, 2L, x_mu, "-"), 2L, x_sd, "/")

  B  <- matrix(stats::rnorm(7 * p), 7, p)   # draws on the standardized scale
  b0 <- stats::rnorm(7)

  bt <- bbqr:::.bbqr_backtransform(B, b0, list(x_mu = x_mu, x_sd = x_sd))

  # eta_std[d, i] and eta_raw[d, i] must agree for every draw and observation
  eta_std <- outer(b0, rep(1, n)) + B %*% t(Xstd)
  eta_raw <- outer(bt$beta0, rep(1, n)) + bt$beta %*% t(Xraw)
  expect_equal(eta_std, eta_raw)
})

test_that("a NULL scaler leaves the draws untouched", {
  B <- matrix(1:6, 3, 2); b0 <- c(1, 2, 3)
  bt <- bbqr:::.bbqr_backtransform(B, b0, NULL)
  expect_identical(bt$beta, B)
  expect_identical(bt$beta0, b0)
})

test_that("standardize rescues covariates on very different scales", {
  s <- sim_binary(n = 400, seed = 21)
  # x1 blown up 100-fold and x2 shrunk and shifted: the fixed priors are then
  # badly matched to the raw scale, which is exactly when standardizing pays.
  d <- s$data
  d$x1 <- d$x1 * 100
  d$x2 <- d$x2 / 50 + 7
  truth <- s$beta / c(100, 1 / 50, 1)   # truth in the rescaled units

  set.seed(3); raw <- coef(bbqr(s$formula, d, ndraw = 2500, burn = 500))
  set.seed(3); std <- coef(bbqr(s$formula, d, ndraw = 2500, burn = 500,
                                standardize = TRUE))

  expect_gt(cos_sim(std[-1], truth), 0.99)
  expect_gt(cos_sim(std[-1], truth), cos_sim(raw[-1], truth))
})

test_that("standardize is near-neutral when covariates share a scale", {
  s <- sim_binary(n = 400, seed = 25)
  set.seed(4); raw <- coef(bbqr(s$formula, s$data, ndraw = 2500, burn = 500))
  set.seed(4); std <- coef(bbqr(s$formula, s$data, ndraw = 2500, burn = 500,
                                standardize = TRUE))
  expect_gt(cos_sim(raw[-1], std[-1]), 0.98)
})

test_that("summary() reports intervals and a selection flag", {
  s <- sim_binary(n = 400, seed = 22)
  sm <- summary(bbqr(s$formula, s$data, ndraw = 2000, burn = 400))
  expect_s3_class(sm, "bbqr.summary")
  tab <- sm$coefficients
  expect_equal(nrow(tab), 4L)
  expect_true(all(tab$lower <= tab$mean))
  expect_true(all(tab$mean <= tab$upper))
  expect_type(tab$selected, "logical")
  # x1 and x2 are genuinely non-zero and should be picked up at n = 400
  expect_true(tab["x1", "selected"])
  expect_true(tab["x2", "selected"])
  expect_equal(sm$level, 0.95)

  sm90 <- summary(bbqr(s$formula, s$data, ndraw = 1000, burn = 200),
                  level = 0.90)
  expect_equal(sm90$level, 0.90)
  expect_error(summary(bbqr(s$formula, s$data, ndraw = 300, burn = 50),
                       level = 1.5))
})

test_that("predict() returns valid probabilities, indices and labels", {
  s <- sim_binary(n = 300, seed = 23)
  fit <- bbqr(s$formula, s$data, ndraw = 1500, burn = 300)

  p <- predict(fit, newdata = s$data, type = "response")
  expect_length(p, 300L)
  expect_true(all(p >= 0 & p <= 1))

  eta <- predict(fit, newdata = s$data, type = "link")
  expect_true(all(is.finite(eta)))

  cl <- predict(fit, newdata = s$data, type = "class")
  expect_true(all(cl %in% c(0L, 1L)))
  expect_identical(cl, as.integer(eta > 0))
  # should beat guessing the majority class
  expect_gt(mean(cl == s$data$y), max(table(s$data$y)) / 300)

  expect_error(predict(fit), "newdata")
  expect_error(predict(fit, newdata = s$data[, c("y", "x1")]))
})

## The reference for the response probability is the model itself, not a second
## copy of the closed form: P(y = 1 | x) is the mass above -eta of an asymmetric
## Laplace density with skewness tau and INVERSE scale sigma, integrated
## numerically. A copy of the formula would agree with any mistake in it, which
## is how an eta / sigma in place of eta * sigma once went unnoticed.
ald_tail_mass <- function(eta, tau, sigma) {
  dens <- function(u)
    sigma * tau * (1 - tau) * exp(-sigma * u * (tau - (u < 0)))
  area <- function(lower, upper)
    stats::integrate(dens, lower, upper, rel.tol = 1e-10, abs.tol = 1e-12)$value
  ## The density has a kink at zero; integrate each side of it separately.
  vapply(eta, function(e)
    if (e > 0) area(-e, 0) + area(0, Inf) else area(-e, Inf), numeric(1))
}

test_that("predict(type = 'response') is the ALD tail mass at sigma != 1", {
  s <- sim_binary(n = 200, seed = 24)
  tau <- 0.3
  fit <- bbqr(s$formula, s$data, quantile = tau, sigma_fixed = 2,
              ndraw = 800, burn = 200)
  expect_equal(unique(fit$sigma), 2)
  rows <- s$data[1:25, ]
  eta <- predict(fit, newdata = rows, type = "link")
  expect_true(any(eta > 0) && any(eta < 0))
  expect_equal(predict(fit, newdata = rows, type = "response"),
               ald_tail_mass(eta, tau, sigma = 2), tolerance = 1e-7)
  # the label switches where the probability crosses 1 - tau
  expect_equal(bbqr:::.bbqr_ald_prob(0, tau), 1 - tau)
})

test_that("the response probability is invariant under the scale map", {
  s <- sim_binary(n = 150, seed = 26)
  fit <- bbqr(s$formula, s$data, quantile = 0.4, anchor = "free",
              ndraw = 600, burn = 150)
  expect_gt(stats::sd(fit$sigma), 0)
  ## (beta0, beta, sigma) -> (c beta0, c beta, sigma / c) leaves the model
  ## unchanged, so it must leave every predicted probability unchanged.
  scaled <- fit
  scaled$beta0 <- 3 * fit$beta0
  scaled$beta  <- 3 * fit$beta
  scaled$sigma <- fit$sigma / 3
  for (est in c("plug_in_mean", "plug_in_median", "draw_average"))
    expect_equal(predict(scaled, newdata = s$data, estimator = est),
                 predict(fit, newdata = s$data, estimator = est),
                 info = est)
})

test_that("coef() returns medians on request and means by default", {
  s <- sim_binary(n = 200, seed = 27)
  fit <- bbqr(s$formula, s$data, ndraw = 600, burn = 100)
  expect_identical(coef(fit), coef(fit, statistic = "mean"))
  expect_identical(coef(fit),
                   c(`(Intercept)` = mean(fit$beta0), colMeans(fit$beta)))
  med <- coef(fit, statistic = "median")
  expect_named(med, names(coef(fit)))
  expect_equal(unname(med),
               c(stats::median(fit$beta0),
                 unname(apply(fit$beta, 2L, stats::median))))
  expect_error(coef(fit, statistic = "trimmed"))

  mfit <- bbqr(s$formula, s$data, quantile = c(0.3, 0.7),
               ndraw = 400, burn = 100)
  mm <- coef(mfit, statistic = "median")
  expect_equal(dim(mm), dim(coef(mfit)))
  expect_equal(mm[, 1L], coef(mfit[[1L]], statistic = "median"))
})

test_that("summary() returns the median last and prints what it always did", {
  s <- sim_binary(n = 200, seed = 28)
  fit <- bbqr(s$formula, s$data, ndraw = 600, burn = 100)
  tab <- summary(fit)$coefficients
  expect_named(tab, c("mean", "sd", "lower", "upper", "selected", "median"))
  expect_equal(tab$median, unname(coef(fit, statistic = "median")))
  out <- capture.output(print(summary(fit)))
  expect_true(any(grepl("Mean +SD +Lower +Upper", out)))
  expect_false(any(grepl("median", out, ignore.case = TRUE)))
})

test_that("predict() estimators do what they say", {
  s <- sim_binary(n = 120, seed = 29)
  tau <- 0.35
  fit <- bbqr(s$formula, s$data, quantile = tau, anchor = "free",
              ndraw = 500, burn = 100)
  X <- cbind(1, as.matrix(s$data[, c("x1", "x2", "x3")]))

  ## the default is the plug-in mean, for every type
  for (ty in c("response", "link", "class"))
    expect_identical(predict(fit, newdata = s$data, type = ty),
                     predict(fit, newdata = s$data, type = ty,
                             estimator = "plug_in_mean"), info = ty)

  ## plug-in median: medians of the coefficients and of sigma
  med <- coef(fit, statistic = "median")
  eta_med <- as.vector(X %*% med)
  expect_equal(predict(fit, newdata = s$data, type = "link",
                       estimator = "plug_in_median"), eta_med)
  expect_identical(predict(fit, newdata = s$data, type = "class",
                           estimator = "plug_in_median"),
                   as.integer(eta_med > 0))
  expect_equal(predict(fit, newdata = s$data[1:10, ],
                       estimator = "plug_in_median"),
               ald_tail_mass(eta_med[1:10], tau, stats::median(fit$sigma)),
               tolerance = 1e-7)

  ## draw average: each draw with its own coefficients and its own sigma
  avg <- predict(fit, newdata = s$data, estimator = "draw_average")
  expect_length(avg, nrow(s$data))
  expect_true(all(avg > 0 & avg < 1))
  by_hand <- vapply(1:6, function(i) {
    eta_m <- fit$beta0 + as.vector(fit$beta %*% X[i, -1])
    mean(mapply(function(e, sg) ald_tail_mass(e, tau, sg), eta_m, fit$sigma))
  }, numeric(1))
  expect_equal(avg[1:6], by_hand, tolerance = 1e-7)

  ## the answer does not depend on how the rows are blocked
  expect_equal(bbqr:::.bbqr_draw_average(fit, X, tau, cells = 1), avg)

  expect_error(predict(fit, newdata = s$data, type = "class",
                       estimator = "draw_average"), "plug_in_median")
  expect_error(predict(fit, newdata = s$data, type = "link",
                       estimator = "draw_average"), "plug_in_median")
})

test_that("multiple quantiles give a bbqr.list with one fit per level", {
  s <- sim_binary()
  fit <- bbqr(s$formula, s$data, quantile = c(0.25, 0.5, 0.75),
              ndraw = 500, burn = 100)
  expect_s3_class(fit, "bbqr.list")
  expect_length(fit, 3L)
  expect_equal(vapply(fit, function(z) z$quantile, numeric(1)),
               c(0.25, 0.5, 0.75), ignore_attr = TRUE)
  expect_equal(dim(coef(fit)), c(4L, 3L))
  expect_length(summary(fit), 3L)
})

test_that("print methods run without error", {
  s <- sim_binary()
  fit  <- bbqr(s$formula, s$data, ndraw = 400, burn = 100)
  mfit <- bbqr(s$formula, s$data, quantile = c(0.4, 0.6),
               ndraw = 400, burn = 100)
  expect_output(print(fit), "Bayesian binary quantile regression")
  expect_output(print(summary(fit)), "credible intervals")
  expect_output(print(mfit), "2 quantiles")
  expect_output(print(prior("alasso")), "alasso")
})

test_that("plot() runs for each type", {
  s <- sim_binary()
  fit <- bbqr(s$formula, s$data, ndraw = 400, burn = 100)
  pdf(NULL)
  on.exit(dev.off(), add = TRUE)
  expect_invisible(plot(fit, which = "x1"))
  expect_silent(plot(fit, which = c("x1", "x2"), type = "trace"))
  expect_silent(plot(fit, which = 1, type = "density"))
  expect_error(plot(fit, which = "nope"))
})