box::use(
  testthat[
    expect_equal,
    expect_error,
    expect_false,
    expect_named,
    expect_null,
    expect_true,
    skip_if_not_installed,
    test_that
  ],
  withr[local_options]
)

box::use(
  artma / econometric / bma[get_bma_data, run_bma],
  artma / econometric / fma[run_fma, unscale_fma_coefs],
  artma / methods / fma[fma, resolve_fma_cluster_ids]
)

make_demo_fma_data <- function() {
  set.seed(123)
  n <- 40L
  data.frame(
    effect = rnorm(n, mean = 0.2, sd = 0.1),
    se = runif(n, min = 0.05, max = 0.15),
    moderator1 = rnorm(n),
    moderator2 = rbinom(n, size = 1, prob = 0.5),
    stringsAsFactors = FALSE
  )
}

test_that("run_fma returns coefficients and weights", {
  skip_if_not_installed("BMS")
  skip_if_not_installed("quadprog")

  df <- make_demo_fma_data()

  var_list <- data.frame(
    var_name = c("effect", "se", "moderator1", "moderator2"),
    var_name_verbose = c("Effect", "SE", "Moderator 1", "Moderator 2"),
    bma = c(TRUE, TRUE, TRUE, TRUE),
    to_log_for_bma = c(FALSE, FALSE, FALSE, FALSE),
    bma_reference_var = c(FALSE, FALSE, FALSE, FALSE),
    stringsAsFactors = FALSE
  )

  bma_data <- get_bma_data(
    df,
    var_list,
    variable_info = c("effect", "se", "moderator1", "moderator2"),
    scale_data = FALSE,
    from_vector = TRUE,
    include_reference_groups = FALSE
  )

  params <- list(
    burn = 100L,
    iter = 500L,
    nmodel = 10L,
    g = "UIP",
    mprior = "uniform",
    mcmc = "bd"
  )

  bma_model <- run_bma(bma_data, params)

  result <- run_fma(
    bma_data = bma_data,
    bma_model = bma_model,
    input_var_list = var_list,
    print_results = "none"
  )

  expect_named(result, c("coefficients", "weights"))
  expect_true(is.data.frame(result$coefficients))
  expect_true(all(c("variable", "coefficient", "se", "p_value") %in% colnames(result$coefficients)))
  expect_true("Intercept" %in% result$coefficients$variable)

  expect_equal(length(result$weights), ncol(bma_data))
  expect_true(abs(sum(result$weights) - 1) < 1e-06)
  expect_true(all(result$weights >= -1e-08))
})

test_that("Mallows penalty favors small models on pure-noise predictors", {
  skip_if_not_installed("BMS")
  skip_if_not_installed("quadprog")

  set.seed(42)
  n <- 200L
  df <- data.frame(
    effect = rnorm(n),
    noise1 = rnorm(n),
    noise2 = rnorm(n),
    noise3 = rnorm(n),
    noise4 = rnorm(n),
    noise5 = rnorm(n),
    stringsAsFactors = FALSE
  )

  var_list <- data.frame(
    var_name = colnames(df),
    var_name_verbose = colnames(df),
    bma = rep(TRUE, ncol(df)),
    to_log_for_bma = rep(FALSE, ncol(df)),
    bma_reference_var = rep(FALSE, ncol(df)),
    stringsAsFactors = FALSE
  )

  bma_data <- get_bma_data(
    df,
    var_list,
    variable_info = colnames(df),
    scale_data = FALSE,
    from_vector = TRUE,
    include_reference_groups = FALSE
  )

  params <- list(
    burn = 100L, iter = 500L, nmodel = 10L,
    g = "UIP", mprior = "uniform", mcmc = "bd"
  )
  bma_model <- run_bma(bma_data, params)

  result <- run_fma(
    bma_data = bma_data,
    bma_model = bma_model,
    input_var_list = var_list,
    print_results = "none"
  )

  # None of the predictors carry signal, so the complexity penalty must pull
  # weight onto the smallest models. The sign-flipped criterion concentrated
  # nearly all weight on the largest model instead.
  m <- length(result$weights)
  small_half <- sum(result$weights[seq_len(ceiling(m / 2))])
  expect_true(small_half > 0.5)
  expect_true(result$weights[m] < 0.5)
})

# Study-clustered DGP: shared study shocks make cluster-robust and iid SEs
# diverge visibly. Every predictor carries signal so the full model keeps
# Mallows weight; a predictor excluded from all weighted models would get an
# exact-zero SE and trip the positivity checks below.
make_clustered_fma_inputs <- function() {
  set.seed(321)
  n_studies <- 12L
  per_study <- 5L
  n <- n_studies * per_study
  study <- rep(seq_len(n_studies), each = per_study)
  study_shock <- stats::rnorm(n_studies, sd = 0.1)

  se_col <- stats::runif(n, min = 0.05, max = 0.15)
  moderator1 <- stats::rnorm(n)
  moderator2 <- stats::rbinom(n, size = 1, prob = 0.5)
  df <- data.frame(
    effect = 0.2 + 2 * se_col + 0.8 * moderator1 + 0.6 * moderator2 +
      study_shock[study] + stats::rnorm(n, sd = 0.05),
    se = se_col,
    moderator1 = moderator1,
    moderator2 = moderator2,
    stringsAsFactors = FALSE
  )

  var_list <- data.frame(
    var_name = c("effect", "se", "moderator1", "moderator2"),
    var_name_verbose = c("effect", "se", "moderator1", "moderator2"),
    bma = rep(TRUE, 4L),
    to_log_for_bma = rep(FALSE, 4L),
    bma_reference_var = rep(FALSE, 4L),
    stringsAsFactors = FALSE
  )

  bma_data <- get_bma_data(
    df,
    var_list,
    variable_info = c("effect", "se", "moderator1", "moderator2"),
    scale_data = FALSE,
    from_vector = TRUE,
    include_reference_groups = FALSE
  )

  params <- list(
    burn = 100L, iter = 500L, nmodel = 10L,
    g = "UIP", mprior = "uniform", mcmc = "bd"
  )

  list(
    study = study,
    var_list = var_list,
    bma_data = bma_data,
    bma_model = run_bma(bma_data, params)
  )
}

test_that("cluster_ids changes standard errors but not coefficients or weights", {
  skip_if_not_installed("BMS")
  skip_if_not_installed("quadprog")

  inputs <- make_clustered_fma_inputs()

  iid <- run_fma(
    bma_data = inputs$bma_data,
    bma_model = inputs$bma_model,
    input_var_list = inputs$var_list,
    print_results = "none"
  )
  clustered <- run_fma(
    bma_data = inputs$bma_data,
    bma_model = inputs$bma_model,
    input_var_list = inputs$var_list,
    cluster_ids = inputs$study,
    print_results = "none"
  )

  # Clustering only reweights the score contributions inside the vcov; the
  # point estimates and the Mallows weights must be untouched.
  expect_equal(clustered$coefficients$coefficient, iid$coefficients$coefficient)
  expect_equal(clustered$weights, iid$weights)

  expect_true(all(is.finite(clustered$coefficients$se)))
  expect_true(all(clustered$coefficients$se > 0))
  expect_false(isTRUE(all.equal(clustered$coefficients$se, iid$coefficients$se)))
})

test_that("clustered FMA standard errors match a sandwich::vcovCL reference", {
  skip_if_not_installed("BMS")
  skip_if_not_installed("quadprog")

  inputs <- make_clustered_fma_inputs()

  result <- run_fma(
    bma_data = inputs$bma_data,
    bma_model = inputs$bma_model,
    input_var_list = inputs$var_list,
    cluster_ids = inputs$study,
    print_results = "none"
  )

  # var_name_verbose equals var_name here, so the output rows reveal the
  # predictor order run_fma actually used.
  predictor_order <- setdiff(result$coefficients$variable, "Intercept")
  x_data <- cbind(1, as.matrix(inputs$bma_data[predictor_order]))
  scale_vector <- apply(abs(x_data), 2, max)
  x_scaled <- sweep(x_data, 2, scale_vector, "/")
  y <- inputs$bma_data$effect

  m <- ncol(x_scaled)
  beta <- matrix(0, nrow = m, ncol = m)
  var_matrix <- matrix(0, nrow = m, ncol = m)
  for (i in seq_len(m)) {
    x_i <- x_scaled[, seq_len(i), drop = FALSE]
    fit_i <- stats::lm(y ~ 0 + x_i)
    beta[seq_len(i), i] <- stats::coef(fit_i)
    var_matrix[seq_len(i), i] <- diag(
      sandwich::vcovCL(fit_i, cluster = inputs$study, type = "HC1")
    )
  }

  weights <- result$weights
  beta_avg <- beta %*% weights
  bias_sq <- (beta - as.numeric(beta_avg))^2
  expected_coef <- as.numeric(beta_avg / scale_vector)
  expected_se <- as.numeric((sqrt(var_matrix + bias_sq) %*% weights) / scale_vector)

  expect_equal(result$coefficients$coefficient, expected_coef, tolerance = 1e-6)
  expect_equal(result$coefficients$se, expected_se, tolerance = 1e-6)
})

test_that("run_fma validates cluster_ids", {
  skip_if_not_installed("BMS")
  skip_if_not_installed("quadprog")

  inputs <- make_clustered_fma_inputs()
  n <- nrow(inputs$bma_data)

  run_with_cluster <- function(cluster_ids) {
    run_fma(
      bma_data = inputs$bma_data,
      bma_model = inputs$bma_model,
      input_var_list = inputs$var_list,
      cluster_ids = cluster_ids,
      print_results = "none"
    )
  }

  expect_error(run_with_cluster(rep("a", n)), "at least two distinct clusters")
  expect_error(run_with_cluster(inputs$study[-1]))
  expect_error(run_with_cluster(c(NA, inputs$study[-1])))
})

test_that("resolve_fma_cluster_ids aligns study_id through na.omit row names", {
  df <- data.frame(
    effect = c(0.1, NA, 0.3, 0.4, 0.5, 0.6),
    se = c(0.05, 0.06, 0.07, NA, 0.09, 0.1),
    study_id = c("s1", "s1", "s2", "s2", "s3", "s3"),
    stringsAsFactors = FALSE
  )
  bma_data <- stats::na.omit(df[c("effect", "se")])

  ids <- resolve_fma_cluster_ids(df, bma_data)
  expect_equal(ids, c("s1", "s2", "s3", "s3"))
  expect_equal(length(ids), nrow(bma_data))
})

test_that("resolve_fma_cluster_ids downgrades to NULL on unusable clusters", {
  bma_data <- data.frame(effect = c(0.1, 0.2), se = c(0.05, 0.06))

  no_study <- data.frame(effect = c(0.1, 0.2), se = c(0.05, 0.06))

  single_cluster <- cbind(no_study, study_id = c("s1", "s1"))

  na_ids <- cbind(no_study, study_id = c("s1", NA))

  # Rows that cannot be matched back to df resolve to NA and downgrade too.
  foreign_rows <- data.frame(effect = 0.1, se = 0.05, study_id = "s1", row.names = "99")

  # Each downgrade warns via cli_alert_warning(), which is the point of this
  # test; the test session's global sink (see setup.R) keeps it off the log.
  expect_null(resolve_fma_cluster_ids(no_study, bma_data))
  expect_null(resolve_fma_cluster_ids(single_cluster, bma_data))
  expect_null(resolve_fma_cluster_ids(na_ids, bma_data))
  expect_null(resolve_fma_cluster_ids(foreign_rows, bma_data))
})

test_that("unscale_fma_coefs inverts the z-scoring of slopes and intercept", {
  bma_data <- data.frame(effect = 1:4, x_cont = c(1, 3, 5, 9), x_bin = c(0, 1, 0, 1))
  attr(bma_data, "bpe_scale_centers") <- c(effect = 2.5, x_cont = 4.5, x_bin = 0)
  attr(bma_data, "bpe_scale_scales") <- c(effect = 2, x_cont = 4, x_bin = 1)

  var_names <- c("(Intercept)", "x_cont", "x_bin")
  coefficient <- c(0.3, 0.5, -0.25)
  se <- c(0.2, 0.1, 0.05)

  out <- unscale_fma_coefs(coefficient, se, var_names, bma_data)

  # Slopes: b * sd(y) / sd(x); the dummy only picks up sd(y).
  expect_equal(out$coefficient[2], 0.5 * 2 / 4)
  expect_equal(out$se[2], 0.1 * 2 / 4)
  expect_equal(out$coefficient[3], -0.25 * 2)
  expect_equal(out$se[3], 0.05 * 2)
  # Intercept: mean(y) + sd(y) * a - sum(raw slope * mean(x)).
  expect_equal(out$coefficient[1], 2.5 + 2 * 0.3 - 0.25 * 4.5)
  expect_equal(out$se[1], 2 * 0.2)

  # A frame that was never scaled is passed through unchanged.
  plain <- unscale_fma_coefs(coefficient, se, var_names, data.frame(effect = 1:4))
  expect_equal(plain$coefficient, coefficient)
  expect_equal(plain$se, se)
})

test_that("run_fma on the z-scored frame matches a fit on the raw data", {
  skip_if_not_installed("BMS")
  skip_if_not_installed("quadprog")

  df <- make_demo_fma_data()
  vars <- c("effect", "se", "moderator1", "moderator2")
  var_list <- data.frame(
    var_name = vars,
    var_name_verbose = vars,
    bma = rep(TRUE, 4L),
    to_log_for_bma = rep(FALSE, 4L),
    bma_reference_var = rep(FALSE, 4L),
    stringsAsFactors = FALSE
  )

  raw <- get_bma_data(df, var_list, vars, scale_data = FALSE, from_vector = TRUE, include_reference_groups = FALSE)
  scaled <- get_bma_data(df, var_list, vars, scale_data = TRUE, from_vector = TRUE, include_reference_groups = FALSE)

  # One model for both runs: it only supplies the predictor ordering, so the
  # two FMA fits average over the same nested sequence.
  bma_model <- run_bma(raw, list(
    burn = 100L, iter = 500L, nmodel = 10L, g = "UIP", mprior = "uniform", mcmc = "bd"
  ))

  from_raw <- run_fma(raw, bma_model, var_list, print_results = "none")
  from_scaled <- run_fma(scaled, bma_model, var_list, print_results = "none")

  expect_equal(from_scaled$coefficients$variable, from_raw$coefficients$variable)
  expect_equal(from_scaled$weights, from_raw$weights, tolerance = 1e-6)
  # The back-transform is exact: the averaging weights are shared and the
  # rescaling is linear, so both the slopes and the corrected intercept land
  # on the raw-data numbers.
  expect_equal(from_scaled$coefficients$coefficient, from_raw$coefficients$coefficient, tolerance = 1e-6)
  # Slope SEs are invariant to centering the regressors; the intercept SE is a
  # plain rescale that carries no correction for the shift, as in BMA.
  slopes <- from_raw$coefficients$variable != "Intercept"
  expect_equal(from_scaled$coefficients$se[slopes], from_raw$coefficients$se[slopes], tolerance = 1e-6)
  # p-values are scale invariant for the slopes.
  expect_equal(from_scaled$coefficients$p_value[slopes], from_raw$coefficients$p_value[slopes], tolerance = 1e-6)

  # The test only means something if the fit on the z-scored frame really did
  # land somewhere else: dropping the attributes reproduces the pre-fix
  # y-SD-per-x-SD reporting.
  bare <- scaled
  attr(bare, "bpe_scale_centers") <- NULL
  attr(bare, "bpe_scale_scales") <- NULL
  standardized <- run_fma(bare, bma_model, var_list, print_results = "none")
  expect_false(isTRUE(all.equal(
    standardized$coefficients$coefficient, from_raw$coefficients$coefficient
  )))
})

test_that("fma reports coefficients on the data scale", {
  skip_if_not_installed("BMS")
  skip_if_not_installed("quadprog")

  df <- make_demo_fma_data()
  vars <- c("effect", "se", "moderator1", "moderator2")
  var_list <- data.frame(
    var_name = vars,
    var_name_verbose = vars,
    bma = rep(TRUE, 4L),
    to_log_for_bma = rep(FALSE, 4L),
    bma_reference_var = rep(FALSE, 4L),
    stringsAsFactors = FALSE
  )

  local_options(list(
    artma.verbose = 0,
    artma.autonomy.level = "autonomous",
    artma.data.columns = list(
      effect = list(var_name = "effect", var_name_verbose = "effect", bma = FALSE),
      se = list(var_name = "se", var_name_verbose = "se", bma = TRUE),
      moderator1 = list(var_name = "moderator1", var_name_verbose = "moderator1", bma = TRUE),
      moderator2 = list(var_name = "moderator2", var_name_verbose = "moderator2", bma = TRUE)
    ),
    artma.output.save_results = FALSE,
    artma.visualization.export_graphics = FALSE,
    artma.methods.bma.burn = 100L,
    artma.methods.bma.iter = 500L,
    artma.methods.bma.nmodel = 10L,
    artma.methods.bma.g = "UIP",
    artma.methods.bma.mprior = "uniform",
    artma.methods.bma.mcmc = "bd",
    artma.methods.fma.cluster = FALSE
  ))

  result <- fma(df)

  raw <- get_bma_data(df, var_list, vars, scale_data = FALSE, from_vector = TRUE, include_reference_groups = FALSE)
  reference <- run_fma(raw, result$meta$model, var_list, print_results = "none")

  # The models still ran on the z-scored frame ...
  expect_true(abs(mean(result$meta$data$effect)) < 1e-10)
  # ... but the reported numbers are on the data scale, intercept included.
  expected <- stats::setNames(reference$coefficients$coefficient, reference$coefficients$variable)
  reported <- stats::setNames(result$estimates$estimate, result$estimates$term)
  expect_equal(reported[names(expected)], expected, tolerance = 1e-6)
  expect_equal(result$tables$coefficients$coefficient, result$estimates$estimate)
})