tests/testthat/test-integration.R

## Integration tests for coresynth
## Run: devtools::test()

library(testthat)

# ── Shared helpers ────────────────────────────────────────────────────────────

# Balanced panel: 1 treated unit (u1), rest control, treatment from T_pre+1 onward.
# DGP: Y_it = alpha_i + f_t * lambda_i + tau*D_it + eps
make_panel <- function(N = 10, T = 20, T_pre = 10, effect = 2.0, seed = 42) {
  set.seed(seed)
  times <- 1:T
  units <- paste0("u", 1:N)
  f <- cumsum(rnorm(T, 0, 0.5)) # common factor
  lam <- rnorm(N, 1, 0.3) # unit loadings
  rows <- expand.grid(time = times, id = units, stringsAsFactors = FALSE)
  rows$y <- as.vector(outer(f, lam)) + rnorm(nrow(rows), 0, 0.3)
  rows$d <- as.integer(rows$id == "u1" & rows$time > T_pre)
  rows$y[rows$d == 1] <- rows$y[rows$d == 1] + effect
  rows[order(rows$id, rows$time), ]
}

panel <- make_panel()

# ── SCM ──────────────────────────────────────────────────────────────────────
test_that("SCM weights lie on the unit simplex", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_s3_class(fit, "coresynth")
  expect_equal(fit$method, "scm")
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
  expect_equal(length(fit$Y_synth), 20)
  expect_equal(length(fit$gap), 20)
})

test_that("SCM estimate is in a reasonable range for effect=2", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_true(abs(fit$estimate - 2.0) < 3.0)
})

test_that("SCM pre-treatment fit is good (low MSPE)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  pre <- seq_len(10)
  mspe <- mean((fit$Y_treat[pre] - fit$Y_synth[pre])^2)
  total_var <- var(fit$Y_treat[pre])
  expect_true(mspe / total_var < 0.5) # synthetic tracks pre-period reasonably
})

# ── SDID ─────────────────────────────────────────────────────────────────────
test_that("SDID unit weights lie on the unit simplex", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_s3_class(fit, "coresynth")
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
})

test_that("SDID time weights lie on the unit simplex", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_equal(sum(fit$time_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$time_weights >= -1e-5))
  expect_equal(length(fit$time_weights), 10L) # T_pre = 10
})

test_that("SDID zeta2 follows Arkhangelsky et al. (2021) Eq.(2.2)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  # zeta2 = sqrt(N_tr * T_post) * sigma2_hat
  # N_tr=1, T_post=10, sigma2_hat from fit
  expected_zeta2 <- sqrt(1 * 10) * fit$sigma2_hat
  expect_equal(fit$zeta2, expected_zeta2, tolerance = 1e-10)
  expect_true(fit$zeta2 > 0)
})

test_that("SDID omega_0 and lambda_0 intercepts are finite", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_true(is.finite(fit$omega0))
  expect_true(is.finite(fit$lambda0))
})

test_that("SDID demeaning shifts unit weights vs. no-intercept baseline", {
  # With intercept (current): demeaned QP allows level shifts
  fit_new <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  # The estimate should still be in a reasonable range
  expect_true(abs(fit_new$estimate - 2.0) < 4.0)
})

test_that("SDID estimate matches synthdid package (California smoke data)", {
  skip_if_not_installed("synthdid")
  data("california_prop99", package = "synthdid")
  setup <- synthdid::panel.matrices(california_prop99)
  ref <- synthdid::synthdid_estimate(setup$Y, setup$N0, setup$T0)
  # Convert data to long format for coresynth
  Y_wide <- setup$Y
  df <- data.frame(
    y = as.vector(Y_wide),
    id = rep(rownames(Y_wide), ncol(Y_wide)),
    time = rep(colnames(Y_wide), each = nrow(Y_wide))
  )
  N0 <- setup$N0
  T0 <- setup$T0
  post_cols <- (T0 + 1):ncol(Y_wide)
  treat_rows <- (N0 + 1):nrow(Y_wide)
  df$d <- as.integer(
    df$id %in%
      rownames(Y_wide)[treat_rows] &
      df$time %in% colnames(Y_wide)[post_cols]
  )
  fit_cs <- scm_fit(y ~ d | id + time, data = df, method = "sdid")
  # Reference: Arkhangelsky et al. (2021) Table 1 reports SDID = -15.6
  expect_equal(fit_cs$estimate, as.numeric(ref), tolerance = 0.5)
})

# ── GSC ──────────────────────────────────────────────────────────────────────
test_that("GSC returns factor matrix with correct dimensions", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  expect_s3_class(fit, "coresynth")
  expect_equal(ncol(fit$F), 2L) # r = 2 factors
  expect_equal(nrow(fit$F), 20L) # T periods
  expect_equal(ncol(fit$L_co), 2L) # r loadings per control unit
  expect_equal(nrow(fit$L_co), 9L) # N_co = 9
})

test_that("GSC counterfactual has correct dimension", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  expect_equal(dim(fit$Y_tr_hat), c(20L, 1L)) # T x N_tr
})

test_that("GSC estimate is finite and in reasonable range", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  expect_true(is.finite(fit$estimate))
  expect_true(abs(fit$estimate - 2.0) < 4.0)
})

test_that("GSC singular values are decreasing", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  sv <- fit$singular_values
  expect_true(all(diff(sv[sv > 0]) <= 1e-10)) # non-increasing
})

test_that("GSC rejects r larger than N_co or T", {
  # R wrapper fires first (r > N_co = 9): "Need at least r control units"
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 20),
    regexp = "at least r"
  )
})

# ── MC ───────────────────────────────────────────────────────────────────────
test_that("MC returns low-rank matrix of correct size", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  expect_s3_class(fit, "coresynth")
  expect_equal(dim(fit$L_hat), c(20L, 10L))
})

test_that("MC observed entries are well-fitted (low residual on mask)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  pan <- coresynth:::panel_to_matrices(
    panel$y,
    panel$d,
    panel$id,
    panel$time
  )
  # Observation mask: post x treated = 0, rest = 1
  O <- matrix(1, nrow = 20, ncol = 10)
  O[11:20, 1] <- 0 # unit u1 post-treatment is "missing"
  Y_obs <- pan$Y
  Y_obs[is.na(Y_obs)] <- 0
  resid_obs <- (Y_obs - fit$L_hat)[O == 1]
  rmse_obs <- sqrt(mean(resid_obs^2))
  # Should be small relative to outcome range
  expect_true(rmse_obs < diff(range(Y_obs[O == 1])) * 0.3)
})

test_that("MC estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  expect_true(is.finite(fit$estimate))
})

# ── TASC ─────────────────────────────────────────────────────────────────────
test_that("TASC runs and returns valid output structure", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel,
    method = "tasc",
    r = 2,
    em_iter = 5
  )
  expect_s3_class(fit, "coresynth")
  expect_true(is.finite(fit$estimate))
  expect_equal(length(fit$times), 20L)
  expect_equal(dim(fit$Y_hat), c(20L, 10L))
})

test_that("TASC general A: ATT is finite", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel,
    method = "tasc",
    r = 2,
    em_iter = 10
  )
  expect_true(is.finite(fit$estimate))
})

test_that("TASC A matrix has correct dimensions (r x r)", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel,
    method = "tasc",
    r = 2,
    em_iter = 5
  )
  expect_equal(dim(fit$A), c(2L, 2L))
})

test_that("TASC A deviates from identity after EM (learning occurred)", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel,
    method = "tasc",
    r = 2,
    em_iter = 20
  )
  expect_true(max(abs(fit$A - diag(2))) > 1e-6)
})

test_that("TASC fix_A=TRUE keeps A as identity", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel,
    method = "tasc",
    r = 2,
    em_iter = 10,
    fix_A = TRUE
  )
  expect_equal(fit$A, diag(2))
  expect_true(is.finite(fit$estimate))
})

test_that("TASC fix_A=FALSE and fix_A=TRUE both return finite ATT", {
  fit_free <- scm_fit(
    y ~ d | id + time,
    data = panel,
    method = "tasc",
    r = 2,
    em_iter = 10,
    fix_A = FALSE
  )
  fit_fix <- scm_fit(
    y ~ d | id + time,
    data = panel,
    method = "tasc",
    r = 2,
    em_iter = 10,
    fix_A = TRUE
  )
  expect_true(is.finite(fit_free$estimate))
  expect_true(is.finite(fit_fix$estimate))
})

# ── print / summary ──────────────────────────────────────────────────────────
test_that("print and summary work for all methods", {
  for (m in c("scm", "sdid", "gsc", "mc")) {
    fit <- scm_fit(y ~ d | id + time, data = panel, method = m)
    expect_output(print(fit), toupper(m))
    expect_output(summary(fit), "ATT estimate")
  }
})

# ── broom tidiers ─────────────────────────────────────────────────────────────
test_that("tidy() and glance() work for SDID", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  td <- broom::tidy(fit)
  gl <- broom::glance(fit)
  expect_s3_class(td, "data.frame")
  expect_true(nrow(td) > 0)
  expect_true(all(c("term", "estimate") %in% names(td)))
  expect_equal(gl$method, "sdid")
})

# ── plot ─────────────────────────────────────────────────────────────────────
test_that("plot.coresynth returns ggplot objects", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_s3_class(plot(fit, type = "trend"), "ggplot")
  expect_s3_class(plot(fit, type = "gap"), "ggplot")
  expect_s3_class(plot(fit, type = "weights"), "ggplot")
})

# ── SI ───────────────────────────────────────────────────────────────────────

test_that("SI returns correct class and structure", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  expect_s3_class(fit, "coresynth")
  expect_s3_class(fit, "coresynth_si")
  expect_equal(fit$method, "si")
  expect_true(is.numeric(fit$k) && fit$k >= 1L)
})

test_that("SI output dimensions are correct", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  # panel: T=20, N=10, N_co=9, N_tr=1, T_pre=10, T_post=10
  expect_equal(dim(fit$Y_treat), c(20L, 1L))
  expect_equal(dim(fit$Y_synth), c(20L, 1L))
  expect_equal(dim(fit$Y_cf), c(10L, 1L))
  expect_equal(dim(fit$weights), c(9L, 1L))
  expect_equal(length(fit$unit_weights), 9L)
  expect_equal(length(fit$times), 20L)
  expect_equal(fit$T_pre, 10L)
})

test_that("SI pre-treatment fit is good (low MSPE)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  pre <- seq_len(10)
  mspe <- mean((fit$Y_treat[pre, ] - fit$Y_synth[pre, ])^2)
  total_var <- var(as.vector(fit$Y_treat[pre, ]))
  expect_true(mspe / total_var < 0.5)
})

test_that("SI estimate is finite and in reasonable range for effect=2", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  expect_true(is.finite(fit$estimate))
  expect_true(abs(fit$estimate - 2.0) < 3.0)
})

test_that("SI respects user-supplied k", {
  fit1 <- scm_fit(y ~ d | id + time, data = panel, method = "si", k = 1L)
  fit3 <- scm_fit(y ~ d | id + time, data = panel, method = "si", k = 3L)
  expect_equal(fit1$k, 1L)
  expect_equal(fit3$k, 3L)
  expect_true(is.finite(fit1$estimate))
  expect_true(is.finite(fit3$estimate))
})

test_that("SI errors when k > min(T_pre, N_co)", {
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "si", k = 100L),
    regexp = "k"
  )
})

test_that("tensor_unfold_cpp produces correct dimensions for all modes", {
  tc <- array(seq_len(3 * 4 * 2), dim = c(3L, 4L, 2L))
  expect_equal(dim(tensor_unfold_cpp(tc, 1L)), c(3L, 8L))
  expect_equal(dim(tensor_unfold_cpp(tc, 2L)), c(4L, 6L))
  expect_equal(dim(tensor_unfold_cpp(tc, 3L)), c(2L, 12L))
  expect_error(tensor_unfold_cpp(tc, 4L), regexp = "Mode")
})

test_that("SI and MC estimates are in similar range", {
  fit_si <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  fit_mc <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  expect_true(abs(fit_si$estimate - fit_mc$estimate) < 3.0)
})

test_that("print and summary work for SI", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  expect_output(print(fit), "SI")
  expect_output(summary(fit), "ATT estimate")
})

# ── edge cases ───────────────────────────────────────────────────────────────
test_that("scm_fit errors on bad formula", {
  expect_error(scm_fit(y ~ d, data = panel, method = "sdid"), "after '|'")
})

test_that("scm_fit errors on missing variable", {
  expect_error(
    scm_fit(z ~ d | id + time, data = panel, method = "sdid"),
    "not found"
  )
})

# ── Phase 4: SCM placebo / MSPE-ratio inference ───────────────────────────────

test_that("fit_scm_cpp returns Y_co_pre and Y_co_post", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_equal(dim(fit$Y_co_pre), c(10L, 9L)) # T_pre x N_co
  expect_equal(dim(fit$Y_co_post), c(10L, 9L)) # T_post x N_co
})

test_that("scm_placebo_cpp returns correct structure", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  plac <- scm_placebo_cpp(fit$Y_co_pre, fit$Y_co_post)
  expect_equal(length(plac$mspe_pre), 9L)
  expect_equal(length(plac$mspe_post), 9L)
  expect_equal(length(plac$effects), 9L)
  expect_true(all(plac$mspe_pre >= 0))
  expect_true(all(plac$mspe_post >= 0))
  expect_true(all(is.finite(plac$effects)))
})

test_that("mspe_ratio_pval returns p_value in [0, 1]", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  expect_true(is.finite(inf$mspe_ratio_treated))
  expect_true(inf$mspe_ratio_treated >= 0)
  expect_equal(length(inf$mspe_ratios_all), 10L) # 1 treated + 9 controls
  expect_equal(length(inf$placebo_effects), 9L)
})

test_that("mspe_ratio_pval mspe_threshold filters all controls", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit, mspe_threshold = 1e10)
  expect_equal(inf$n_placebo_used, 0L)
})

test_that("mspe_ratio_pval rejects non-SCM fits", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_error(mspe_ratio_pval(fit), "method = 'scm'")
})

test_that("mspe_ratio_pval default mirrors the fit's predictor spec", {
  fit_cov <- scm_fit(
    y ~ d | id + time,
    data = panel, method = "scm",
    predictors = list(pred("y", 1:10, op = "mean"))
  )
  auto     <- mspe_ratio_pval(fit_cov)
  explicit <- mspe_ratio_pval(fit_cov, use_covariates = TRUE)
  expect_equal(auto$p_value, explicit$p_value)
  expect_equal(auto$placebo_effects, explicit$placebo_effects)
  expect_equal(auto$mspe_pre_placebo, explicit$mspe_pre_placebo)

  fit_out       <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  auto_out      <- mspe_ratio_pval(fit_out)
  explicit_out  <- mspe_ratio_pval(fit_out, use_covariates = FALSE)
  expect_equal(auto_out$p_value, explicit_out$p_value)
  expect_equal(auto_out$placebo_effects, explicit_out$placebo_effects)
})

test_that("explicit use_covariates=TRUE without covariates warns and falls back", {
  fit_out <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_warning(
    inf <- mspe_ratio_pval(fit_out, use_covariates = TRUE),
    "X0_mat is NULL"
  )
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
})

test_that("outcomes-only placebo solver agrees with scm_weights_cpp per unit", {
  set.seed(7)
  T_pre <- 40L; T_post <- 8L; N <- 20L
  f <- cumsum(rnorm(T_pre + T_post, 0, 0.5))
  Y <- sapply(seq_len(N), function(j) {
    100 + rnorm(1, 0, 5) + f * rnorm(1, 1, 0.3) +
      rnorm(T_pre + T_post, 0, 0.5)
  })
  Y_pre  <- Y[seq_len(T_pre), ]
  Y_post <- Y[(T_pre + 1):(T_pre + T_post), ]

  plac <- scm_placebo_cpp(Y_pre, Y_post)
  for (i in c(1L, 10L, 20L)) {
    ref <- scm_weights_cpp(
      Y_pre[, -i, drop = FALSE], Y_pre[, i],
      Y_pre[, -i, drop = FALSE], Y_pre[, i]
    )
    mspe_ref <- mean((Y_pre[, i] - Y_pre[, -i, drop = FALSE] %*% ref$W)^2)
    # Same nested V/W estimator, different inner QP solver: results agree
    # within the coordinate-descent convergence tolerance.
    expect_lt(
      abs(plac$mspe_pre[i] - mspe_ref) / max(mspe_ref, 1e-12),
      5e-3
    )
  }
})

test_that("realistic-size outcomes-only placebo completes quickly", {
  skip_on_cran()
  set.seed(2)
  T_pre <- 139L; T_post <- 19L; N <- 75L  # ~12y monthly, 75 donors
  tt <- seq_len(T_pre + T_post)
  common <- 200 + 10 * sin(tt / 12) + 0.05 * tt
  Y <- sapply(seq_len(N), function(j) {
    common + rnorm(1, 0, 15) +
      as.numeric(arima.sim(list(ar = 0.8), length(tt), sd = 3))
  })
  t0 <- Sys.time()
  plac <- scm_placebo_cpp(
    Y[seq_len(T_pre), ],
    Y[(T_pre + 1):(T_pre + T_post), ]
  )
  elapsed <- as.numeric(difftime(Sys.time(), t0, units = "secs"))
  expect_true(all(is.finite(plac$mspe_pre)))
  expect_true(all(is.finite(plac$mspe_post)))
  expect_true(all(is.finite(plac$effects)))
  expect_equal(dim(plac$gaps), c(T_pre + T_post, N))
  # Generous bound (runs in ~2 s on a laptop; was ~an hour before the
  # active-set inner solver).
  expect_lt(elapsed, 60)
})

# ── Phase 4: GSC bootstrap inference ─────────────────────────────────────────

test_that("fit_gsc_cpp returns Y_co_all", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  expect_equal(dim(fit$Y_co_all), c(20L, 9L)) # T x N_co
})

test_that("gsc_boot returns correct structure", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  boot <- gsc_boot(fit, B = 29L, seed = 1L)
  expect_named(
    boot,
    c("p_value", "ci_lower", "ci_upper", "se", "boot_dist", "att_obs")
  )
  expect_equal(length(boot$boot_dist), 29L)
  expect_true(boot$p_value >= 0 && boot$p_value <= 1)
  expect_true(boot$ci_lower <= boot$ci_upper)
  expect_true(is.finite(boot$se) && boot$se >= 0)
  expect_equal(boot$att_obs, fit$estimate)
})

test_that("gsc_boot seed ensures reproducibility", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  boot1 <- gsc_boot(fit, B = 29L, seed = 42L)
  boot2 <- gsc_boot(fit, B = 29L, seed = 42L)
  expect_equal(boot1$boot_dist, boot2$boot_dist)
})

test_that("gsc_boot rejects non-GSC fits", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(gsc_boot(fit), "method = 'gsc'")
})

# ── Phase 4: export_json ──────────────────────────────────────────────────────

test_that("export_json NULL file returns list without writing", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  res <- export_json(fit, file = NULL)
  expect_type(res, "list")
  expect_equal(res$meta$method, "sdid")
  expect_equal(res$meta$package, "coresynth")
})

test_that("export_json has required top-level sections", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  res <- export_json(fit, file = NULL)
  expect_true(all(c("meta", "data", "estimate", "time_series") %in% names(res)))
})

test_that("export_json estimate section contains att", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  res <- export_json(fit, file = NULL)
  expect_equal(res$estimate$att, signif(fit$estimate, 6))
})

test_that("export_json time_series has required fields", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  res <- export_json(fit, file = NULL)
  expect_true(all(
    c("times", "Y_treat", "Y_synth", "gap") %in% names(res$time_series)
  ))
  expect_equal(length(res$time_series$times), 20L)
})

test_that("export_json writes valid JSON for SCM, SDID, GSC, MC", {
  skip_if_not_installed("jsonlite")
  tmp <- tempfile(fileext = ".json")
  on.exit(unlink(tmp), add = TRUE)
  for (m in c("scm", "sdid", "gsc", "mc")) {
    fit <- scm_fit(y ~ d | id + time, data = panel, method = m)
    export_json(fit, file = tmp)
    dat <- jsonlite::read_json(tmp)
    expect_equal(dat$meta$method, m)
    expect_false(is.null(dat$estimate$att))
    expect_equal(dat$data$T_pre, 10L)
  }
})

test_that("export_json with inference includes inference section (SCM)", {
  skip_if_not_installed("jsonlite")
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  res <- export_json(fit, file = NULL, inference = inf)
  expect_equal(res$inference$type, "mspe_permutation")
  expect_false(is.null(res$inference$p_value))
  expect_false(is.null(res$inference$mspe_ratios_all))
  expect_equal(res$estimate$p_value, signif(inf$p_value, 6))
})

test_that("export_json with inference includes inference section (GSC)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  boot <- gsc_boot(fit, B = 29L, seed = 1L)
  res <- export_json(fit, file = NULL, inference = boot)
  expect_equal(res$inference$type, "parametric_bootstrap")
  expect_equal(res$inference$B, 29L)
  expect_false(is.null(res$inference$ci_lower))
})

# ── Phase 5a: Staggered adoption ──────────────────────────────────────────────

make_staggered_panel <- function(N = 12, TT = 24, seed = 7) {
  set.seed(seed)
  f <- cumsum(rnorm(TT, 0, 0.5))
  lam <- rnorm(N, 1, 0.3)
  rows <- expand.grid(
    time = seq_len(TT),
    id = paste0("u", seq_len(N)),
    stringsAsFactors = FALSE
  )
  rows$y <- as.vector(outer(f, lam)) + rnorm(nrow(rows), 0, 0.3)
  # u1 採用 t=9、u2 採用 t=15
  rows$d <- as.integer(
    (rows$id == "u1" & rows$time >= 9) |
      (rows$id == "u2" & rows$time >= 15)
  )
  rows$y[rows$d == 1] <- rows$y[rows$d == 1] + 1.5
  rows[order(rows$id, rows$time), ]
}

staggered <- make_staggered_panel()

test_that("panel_to_matrices detects staggered adoption correctly", {
  pan <- coresynth:::panel_to_matrices(
    staggered$y,
    staggered$d,
    staggered$id,
    staggered$time
  )
  expect_false(pan$is_sharp)
  expect_equal(pan$T_pre, 8L)
  u1c <- which(pan$units == "u1")
  u2c <- which(pan$units == "u2")
  expect_equal(pan$T_adopt[u1c], 9L)
  expect_equal(pan$T_adopt[u2c], 15L)
})

test_that("panel_to_matrices is_sharp TRUE for non-staggered panel", {
  pan <- coresynth:::panel_to_matrices(
    panel$y,
    panel$d,
    panel$id,
    panel$time
  )
  expect_true(pan$is_sharp)
})

test_that("MC handles staggered adoption and returns finite estimate", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "mc")
  expect_true(is.finite(fit$estimate))
})

test_that("MC staggered: L_hat dimensions are correct", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "mc")
  expect_equal(dim(fit$L_hat), c(24L, 12L))
})

test_that("TASC handles staggered adoption and returns finite estimate", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = staggered,
    method = "tasc",
    r = 2L,
    em_iter = 5L
  )
  expect_true(is.finite(fit$estimate))
})

test_that("SCM staggered with predictors arg errors informatively", {
  expect_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm",
            predictors = list(pred("y", 1:3))),
    regexp = "predictors"
  )
})

test_that("GSC handles staggered adoption (Phase 16)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_true(isTRUE(fit$staggered))
  expect_true(is.finite(fit$estimate))
})

test_that("SI handles staggered adoption (Phase 16)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expect_true(isTRUE(fit$staggered))
  expect_true(is.finite(fit$estimate))
})

# ── SCM with covariates (pred() interface) ────────────────────────────────────

# Helper: panel with two extra covariate columns
make_panel_with_cov <- function(
  N = 10,
  T = 20,
  T_pre = 10,
  effect = 2.0,
  seed = 42
) {
  set.seed(seed)
  p <- make_panel(N = N, T = T, T_pre = T_pre, effect = effect, seed = seed)
  units_all <- unique(p$id)
  cov1_vals <- setNames(rnorm(N, mean = 5, sd = 1), units_all)
  cov2_vals <- setNames(rnorm(N, mean = 10, sd = 2), units_all)
  p$cov1 <- cov1_vals[p$id]
  p$cov2 <- cov2_vals[p$id]
  p
}

panel_cov <- make_panel_with_cov()

test_that("SCM with pred() returns valid fit object", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  expect_s3_class(fit, "coresynth")
  expect_equal(fit$method, "scm")
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
  expect_equal(length(fit$Y_synth), 20)
})

test_that("SCM with pred(): v_weights have predictor names", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  expect_named(fit$v_weights, c("cov1[1:10]", "cov2[1:10]"))
})

test_that("SCM with pred(): predictor_table is a data.frame with correct cols", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  expect_s3_class(fit$predictor_table, "data.frame")
  expect_true(all(
    c("predictor", "treated", "synthetic") %in% names(fit$predictor_table)
  ))
  expect_equal(nrow(fit$predictor_table), 2L)
})

test_that("SCM with multiple pred() returns valid fit object", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(
      pred("y", c(3, 5, 7)),
      pred("cov1", 1:10)
    )
  )
  expect_s3_class(fit, "coresynth")
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_equal(length(fit$v_weights), 2L)
  expect_true(all(is.finite(fit$Y_synth)))
})

test_that("SCM with mixed pred() entries stacks rows", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(
      pred("cov1", 1:10),
      pred("y", c(5, 8, 10))
    )
  )
  expect_equal(length(fit$v_weights), 2L)
  expect_equal(nrow(fit$predictor_table), 2L)
})

test_that("SCM without predictors matches previous outcome-only behaviour", {
  fit_plain <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  fit_cov <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm")
  expect_equal(fit_plain$unit_weights, fit_cov$unit_weights, tolerance = 1e-4)
  expect_null(fit_cov$predictor_table)
  expect_equal(
    names(fit_plain$v_weights),
    paste0("V", seq_along(fit_plain$v_weights))
  )
})

test_that("SCM with pred(): predictor_table synthetic matches W dot X0", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  # Synthetic predictor values should be a weighted combination of control values
  expect_true(all(is.finite(fit$predictor_table$synthetic)))
})

test_that("SCM errors when pred() variable not in data", {
  expect_error(
    scm_fit(
      y ~ d | id + time,
      data = panel_cov,
      method = "scm",
      predictors = list(pred("nonexistent_var", 1:5))
    ),
    regexp = "not found in data"
  )
})

test_that("SCM errors when predictors is a character vector (old API)", {
  expect_error(
    scm_fit(
      y ~ d | id + time,
      data = panel_cov,
      method = "scm",
      predictors = "cov1"
    ),
    regexp = "list of pred\\(\\) objects"
  )
})

# ── Phase 7a: GSC covariates (time-varying, EM loop) ─────────────────────────

# Helper: panel with one time-varying covariate column
make_panel_with_timecov <- function(
  N = 10,
  T = 20,
  T_pre = 10,
  effect = 2.0,
  seed = 42
) {
  p <- make_panel(N = N, T = T, T_pre = T_pre, effect = effect, seed = seed)
  set.seed(seed + 1L)
  units_all <- unique(p$id)
  for (u in units_all) {
    base <- rnorm(1, 5, 1)
    p$cov1[p$id == u] <- base +
      seq(0, by = 0.1, length.out = T) +
      rnorm(T, 0, 0.2)
  }
  p
}

panel_timecov <- make_panel_with_timecov()

test_that("GSC with covariates: fit$beta has length p", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_timecov,
    method = "gsc",
    r = 2,
    covariates = "cov1"
  )
  expect_equal(length(fit$beta), 1L)
  expect_true(is.numeric(fit$beta))
})

test_that("GSC with covariates: ATT is finite and in reasonable range", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_timecov,
    method = "gsc",
    r = 2,
    covariates = "cov1"
  )
  expect_true(is.finite(fit$estimate))
  expect_true(abs(fit$estimate - 2.0) < 5.0)
})

test_that("GSC with covariates: Y_hat_co has correct dimensions (T x N_co)", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_timecov,
    method = "gsc",
    r = 2,
    covariates = "cov1"
  )
  expect_equal(dim(fit$Y_hat_co), c(20L, 9L))
})

test_that("GSC with covariates via scm_fit formula interface", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_timecov,
    method = "gsc",
    r = 2,
    covariates = "cov1"
  )
  expect_s3_class(fit, "coresynth")
  expect_equal(fit$method, "gsc")
})

test_that("GSC without covariates backward compat: beta is empty, covariates NULL", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2)
  expect_equal(length(fit$beta), 0L)
  expect_null(fit$covariates)
  expect_true(is.finite(fit$estimate))
})

test_that("gsc_boot with covariates: p_value in [0,1], CI valid", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_timecov,
    method = "gsc",
    r = 2,
    covariates = "cov1"
  )
  boot <- gsc_boot(fit, B = 19L, seed = 1L)
  expect_true(boot$p_value >= 0 && boot$p_value <= 1)
  expect_true(boot$ci_lower <= boot$ci_upper)
  expect_equal(length(boot$boot_dist), 19L)
})

test_that("build_covariate_array returns correct dimensions (T x N x p)", {
  arr <- coresynth:::build_covariate_array(
    panel_timecov,
    "id",
    "time",
    "cov1",
    units = unique(panel_timecov$id),
    times = 1:20
  )
  expect_equal(dim(arr), c(20L, 10L, 1L))
  expect_true(all(is.finite(arr)))
})

test_that("GSC with covariates: beta deviates from zero (EM learned)", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_timecov,
    method = "gsc",
    r = 2,
    covariates = "cov1"
  )
  expect_true(abs(fit$beta[1]) > 1e-10)
})

# ── Phase 7b: SCM placebo test covariate consistency ─────────────────────────

test_that("fit_scm_cpp with covariates stores X0_mat (k x N_co)", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  expect_false(is.null(fit$X0_mat))
  expect_equal(dim(fit$X0_mat), c(2L, 9L))
})

test_that("mspe_ratio_pval use_covariates=TRUE returns p_value in [0,1]", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  inf <- mspe_ratio_pval(fit, use_covariates = TRUE)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  expect_equal(length(inf$mspe_ratios_all), 10L)
})

test_that("mspe_ratio_pval use_covariates=TRUE: mspe_ratios_all correct length", {
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(pred("cov1", 1:10))
  )
  inf <- mspe_ratio_pval(fit, use_covariates = TRUE)
  expect_equal(length(inf$mspe_ratios_all), 10L) # treated + 9 controls
  expect_true(is.numeric(inf$n_placebo_used))
})

test_that("mspe_ratio_pval covariate placebo matches a per-donor scm_weights_cpp loop", {
  # The covariate placebo path runs scm_placebo_x_cpp (OpenMP batch). Each
  # leave-one-out problem must be identical to an individual scm_weights_cpp
  # call on the same submatrices -- same solver core, so results agree to
  # machine precision regardless of parallel execution order.
  fit <- scm_fit(
    y ~ d | id + time,
    data = panel_cov,
    method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  inf <- mspe_ratio_pval(fit, use_covariates = TRUE)

  N_co <- ncol(fit$Y_co_pre)
  for (i in seq_len(N_co)) {
    # The default predictor-path fit runs the multi-start outer search, and
    # the placebo batch mirrors it -- so the reference loop must too.
    res_i <- scm_weights_cpp(
      fit$X0_mat[, -i, drop = FALSE], fit$X0_mat[, i],
      fit$Y_co_pre[, -i, drop = FALSE], fit$Y_co_pre[, i],
      multistart = identical(fit$v_optim_effective, "multistart")
    )
    synth_post <- fit$Y_co_post[, -i, drop = FALSE] %*% res_i$W
    expect_equal(unname(inf$placebo_effects[i]),
                 mean(fit$Y_co_post[, i] - synth_post),
                 tolerance = 1e-10)
    synth_pre <- fit$Y_co_pre[, -i, drop = FALSE] %*% res_i$W
    expect_equal(unname(inf$mspe_pre_placebo[i]),
                 mean((fit$Y_co_pre[, i] - synth_pre)^2),
                 tolerance = 1e-10)
  }
})

test_that("mspe_ratio_pval use_covariates=TRUE warns and falls back for no-cov fit", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_null(fit$X0_mat)
  expect_warning(
    inf <- mspe_ratio_pval(fit, use_covariates = TRUE),
    regexp = "no covariates"
  )
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
})

# ── Phase 9: Abadie (2021) — OOS V selection & one-sided inference ────────────

test_that("SCM v_selection='oos' returns valid fit (Abadie 2021 §3.2)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 v_selection = "oos")
  expect_s3_class(fit, "coresynth")
  expect_equal(fit$method, "scm")
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
  expect_true(all(is.finite(fit$Y_synth)))
})

test_that("SCM v_selection='oos' estimate is in reasonable range", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 v_selection = "oos")
  expect_true(abs(fit$estimate - 2.0) < 3.0)
})

test_that("SCM v_selection='insample' and 'oos' give valid V weights", {
  fit_is  <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                     v_selection = "insample")
  fit_oos <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                     v_selection = "oos")
  # In-sample: V spans all T_pre outcome rows. OOS (Abadie 2021 §3.2 / Phase
  # 26): V is selected on the training half and applied to the last
  # floor(T_pre/2) rows, so it has floor(T_pre/2) entries.
  expect_equal(length(fit_is$v_weights),  10L)
  expect_equal(length(fit_oos$v_weights), 5L)
  expect_equal(fit_oos$v_rows, 6:10)
  # Both should be valid simplex weights for V
  expect_equal(sum(fit_is$v_weights),  1, tolerance = 1e-4)
  expect_equal(sum(fit_oos$v_weights), 1, tolerance = 1e-4)
})

test_that("mspe_ratio_pval alternative='two.sided' matches default (Abadie 2021 §3.5)", {
  fit  <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf1 <- mspe_ratio_pval(fit)
  inf2 <- mspe_ratio_pval(fit, alternative = "two.sided")
  expect_equal(inf1$p_value, inf2$p_value)
  expect_equal(inf1$mspe_ratio_treated, inf2$mspe_ratio_treated)
})

test_that("mspe_ratio_pval alternative='greater' returns valid one-sided p-value", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit, alternative = "greater")
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  expect_true(is.numeric(inf$treated_effect))
  expect_equal(inf$treated_effect, fit$estimate)
})

test_that("mspe_ratio_pval alternative='less' returns valid one-sided p-value", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf_g <- mspe_ratio_pval(fit, alternative = "greater")
  inf_l <- mspe_ratio_pval(fit, alternative = "less")
  # greater and less p-values should sum to approximately 1 + 1/N (both include treated)
  expect_true(inf_l$p_value >= 0 && inf_l$p_value <= 1)
  expect_true(inf_g$p_value + inf_l$p_value <= 1 + 1/10 + 1e-6)
})

test_that("mspe_ratio_pval one-sided test is more powerful for correct direction", {
  # With a positive effect (effect=2), 'greater' p-value should be <= two-sided
  fit    <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf_2s <- mspe_ratio_pval(fit, alternative = "two.sided")
  inf_gt <- mspe_ratio_pval(fit, alternative = "greater")
  # One-sided can only be <= two-sided for the correct direction (not always guaranteed
  # due to different test statistics, but treated_effect should be accessible)
  expect_true(is.numeric(inf_gt$p_value))
  expect_true(is.numeric(inf_2s$p_value))
})

# ── Phase 10: Abadie & Zhao (2026) — Experimental SCM Design ─────────────────

# Shared helper: balanced panel without a treatment indicator (design setting)
make_design_panel <- function(J = 10, TT = 20, seed = 42) {
  set.seed(seed)
  f_t  <- cumsum(rnorm(TT, 0, 0.5))
  lam  <- rnorm(J, 1, 0.3)
  df   <- expand.grid(unit = paste0("m", seq_len(J)),
                      time = seq_len(TT),
                      stringsAsFactors = FALSE)
  df$outcome <- as.vector(outer(f_t, lam)) + rnorm(nrow(df), 0, 0.2)
  df
}

test_that("scm_design returns scm_design class with correct structure", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L)
  expect_s3_class(res, "scm_design")
  expect_equal(res$method, "scm_design")
  expect_length(res$treated_units, 1L)
  expect_length(res$tau_hat, 5L)          # 20 - 15 post periods
  expect_equal(res$T_pre, 15L)
})

test_that("scm_design w sums to 1 and v sums to 1", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L)
  expect_equal(sum(res$w), 1, tolerance = 1e-6)
  expect_equal(sum(res$v), 1, tolerance = 1e-6)
})

test_that("scm_design w and v are non-negative and disjoint", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L)
  expect_true(all(res$w >= -1e-8))
  expect_true(all(res$v >= -1e-8))
  # No unit has positive weight in both w and v
  expect_true(all(res$w * res$v < 1e-10))
  # treated and control sets are disjoint
  expect_length(intersect(res$treated_units, res$control_units), 0L)
})

test_that("scm_design inference: blank periods produce valid p_value and CI", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L, T_fit = 10L)
  expect_equal(length(res$u_blank), 5L)    # 15 - 10 blank periods
  expect_true(!is.na(res$p_value))
  expect_true(res$p_value >= 0 && res$p_value <= 1)
  expect_length(res$ci_lower, 5L)
  expect_length(res$ci_upper, 5L)
  expect_true(all(!is.na(res$ci_lower)))
  expect_true(all(res$ci_upper >= res$ci_lower))
})

test_that("scm_design default (no T_fit) gives NA inference without warning", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L)  # no warning for default
  expect_true(is.na(res$p_value))
  expect_true(all(is.na(res$ci_lower)))
})

test_that("scm_design explicit T_fit consuming all pre-periods warns", {
  df <- make_design_panel()  # TT=20, T0=15 -> T0_idx=15
  expect_warning(
    scm_design(df, "outcome", "unit", "time", T0 = 15L, T_fit = 15L),
    "No blank periods"
  )
})

test_that("scm_design m_max=2 selects 2 treated units", {
  df  <- make_design_panel(J = 10L)
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L,
                    m_min = 2L, m_max = 2L)
  expect_length(res$treated_units, 2L)
  expect_equal(sum(res$w > 1e-6), 2L)    # exactly 2 units have positive w
})

test_that("scm_design m_min=1, m_max=2 returns between 1 and 2 treated units", {
  df  <- make_design_panel(J = 8L)
  res <- scm_design(df, "outcome", "unit", "time", T0 = 12L,
                    m_min = 1L, m_max = 2L)
  n_treated <- length(res$treated_units)
  expect_true(n_treated >= 1L && n_treated <= 2L)
})

test_that("scm_design weakly_targeted variant works", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L,
                    design = "weakly_targeted", beta = 1)
  expect_equal(res$design, "weakly_targeted")
  expect_equal(sum(res$w), 1, tolerance = 1e-6)
  expect_equal(sum(res$v), 1, tolerance = 1e-6)
})

test_that("scm_design unit_level variant works and returns disjoint weights", {
  df  <- make_design_panel(J = 8L)
  res <- scm_design(df, "outcome", "unit", "time", T0 = 12L,
                    m_max = 2L, design = "unit_level", xi = 1)
  expect_equal(res$design, "unit_level")
  expect_true(all(res$w * res$v < 1e-10))
})

# Rebuild the normalised predictor matrix scm_design() uses internally
# (fitting-period outcomes, row-wise SD scaling) to evaluate the design
# objectives of eq. (9)/(10) from a fitted object.
make_design_X <- function(df, TE) {
  units <- sort(unique(df$unit))
  times <- sort(unique(df$time))
  Y <- matrix(NA_real_, length(times), length(units),
              dimnames = list(NULL, units))
  Y[cbind(match(df$time, times), match(df$unit, units))] <- df$outcome
  X <- Y[seq_len(TE), , drop = FALSE]
  row_sd <- apply(X, 1L, sd)
  row_sd[row_sd < 1e-10] <- 1
  X / row_sd
}

test_that("scm_design weakly_targeted: (w, v) are jointly optimal in eq. (9)", {
  df    <- make_design_panel(J = 8L, TT = 20L, seed = 11)
  X     <- make_design_X(df, TE = 12L)
  X_bar <- rowMeans(X)
  Vk    <- rep(1, nrow(X))

  base_term  <- function(res) {
    Xw <- drop(X[, res$treated_units, drop = FALSE] %*% res$w[res$treated_units])
    sum((X_bar - Xw)^2)
  }
  match_term <- function(res) {
    Xw <- drop(X[, res$treated_units, drop = FALSE] %*% res$w[res$treated_units])
    Xv <- drop(X[, res$control_units, drop = FALSE] %*% res$v[res$control_units])
    sum((Xw - Xv)^2)
  }

  fit_lo <- scm_design(df, "outcome", "unit", "time", T0 = 12L,
                       m_min = 2L, m_max = 2L,
                       design = "weakly_targeted", beta = 0.1)
  fit_hi <- scm_design(df, "outcome", "unit", "time", T0 = 12L,
                       m_min = 2L, m_max = 2L,
                       design = "weakly_targeted", beta = 50)

  # comparative statics: the treated/control matching term is non-increasing
  # in beta across optimal solutions
  expect_lte(match_term(fit_hi), match_term(fit_lo) + 1e-8)

  # the returned objective equals eq. (9) recomputed from the returned weights
  expect_equal(fit_hi$objective, base_term(fit_hi) + 50 * match_term(fit_hi),
               tolerance = 1e-6)

  # the joint solution is at least as good as the decoupled one (w fit to
  # X_bar alone, v best-responding) on the same treated set
  Xs <- X[, fit_hi$treated_units, drop = FALSE]
  Xc <- X[, fit_hi$control_units, drop = FALSE]
  w0 <- scm_inner_weights_cpp(Xs, X_bar, Vk)
  v0 <- scm_inner_weights_cpp(Xc, drop(Xs %*% w0), Vk)
  obj0 <- sum((X_bar - drop(Xs %*% w0))^2) +
    50 * sum((drop(Xs %*% w0) - drop(Xc %*% v0))^2)
  expect_lte(fit_hi$objective, obj0 + 1e-8)
})

test_that("scm_design unit_level: objective matches eq. (10) and responds to xi", {
  df    <- make_design_panel(J = 8L, TT = 20L, seed = 11)
  X     <- make_design_X(df, TE = 12L)
  X_bar <- rowMeans(X)
  Vk    <- rep(1, nrow(X))

  loss_term <- function(res) {
    Xc <- X[, res$control_units, drop = FALSE]
    sum(vapply(res$treated_units, function(u) {
      v_j <- scm_inner_weights_cpp(Xc, X[, u], Vk)
      res$w[[u]] * sum((X[, u] - drop(Xc %*% v_j))^2)
    }, numeric(1)))
  }

  fit_lo <- scm_design(df, "outcome", "unit", "time", T0 = 12L,
                       m_min = 2L, m_max = 2L, design = "unit_level", xi = 0.1)
  fit_hi <- scm_design(df, "outcome", "unit", "time", T0 = 12L,
                       m_min = 2L, m_max = 2L, design = "unit_level", xi = 50)

  # comparative statics: the w-weighted per-unit losses are non-increasing in xi
  expect_lte(loss_term(fit_hi), loss_term(fit_lo) + 1e-8)

  # the returned objective equals eq. (10) recomputed from the returned weights
  Xw <- drop(X[, fit_hi$treated_units, drop = FALSE] %*%
               fit_hi$w[fit_hi$treated_units])
  expect_equal(fit_hi$objective, sum((X_bar - Xw)^2) + 50 * loss_term(fit_hi),
               tolerance = 1e-6)

  # the aggregate control weights follow eq. (11): v = sum_j w_j v_j
  Xc    <- X[, fit_hi$control_units, drop = FALSE]
  v_agg <- Reduce(`+`, lapply(fit_hi$treated_units, function(u) {
    fit_hi$w[[u]] * as.numeric(scm_inner_weights_cpp(Xc, X[, u], Vk))
  }))
  expect_equal(unname(fit_hi$v[fit_hi$control_units]), v_agg, tolerance = 1e-6)
})

test_that("scm_design unit_level: large xi shifts the selected treated set", {
  # Panel where the eq. (10) trade-off flips the optimal treated set: at small
  # xi the pair best matching X_bar wins; at large xi the per-unit synthetic
  # fits dominate and a different pair is selected. The seed matters: the
  # margin between candidate subsets is small, so a panel has to be picked
  # where the flip is a genuine property of the objective rather than of how
  # accurately the per-unit QPs happen to be solved.
  set.seed(20260713)
  J <- 8L; TT <- 20L
  Fmat    <- matrix(rnorm(TT * 2L), TT, 2L)
  Lam     <- matrix(rnorm(J * 2L, 1, 0.5), J, 2L)
  delta_t <- cumsum(rnorm(TT, 0, 0.3))
  Y  <- outer(delta_t, rep(1, J)) + Fmat %*% t(Lam) +
    matrix(rnorm(TT * J, 0, 1), TT, J)
  df <- data.frame(unit = rep(paste0("u", 1:J), each = TT),
                   time = rep(1:TT, times = J),
                   y    = as.vector(Y))

  fit_lo <- scm_design(df, "y", "unit", "time", T0 = 14L, T_fit = 10L,
                       m_min = 2L, m_max = 2L, design = "unit_level", xi = 0.1)
  fit_hi <- scm_design(df, "y", "unit", "time", T0 = 14L, T_fit = 10L,
                       m_min = 2L, m_max = 2L, design = "unit_level", xi = 100)
  expect_false(setequal(fit_lo$treated_units, fit_hi$treated_units))
})

test_that("scm_design custom population weights f are normalised and applied", {
  J     <- 8L
  units <- paste0("m", seq_len(J))
  df    <- make_design_panel(J = J)
  f_raw <- setNames(c(0.4, 0.2, 0.1, 0.1, 0.1, 0.05, 0.03, 0.02), units)
  res   <- scm_design(df, "outcome", "unit", "time", T0 = 12L, f = f_raw)
  expect_equal(sum(res$f), 1, tolerance = 1e-10)
  # f stored in the object should be normalised
  expect_equal(res$f[["m1"]], 0.4 / sum(f_raw), tolerance = 1e-10)
})

test_that("scm_design rejects a pred() window with no data", {
  df <- make_design_panel(J = 8L, TT = 20L)
  expect_error(
    scm_design(df, "outcome", "unit", "time", T0 = 15L,
               predictors = list(pred("outcome", 100:110))),
    regexp = "non-finite"
  )
})

test_that("scm_design with pred() predictors works", {
  df <- make_design_panel(J = 8L, TT = 20L)
  # Add a covariate
  df$cov1 <- rnorm(nrow(df), 5, 1)
  res <- scm_design(
    df, "outcome", "unit", "time", T0 = 15L,
    predictors = list(
      pred("outcome", 1:10),
      pred("cov1",   1:15)
    )
  )
  expect_s3_class(res, "scm_design")
  expect_true(grepl("outcome", res$pred_names[1]))
  expect_true(grepl("cov1", res$pred_names[length(res$pred_names)]))
})

test_that("scm_design estimate and tau_hat are finite", {
  df  <- make_design_panel(J = 10L, TT = 25L)
  res <- scm_design(df, "outcome", "unit", "time", T0 = 18L, T_fit = 12L)
  expect_true(is.finite(res$estimate))
  expect_true(all(is.finite(res$tau_hat)))
  expect_true(all(is.finite(res$u_blank)))
})

test_that("scm_design print output mentions scm_design", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L, T_fit = 10L)
  expect_output(print(res), "scm_design")
  expect_output(print(res), "Treated units")
})

test_that("scm_design summary output shows design variant and CI", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L, T_fit = 10L)
  expect_output(summary(res), "Design variant")
  expect_output(summary(res), "Confidence intervals")
})

test_that("scm_design plot outcome type returns ggplot", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L)
  p   <- plot(res, type = "outcome")
  expect_s3_class(p, "gg")
})

test_that("scm_design plot gap type with blank periods returns ggplot with ribbon", {
  df  <- make_design_panel()
  res <- scm_design(df, "outcome", "unit", "time", T0 = 15L, T_fit = 10L)
  p   <- plot(res, type = "gap")
  expect_s3_class(p, "gg")
})

test_that("scm_design error on invalid m_max >= J", {
  df <- make_design_panel(J = 5L)
  expect_error(
    scm_design(df, "outcome", "unit", "time", T0 = 10L, m_max = 5L),
    "m_max < J"
  )
})

test_that("scm_design error on missing outcome column", {
  df <- make_design_panel()
  expect_error(
    scm_design(df, "no_col", "unit", "time", T0 = 15L),
    "not found"
  )
})

test_that("scm_design error when T0 not in time column", {
  df <- make_design_panel()
  expect_error(
    scm_design(df, "outcome", "unit", "time", T0 = 999L),
    "present in the time column"
  )
})

# ── Phase 11a: Donor pool filtering (Abadie 2021 §4) ─────────────────────────

# Helper: panel where two donors have extreme pre-treatment trajectories
make_panel_outliers <- function() {
  df <- make_panel(N = 8L, T = 20L, T_pre = 10L, seed = 42L)
  df$y[df$id == "u5"] <- df$y[df$id == "u5"] * 20  # clear outlier
  df$y[df$id == "u6"] <- df$y[df$id == "u6"] * 20  # clear outlier
  df
}

test_that("donor_mspe_threshold = Inf (default) keeps all donors", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_equal(fit$excluded_donors, character(0L))
})

test_that("donor_mspe_threshold too strict raises warning and keeps all donors", {
  # threshold = 1 means keep only donors with exactly min MSPE → < 2 donors → warning
  expect_warning(
    fit <- scm_fit(
      y ~ d | id + time, data = panel, method = "scm",
      donor_mspe_threshold = 1.0
    ),
    "ignoring filter"
  )
  expect_equal(fit$excluded_donors, character(0L))
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
})

test_that("donor_mspe_threshold excludes clear outliers", {
  df  <- make_panel_outliers()
  fit <- scm_fit(y ~ d | id + time, data = df, method = "scm",
                 donor_mspe_threshold = 5)
  expect_true(length(fit$excluded_donors) > 0L)
  expect_true("u5" %in% fit$excluded_donors || "u6" %in% fit$excluded_donors)
})

test_that("excluded donors are absent from unit_weights", {
  df  <- make_panel_outliers()
  fit <- scm_fit(y ~ d | id + time, data = df, method = "scm",
                 donor_mspe_threshold = 5)
  if (length(fit$excluded_donors) > 0L) {
    expect_false(any(fit$excluded_donors %in% names(fit$unit_weights)))
  }
})

test_that("weights remain on the simplex after donor filtering", {
  df  <- make_panel_outliers()
  fit <- scm_fit(y ~ d | id + time, data = df, method = "scm",
                 donor_mspe_threshold = 5)
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
})

test_that("ATT is finite after donor filtering", {
  df  <- make_panel_outliers()
  fit <- scm_fit(y ~ d | id + time, data = df, method = "scm",
                 donor_mspe_threshold = 5)
  expect_true(is.finite(fit$estimate))
})

test_that("donor_mspe_threshold via scm_fit() passes through correctly", {
  df  <- make_panel_outliers()
  fit <- scm_fit(y ~ d | id + time, data = df, method = "scm",
                 donor_mspe_threshold = 5)
  expect_type(fit$excluded_donors, "character")
})

# ── Phase 11b: Penalised SCM (Abadie & L'Hour 2021) ─────────────────────────

test_that("lambda_pen = NULL gives standard SCM (lambda_pen field is NA)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 lambda_pen = NULL)
  expect_true(is.na(fit$lambda_pen))
})

test_that("lambda_pen = 0 weights lie on the unit simplex", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 lambda_pen = 0)
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
})

test_that("lambda_pen = 0 records lambda_pen = 0 in fit object", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 lambda_pen = 0)
  expect_equal(fit$lambda_pen, 0)
})

test_that("lambda_pen = 1 weights lie on the unit simplex", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 lambda_pen = 1)
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
  expect_true(is.finite(fit$estimate))
})

test_that("large lambda_pen concentrates weights (fewer donors used)", {
  fit_pen  <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                      lambda_pen = 1e6)
  fit_std  <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  # Very large penalty pushes toward nearest-neighbor → max weight should increase
  expect_true(max(fit_pen$unit_weights) >= max(fit_std$unit_weights) - 1e-4)
})

test_that("lambda_pen = 'auto' returns a non-negative finite lambda", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 lambda_pen = "auto")
  expect_true(is.numeric(fit$lambda_pen))
  expect_true(is.finite(fit$lambda_pen))
  expect_true(fit$lambda_pen >= 0)
})

test_that("lambda_pen = 'auto' gives finite ATT and valid weights", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 lambda_pen = "auto")
  expect_true(is.finite(fit$estimate))
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
})

test_that("invalid lambda_pen raises an error", {
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            lambda_pen = -1),
    "non-negative"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            lambda_pen = "bad_value"),
    "NULL"
  )
})

# ── Phase 11c: Augmented SCM (Ben-Michael, Feller & Rothstein 2021) ──────────

test_that("augment_scm returns a list with required fields", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- augment_scm(fit)
  expect_type(aug, "list")
  expect_true(all(c("att_aug", "delta", "att_scm", "lambda_ridge", "beta_hat") %in% names(aug)))
})

test_that("augment_scm att_aug and delta are finite", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- augment_scm(fit)
  expect_true(is.finite(aug$att_aug))
  expect_true(is.finite(aug$delta))
})

test_that("augment_scm att_scm equals fit$estimate", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- augment_scm(fit)
  expect_equal(aug$att_scm, fit$estimate)
})

test_that("augment_scm beta_hat has length equal to T_pre", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- augment_scm(fit)
  expect_length(aug$beta_hat, fit$T_pre)
})

test_that("augment_scm auto-selected lambda_ridge is positive", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- augment_scm(fit)
  expect_true(aug$lambda_ridge > 0)
})

test_that("augment_scm with very large lambda_ridge gives delta near zero", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- augment_scm(fit, lambda_ridge = 1e12)
  # With huge lambda, beta -> 0, so delta = m_tr - W'm_co -> 0
  expect_true(abs(aug$delta) < 1e-3)
  expect_equal(aug$att_aug, aug$att_scm, tolerance = 1e-3)
})

test_that("augment_scm explicit lambda_ridge is honoured", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- augment_scm(fit, lambda_ridge = 0.5)
  expect_equal(aug$lambda_ridge, 0.5)
  expect_true(is.finite(aug$att_aug))
})

test_that("augment_scm rejects wrong method or missing Y_co_pre", {
  fit_sdid <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_error(augment_scm(fit_sdid), "method = 'scm'")

  fit_scm <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  fit_bad <- fit_scm
  fit_bad$Y_co_pre <- NULL
  expect_error(augment_scm(fit_bad), "Y_co_pre")
})

# ── Phase 12: SCM 高速化 ──────────────────────────────────────────────────────

test_that("Phase 12b: low-rank path returns valid simplex weights (k=7, N_co=46)", {
  set.seed(42L)
  k <- 7L; N_co <- 46L
  X0     <- matrix(rnorm(k * N_co), k, N_co)
  X1     <- rnorm(k)
  V_diag <- rep(1 / k, k)
  w <- scm_inner_weights_cpp(X0, X1, V_diag)
  expect_equal(length(w), N_co)
  expect_equal(sum(w), 1, tolerance = 1e-4)
  expect_true(all(w >= -1e-5))
})

test_that("Phase 12b: low-rank and full QP give algebraically equivalent weights", {
  set.seed(7L)
  k <- 5L; N_co <- 20L    # 2*k=10 < N_co=20 → low-rank path
  X0     <- matrix(rnorm(k * N_co), k, N_co)
  X1     <- rnorm(k)
  V_diag <- runif(k); V_diag <- V_diag / sum(V_diag)
  w_lr   <- scm_inner_weights_cpp(X0, X1, V_diag)
  # Reference via explicit full QP
  Q_full <- t(X0) %*% diag(V_diag) %*% X0
  c_full <- drop(t(X0) %*% diag(V_diag) %*% X1)
  w_full <- solve_simplex_qp(Q_full, c_full)
  expect_equal(w_lr, w_full, tolerance = 1e-3)
})

test_that("Phase 12b: full path used when k >= N_co/2 (k=10, N_co=16)", {
  set.seed(11L)
  k <- 10L; N_co <- 16L   # 2*k=20 >= N_co=16 → full path
  X0 <- matrix(rnorm(k * N_co), k, N_co)
  X1 <- rnorm(k)
  w  <- scm_inner_weights_cpp(X0, X1, rep(1 / k, k))
  expect_equal(sum(w), 1, tolerance = 1e-4)
  expect_true(all(w >= -1e-5))
})

test_that("Phase 12a: v_optim='bfgs' returns valid simplex weights", {
  # bfgs is deprecated (Phase 37): still functional, warns.
  expect_warning(
    fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                   v_optim = "bfgs"),
    regexp = "deprecated"
  )
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-5))
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 12a: v_optim='bfgs' + v_selection='oos' works", {
  fit <- suppressWarnings(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            v_optim = "bfgs", v_selection = "oos"))
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 12a: v_optim='coord_descent' is identical to default", {
  fit1 <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  fit2 <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                  v_optim = "coord_descent")
  expect_equal(fit1$unit_weights, fit2$unit_weights, tolerance = 1e-6)
})

test_that("Phase 12a: v_optim='auto' returns valid weights", {
  # panel is outcomes-only, so auto resolves to coord_descent.
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 v_optim = "auto")
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 12a: v_optim='bfgs' with lambda_pen works", {
  fit <- suppressWarnings(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            v_optim = "bfgs", lambda_pen = 0.5))
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(is.finite(fit$estimate))
})

# ── Phase 13a: SDID 共変量対応 ────────────────────────────────────────────────

# Time-varying covariate panel: panel_cov already has cov1/cov2 (time-invariant here;
# sufficient to test the partial-out API)
test_that("Phase 13a: SDID with one covariate returns valid fit", {
  fit <- scm_fit(y ~ d | id + time, data = panel_cov, method = "sdid",
                 covariates = "cov1")
  expect_s3_class(fit, "coresynth")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 13a: beta_hat length equals number of covariates", {
  fit <- scm_fit(y ~ d | id + time, data = panel_cov, method = "sdid",
                 covariates = c("cov1", "cov2"))
  expect_length(fit$beta_hat, 2L)
})

test_that("Phase 13a: SDID covariates — unit weights stay on simplex", {
  fit <- scm_fit(y ~ d | id + time, data = panel_cov, method = "sdid",
                 covariates = "cov1")
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(fit$unit_weights >= -1e-8))
})

test_that("Phase 13a: SDID without covariates — beta_hat is length-0", {
  fit <- scm_fit(y ~ d | id + time, data = panel_cov, method = "sdid")
  expect_length(fit$beta_hat, 0L)
})

test_that("Phase 13a: covariate partial-out changes the ATT estimate", {
  # cov1/cov2 in panel_cov are time-invariant, so they are absorbed by the
  # two-way fixed effects and (correctly) leave the SDID ATT essentially
  # unchanged. A meaningful partial-out test needs a genuinely time-varying
  # covariate that drives the outcome.
  set.seed(99)
  pc <- panel_cov
  pc$covtv <- rnorm(nrow(pc), 0, 1) + 0.5 * pc$time
  pc$y     <- pc$y + 0.8 * pc$covtv
  fit_no  <- scm_fit(y ~ d | id + time, data = pc, method = "sdid")
  fit_cov <- scm_fit(y ~ d | id + time, data = pc, method = "sdid",
                     covariates = "covtv")
  expect_false(isTRUE(all.equal(fit_no$estimate, fit_cov$estimate, tol = 1e-6)))
})

test_that("Phase 13a: staggered=FALSE for sharp SDID", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_false(fit$staggered)
})

# ── Phase 13b: SDID Staggered Adoption ───────────────────────────────────────

test_that("Phase 13b: SDID handles staggered adoption without error", {
  expect_no_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  )
})

test_that("Phase 13b: SDID staggered — staggered field is TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  expect_true(fit$staggered)
})

test_that("Phase 13b: SDID staggered — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 13b: SDID staggered — cohort_estimates has correct columns", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  expected_cols <- c("cohort", "n_treated", "T_pre", "T_post", "estimate", "weight")
  expect_true(all(expected_cols %in% names(fit$cohort_estimates)))
})

test_that("Phase 13b: SDID staggered — cohort weights sum to 1", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  expect_equal(sum(fit$cohort_estimates$weight), 1, tolerance = 1e-8)
})

test_that("Phase 13b: SDID staggered — each cohort estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  expect_true(all(is.finite(fit$cohort_estimates$estimate)))
})

test_that("Phase 13b: SDID staggered control_group='never_treated' works", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid",
                 control_group = "never_treated")
  expect_true(fit$staggered)
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 13b: SDID staggered control_group='clean' is default", {
  fit_clean  <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid",
                        control_group = "clean")
  fit_never  <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid",
                        control_group = "never_treated")
  # Different control groups should (in general) give different estimates
  expect_true(is.finite(fit_clean$estimate))
  expect_true(is.finite(fit_never$estimate))
})

test_that("Phase 13b: SDID staggered — Y_treat dimensions match panel", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  pan <- coresynth:::panel_to_matrices(staggered$y, staggered$d,
                                       staggered$id, staggered$time)
  expect_equal(nrow(fit$Y_treat), nrow(pan$Y))
})

# ── Phase 14a: SDID staggered + covariates ────────────────────────────────────

make_staggered_cov_panel <- function(N = 12, TT = 24, seed = 7) {
  p <- make_staggered_panel(N = N, TT = TT, seed = seed)
  set.seed(seed + 10L)
  p$cov1 <- NA_real_
  for (u in unique(p$id)) {
    base_val <- rnorm(1, 5, 1)
    p$cov1[p$id == u] <- base_val +
      seq(0, by = 0.05, length.out = TT) + rnorm(TT, 0, 0.1)
  }
  p
}

stag_cov <- make_staggered_cov_panel()

test_that("Phase 14a: SDID staggered + covariates — no error", {
  expect_no_error(
    scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid",
            covariates = "cov1")
  )
})

test_that("Phase 14a: SDID staggered + covariates — beta_hat length equals p", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid",
                 covariates = "cov1")
  expect_length(fit$beta_hat, 1L)
})

test_that("Phase 14a: SDID staggered + covariates — staggered=TRUE preserved", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid",
                 covariates = "cov1")
  expect_true(fit$staggered)
})

test_that("Phase 14a: SDID staggered + covariates — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid",
                 covariates = "cov1")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 14a: SDID staggered + covariates — cohort_estimates columns correct", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid",
                 covariates = "cov1")
  expected_cols <- c("cohort", "n_treated", "T_pre", "T_post", "estimate", "weight")
  expect_true(all(expected_cols %in% names(fit$cohort_estimates)))
})

test_that("Phase 14a: SDID staggered + covariates — cohort weights sum to 1", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid",
                 covariates = "cov1")
  expect_equal(sum(fit$cohort_estimates$weight), 1, tolerance = 1e-8)
})

test_that("Phase 14a: SDID staggered + covariates — covariate changes estimate", {
  fit_no_cov <- scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid")
  fit_cov    <- scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid",
                        covariates = "cov1")
  expect_false(isTRUE(all.equal(fit_no_cov$estimate, fit_cov$estimate,
                                tolerance = 1e-10)))
})

test_that("Phase 14a: SDID staggered + covariates — control_group='never_treated' works", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "sdid",
                 covariates = "cov1", control_group = "never_treated")
  expect_true(fit$staggered)
  expect_true(is.finite(fit$estimate))
})

# ── Phase 14b: SCM staggered adoption ─────────────────────────────────────────

test_that("Phase 14b: SCM handles staggered adoption without error", {
  expect_no_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  )
})

test_that("Phase 14b: SCM staggered — staggered = TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_true(fit$staggered)
})

test_that("Phase 14b: SCM staggered — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 14b: SCM staggered — cohort_estimates has required columns", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expected_cols <- c("cohort", "n_treated", "T_pre", "T_post", "estimate", "weight")
  expect_true(all(expected_cols %in% names(fit$cohort_estimates)))
})

test_that("Phase 14b: SCM staggered — cohort weights sum to 1", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_equal(sum(fit$cohort_estimates$weight), 1, tolerance = 1e-8)
})

test_that("Phase 14b: SCM staggered — each cohort unit_weights on simplex", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  for (cf in fit$cohort_fits) {
    w <- cf$unit_weights
    expect_equal(sum(w), 1, tolerance = 1e-5)
    expect_true(all(w >= -1e-5))
  }
})

test_that("Phase 14b: SCM staggered — v_selection='oos' gives finite ATT", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm",
                 v_selection = "oos")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 14b: SCM staggered — predictors != NULL errors", {
  expect_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm",
            predictors = list(pred("y", 1:3))),
    regexp = "predictors"
  )
})

test_that("Phase 14b: SCM staggered — estimate within ±4 of true ATT (1.5)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_lt(abs(fit$estimate - 1.5), 4.0)
})

test_that("Phase 14b: SCM staggered — Y_treat nrow matches panel T", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_equal(nrow(fit$Y_treat), 24L)
})

test_that("Phase 14b: SCM staggered — control_group='never_treated' gives finite ATT", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm",
                 control_group = "never_treated")
  expect_true(is.finite(fit$estimate))
})

# ── Phase 14c: augment.coresynth() ───────────────────────────────────────────

test_that("Phase 14c: augment() returns data.frame for SCM", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit)
  expect_s3_class(aug, "data.frame")
})

test_that("Phase 14c: augment() column names are correct", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit)
  expected <- c(".time", ".observed", ".fitted", ".resid", ".treated", ".period")
  expect_true(all(expected %in% names(aug)))
})

test_that("Phase 14c: augment() nrow equals T", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit)
  expect_equal(nrow(aug), 20L)
})

test_that("Phase 14c: augment() .resid == .observed - .fitted", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit)
  expect_equal(aug$.resid, aug$.observed - aug$.fitted, tolerance = 1e-10)
})

test_that("Phase 14c: augment() .treated is FALSE pre, TRUE post", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit)
  T_pre <- fit$T_pre
  expect_true(all(!aug$.treated[seq_len(T_pre)]))
  expect_true(all(aug$.treated[(T_pre + 1L):nrow(aug)]))
})

test_that("Phase 14c: augment() .period factor levels correct", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit)
  expect_equal(levels(aug$.period), c("pre", "post"))
})

test_that("Phase 14c: augment() works for SDID sharp", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  aug <- broom::augment(fit)
  expect_s3_class(aug, "data.frame")
  expect_true(all(c(".time", ".resid") %in% names(aug)))
})

test_that("Phase 14c: augment() works for GSC", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2L)
  aug <- broom::augment(fit)
  expect_s3_class(aug, "data.frame")
  expect_equal(nrow(aug), 20L)
})

test_that("Phase 14c (Phase 24-updated): augment() staggered SDID returns cohort-long df", {
  # Phase 24: augment.coresynth() now supports staggered fits and returns a
  # long-format data.frame stacking per-cohort treated paths (and donors when
  # include_donors=TRUE). The old behavior (warning + empty df) is replaced.
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  aug <- broom::augment(fit)
  expect_s3_class(aug, "data.frame")
  expect_gt(nrow(aug), 0L)
  expect_true(".cohort" %in% names(aug))
})

test_that("Phase 14c: tidy() returns cohort_estimates for staggered SCM", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  td <- broom::tidy(fit)
  expect_s3_class(td, "data.frame")
  expect_true(nrow(td) > 0L)
  expect_true("cohort" %in% names(td) || grepl("^cohort_", td$term[1]))
})

# ── Phase 15a: SDID 推論 API (sdid_inference) ─────────────────────────────────

test_that("Phase 15a: SDID sharp fit stores Y_co_pre and Y_co_post", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_equal(dim(fit$Y_co_pre),  c(10L, 9L))  # T_pre x N_co
  expect_equal(dim(fit$Y_co_post), c(10L, 9L))  # T_post x N_co
})

test_that("Phase 15a: SDID sharp fit stores Y_tr_pre_mean, Y_tr_post_mean, N_tr", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_length(fit$Y_tr_pre_mean,  10L)  # T_pre
  expect_length(fit$Y_tr_post_mean, 10L)  # T_post
  expect_equal(fit$N_tr, 1L)
})

test_that("Phase 15a: sdid_inference(placebo) returns valid structure", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  expect_equal(inf$estimate, fit$estimate)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  expect_length(inf$placebo_effects, 9L)  # N_co = 9
  # Phase 34: placebo now reports the placebo-distribution SE and a normal CI
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(is.finite(inf$ci_lower))
  expect_equal(inf$n_controls, 9L)
})

test_that("Phase 15a: sdid_inference(placebo) named placebo_effects", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  expect_equal(names(inf$placebo_effects), colnames(fit$Y_co_pre))
})

test_that("Phase 15a: sdid_inference(bootstrap) returns se and CI", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "bootstrap", n_boot = 50L, seed = 1L)
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  expect_true(inf$ci_lower < inf$ci_upper)
  expect_length(inf$boot_ests, 50L)
})

test_that("Phase 15a: sdid_inference(bootstrap) is reproducible with seed", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf1 <- sdid_inference(fit, method = "bootstrap", n_boot = 30L, seed = 99L)
  inf2 <- sdid_inference(fit, method = "bootstrap", n_boot = 30L, seed = 99L)
  expect_equal(inf1$se, inf2$se)
  expect_equal(inf1$boot_ests, inf2$boot_ests)
})

test_that("Phase 15a: sdid_inference(jackknife) returns se and CI", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "jackknife")
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  expect_true(inf$ci_lower < inf$ci_upper)
  expect_null(inf$boot_ests)
})

test_that("Phase 15a: sdid_inference alternative='greater' gives one-sided p-value", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf_g <- sdid_inference(fit, method = "placebo", alternative = "greater")
  inf_l <- sdid_inference(fit, method = "placebo", alternative = "less")
  inf_t <- sdid_inference(fit, method = "placebo", alternative = "two.sided")
  expect_equal(inf_g$alternative, "greater")
  expect_true(inf_g$p_value >= 0 && inf_g$p_value <= 1)
  # sum of one-sided p-values cannot exceed 1 + 1/N
  N_co <- inf_t$n_controls
  expect_true(inf_g$p_value + inf_l$p_value <= 1 + 1 / (N_co + 1) + 1e-9)
})

test_that("Phase 15a: sdid_inference errors on non-SDID fit", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(sdid_inference(fit), regexp = "method = 'sdid'")
})

test_that("Phase 15a: sdid_inference staggered placebo no longer errors (Phase 20)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  expect_no_error(sdid_inference(fit, method = "placebo"))
})

test_that("Phase 15a: print.sdid_inference produces output", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "jackknife")
  expect_output(print(inf), regexp = "Estimate")
  expect_output(print(inf), regexp = "p-value")
})

test_that("Phase 15a: SDID covariates fit also stores Y_co_pre/Y_co_post", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov[stag_cov$id %in%
    c("u1", paste0("u", 3:12)), ],
    method = "sdid", covariates = "cov1")
  # Sharp fit with covariates should store Y_co_pre
  if (!isTRUE(fit$staggered)) {
    expect_false(is.null(fit$Y_co_pre))
  }
})

# ── Phase 15b: SCM staggered covariates ────────────────────────────────────────

test_that("Phase 15b: SCM staggered + covariates — no error", {
  expect_no_error(
    scm_fit(y ~ d | id + time, data = stag_cov, method = "scm",
            covariates = "cov1")
  )
})

test_that("Phase 15b: SCM staggered + covariates — staggered=TRUE preserved", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "scm",
                 covariates = "cov1")
  expect_true(fit$staggered)
})

test_that("Phase 15b: SCM staggered + covariates — beta_hat length equals p", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "scm",
                 covariates = "cov1")
  expect_length(fit$beta_hat, 1L)
})

test_that("Phase 15b: SCM staggered + covariates — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "scm",
                 covariates = "cov1")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 15b: SCM staggered + covariates — cohort_estimates exist", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "scm",
                 covariates = "cov1")
  expect_s3_class(fit$cohort_estimates, "data.frame")
  expect_true(nrow(fit$cohort_estimates) > 0L)
})

test_that("Phase 15b: SCM staggered without covariates unchanged (regression)", {
  fit_no  <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  fit_cov <- scm_fit(y ~ d | id + time, data = stag_cov, method = "scm",
                     covariates = "cov1")
  # Estimates differ when covariates are included
  expect_false(isTRUE(all.equal(fit_no$estimate, fit_cov$estimate)))
})

test_that("Phase 15b: SCM staggered + predictors still errors", {
  fit_sharp <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  T_pre <- fit_sharp$T_pre
  times <- sort(unique(panel$time))
  expect_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm",
            predictors = list(pred("y", times[1:3]))),
    regexp = "predictors"
  )
})

test_that("Phase 15b: SCM sharp + covariates forwarded (no error)", {
  # Sharp SCM doesn't use covariates the same way, but passing it shouldn't error
  # (covariates is silently ignored for sharp SCM unless wired up in fit_scm_cpp)
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_s3_class(fit, "coresynth")
})

# ── Phase 15c: augment() include_donors ──────────────────────────────────────

test_that("Phase 15c: augment(include_donors=FALSE) is same as before", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug_old <- broom::augment(fit)
  aug_new <- broom::augment(fit, include_donors = FALSE)
  expect_equal(nrow(aug_new), nrow(aug_old))
  expect_equal(names(aug_new), names(aug_old))
})

test_that("Phase 15c: augment(include_donors=TRUE) adds control rows for SCM", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit, include_donors = TRUE)
  T_total <- length(fit$times)
  N_co    <- ncol(fit$Y_co_pre)
  expect_equal(nrow(aug), T_total * (1L + N_co))
})

test_that("Phase 15c: augment(include_donors=TRUE) has .unit and .type columns", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit, include_donors = TRUE)
  expect_true(".unit" %in% names(aug))
  expect_true(".type" %in% names(aug))
})

test_that("Phase 15c: augment control rows have NA .fitted and .resid", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit, include_donors = TRUE)
  ctrl <- aug[aug$.type == "control", ]
  expect_true(all(is.na(ctrl$.fitted)))
  expect_true(all(is.na(ctrl$.resid)))
})

test_that("Phase 15c: augment .type values are 'treated' and 'control'", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit, include_donors = TRUE)
  expect_true(all(aug$.type %in% c("treated", "control")))
  expect_true("treated" %in% aug$.type)
  expect_true("control" %in% aug$.type)
})

test_that("Phase 15c: augment(include_donors=TRUE) works for SDID", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  aug <- broom::augment(fit, include_donors = TRUE)
  T_total <- length(fit$times)
  N_co    <- ncol(fit$Y_co_pre)
  expect_equal(nrow(aug), T_total * (1L + N_co))
  expect_true("control" %in% aug$.type)
})

test_that("Phase 15c: augment(include_donors=TRUE) .unit names match donor names", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit, include_donors = TRUE)
  ctrl_units <- unique(aug$.unit[aug$.type == "control"])
  expect_setequal(ctrl_units, colnames(fit$Y_co_pre))
})

# ── Phase 16a: GSC staggered adoption ─────────────────────────────────────────

make_staggered_panel_multi <- function(N = 12, TT = 24, seed = 7) {
  set.seed(seed)
  f   <- cumsum(rnorm(TT, 0, 0.5))
  lam <- rnorm(N, 1, 0.3)
  rows <- expand.grid(time = seq_len(TT), id = paste0("u", seq_len(N)),
                      stringsAsFactors = FALSE)
  rows$y <- as.vector(outer(f, lam)) + rnorm(nrow(rows), 0, 0.3)
  rows$d <- as.integer(
    (rows$id %in% c("u1", "u2", "u3") & rows$time >= 9) |
    (rows$id %in% c("u4", "u5", "u6") & rows$time >= 15)
  )
  rows$y[rows$d == 1] <- rows$y[rows$d == 1] + 1.5
  rows[order(rows$id, rows$time), ]
}
staggered_multi <- make_staggered_panel_multi()

test_that("Phase 16a: GSC staggered — no error", {
  expect_no_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  )
})

test_that("Phase 16a: GSC staggered — staggered = TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_true(fit$staggered)
})

test_that("Phase 16a: GSC staggered — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 16a: GSC staggered — cohort_estimates has required columns", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expected_cols <- c("cohort", "n_treated", "T_pre", "T_post", "estimate", "weight")
  expect_true(all(expected_cols %in% names(fit$cohort_estimates)))
})

test_that("Phase 16a: GSC staggered — cohort_estimates weights sum to 1", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_equal(sum(fit$cohort_estimates$weight), 1, tolerance = 1e-8)
})

test_that("Phase 16a: GSC staggered — 2 cohorts present", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_equal(nrow(fit$cohort_estimates), 2L)
})

test_that("Phase 16a: GSC staggered — estimate within ±4 of true ATT (1.5)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_lt(abs(fit$estimate - 1.5), 4.0)
})

test_that("Phase 16a: GSC staggered — Y_treat nrow matches panel T (24)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_equal(nrow(fit$Y_treat), 24L)
})

test_that("Phase 16a: GSC staggered — per-cohort F is T x r matrix", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc", r = 2L)
  for (cf in fit$cohort_fits) {
    expect_equal(ncol(cf$F), 2L)
    expect_equal(nrow(cf$F), 24L)
  }
})

test_that("Phase 16a: GSC staggered — per-cohort L_co is N_co x r matrix", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc", r = 2L)
  for (cf in fit$cohort_fits) {
    expect_equal(ncol(cf$L_co), 2L)
    expect_true(nrow(cf$L_co) >= 2L)
  }
})

test_that("Phase 16a: GSC staggered — control_group='never_treated' gives finite ATT", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc",
                 control_group = "never_treated")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 16a: GSC staggered — multi-unit cohort gives finite ATT", {
  fit <- scm_fit(y ~ d | id + time, data = staggered_multi, method = "gsc")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 16a: GSC staggered — r=3 changes F dimensions", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc", r = 3L)
  for (cf in fit$cohort_fits) {
    expect_equal(ncol(cf$F), 3L)
  }
})

test_that("Phase 16a: GSC staggered — gsc_boot() throws error", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_error(gsc_boot(fit), regexp = "staggered")
})

test_that("Phase 16a: GSC staggered — tidy() returns cohort rows", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  td <- broom::tidy(fit)
  expect_s3_class(td, "data.frame")
  expect_true(nrow(td) == 2L)
  expect_true(grepl("^cohort_", td$term[1]))
})

# ── Phase 16b: SI staggered adoption ──────────────────────────────────────────

test_that("Phase 16b: SI staggered — no error", {
  expect_no_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "si")
  )
})

test_that("Phase 16b: SI staggered — staggered = TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expect_true(fit$staggered)
})

test_that("Phase 16b: SI staggered — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 16b: SI staggered — cohort_estimates has required columns", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expected_cols <- c("cohort", "n_treated", "T_pre", "T_post", "estimate", "weight")
  expect_true(all(expected_cols %in% names(fit$cohort_estimates)))
})

test_that("Phase 16b: SI staggered — cohort_estimates weights sum to 1", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expect_equal(sum(fit$cohort_estimates$weight), 1, tolerance = 1e-8)
})

# ── Phase 21: SI multi-arm (K>1, Agarwal et al. 2025) ────────────────────────

# DGP: tensor factor model  Y^(d)_{ti} = sum_l u_tl * Lambda[d,l] * v_il + eps
# Lambda: arm 0 = c(1,1), arm 1 = c(1.5,0.8), arm 2 = c(0.7,1.3)
make_tensor_panel <- function(N_co = 10L, N1 = 3L, N2 = 3L,
                               T = 24L, T_pre = 12L, seed = 42L) {
  set.seed(seed)
  r <- 2L
  N <- N_co + N1 + N2
  U      <- matrix(cumsum(rnorm(T * r)), T, r)
  V      <- matrix(rnorm(N * r, 0, 0.5), N, r)
  Lambda <- rbind(c(1.0, 1.0), c(1.5, 0.8), c(0.7, 1.3))
  arm_unit <- c(rep(0L, N_co), rep(1L, N1), rep(2L, N2))
  ids  <- paste0("u", seq_len(N))
  rows <- expand.grid(time = seq_len(T), id = ids, stringsAsFactors = FALSE)
  rows <- rows[order(rows$id, rows$time), ]
  i_idx <- match(rows$id, ids)
  t_idx <- rows$time
  arm_i <- arm_unit[i_idx]
  rows$d <- ifelse(arm_i == 0L, 0L, ifelse(t_idx > T_pre, arm_i, 0L))
  rows$y <- mapply(function(t, i, arm) {
    sum(U[t, ] * Lambda[arm + 1L, ] * V[i, ]) + rnorm(1L, 0, 0.05)
  }, t_idx, i_idx, arm_i)
  rows
}

tensor_panel <- make_tensor_panel()

# § A: panel_to_tensor() ──────────────────────────────────────────────────────

test_that("Phase 21: panel_to_tensor — arm_levels is c(0,1,2)", {
  tp <- coresynth:::panel_to_tensor(tensor_panel$y, tensor_panel$d,
                                     tensor_panel$id, tensor_panel$time)
  expect_equal(tp$arm_levels, c(0L, 1L, 2L))
})

test_that("Phase 21: panel_to_tensor — idx_by_arm unit counts correct", {
  tp <- coresynth:::panel_to_tensor(tensor_panel$y, tensor_panel$d,
                                     tensor_panel$id, tensor_panel$time)
  expect_equal(length(tp$idx_by_arm[["0"]]), 10L)  # N_co
  expect_equal(length(tp$idx_by_arm[["1"]]), 3L)   # N1
  expect_equal(length(tp$idx_by_arm[["2"]]), 3L)   # N2
})

test_that("Phase 21: panel_to_tensor — error when arm 0 absent", {
  bad <- tensor_panel
  bad$d <- bad$d + 1L  # shift: now arms are 1, 2, 3 (no 0)
  # The shift makes every unit "treated" from period 1, so the shared
  # no-control-units check fires before panel_to_tensor's own arm-0 check.
  expect_error(
    coresynth:::panel_to_tensor(bad$y, bad$d, bad$id, bad$time),
    "arm 0|No control units"
  )
})

test_that("Phase 21: panel_to_tensor — works on binary d (K=1)", {
  tp <- coresynth:::panel_to_tensor(panel$y, panel$d, panel$id, panel$time)
  expect_equal(tp$arm_levels, c(0L, 1L))
  expect_true(!is.null(tp$idx_by_arm[["0"]]))
  expect_true(!is.null(tp$idx_by_arm[["1"]]))
})

test_that("Phase 21: panel_to_tensor — return has required fields", {
  tp <- coresynth:::panel_to_tensor(tensor_panel$y, tensor_panel$d,
                                     tensor_panel$id, tensor_panel$time)
  for (fld in c("Y", "T_pre", "is_sharp", "arm_levels", "idx_by_arm"))
    expect_true(fld %in% names(tp), info = paste("missing field:", fld))
})

# § B: fit_si_cpp() multi-arm basic behaviour ─────────────────────────────────

test_that("Phase 21: SI multi-arm — class is coresynth", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_s3_class(fit, "coresynth")
  expect_s3_class(fit, "coresynth_si")
})

test_that("Phase 21: SI multi-arm — multi_arm flag is TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_true(isTRUE(fit$multi_arm))
})

test_that("Phase 21: SI multi-arm — arm_levels is c(1L, 2L)", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_equal(fit$arm_levels, c(1L, 2L))
})

test_that("Phase 21: SI multi-arm — arm_estimates length equals K", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_equal(length(fit$arm_estimates), 2L)
})

test_that("Phase 21: SI multi-arm — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 21: SI multi-arm — arm_fits has required sub-fields", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  for (arm in c("1", "2")) {
    af <- fit$arm_fits[[arm]]
    for (fld in c("weights", "Y_cf", "Y_treat", "Y_synth", "gap"))
      expect_true(fld %in% names(af),
                  info = paste("arm", arm, "missing field:", fld))
  }
})

test_that("Phase 21: SI multi-arm — arm 1 Y_cf dimension is T_post x N1", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  T_post <- 24L - 12L
  expect_equal(dim(fit$arm_fits[["1"]]$Y_cf), c(T_post, 3L))
})

test_that("Phase 21: SI multi-arm — arm 1 weights dimension is N_co x N1", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_equal(dim(fit$arm_fits[["1"]]$weights), c(10L, 3L))
})

# § C: estimation accuracy ────────────────────────────────────────────────────

test_that("Phase 21: SI multi-arm — arm 1 ATT is finite and bounded", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_true(is.finite(fit$arm_estimates[["1"]]))
  expect_lt(abs(fit$arm_estimates[["1"]]), 20)
})

test_that("Phase 21: SI multi-arm — arm 2 ATT is finite and bounded", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_true(is.finite(fit$arm_estimates[["2"]]))
  expect_lt(abs(fit$arm_estimates[["2"]]), 20)
})

test_that("Phase 21: SI multi-arm — manual k is respected", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si", k = 2L)
  expect_equal(fit$k, 2L)
  expect_true(is.finite(fit$estimate))
})

# § D: scm_fit() interface ───────────────────────────────────────────────────

test_that("Phase 21: SI multi-arm via scm_fit — no error with multi-valued d", {
  expect_no_error(scm_fit(y ~ d | id + time, data = tensor_panel, method = "si"))
})

test_that("Phase 21: SI multi-arm via scm_fit — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 21: SI multi-arm — print() outputs Multi-arm", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  out <- capture.output(print(fit))
  expect_true(any(grepl("Multi-arm", out)))
})

# § E: si_inference() multi-arm ───────────────────────────────────────────────

test_that("Phase 21: si_inference multi-arm bootstrap — class is coresynth_inference", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 50L, seed = 1L)
  expect_s3_class(inf, "coresynth_inference")
})

test_that("Phase 21: si_inference multi-arm bootstrap — SE finite and positive", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 50L, seed = 1L)
  expect_true(is.finite(inf$se))
  expect_gt(inf$se, 0)
})

test_that("Phase 21: si_inference multi-arm bootstrap — CI contains estimate", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 50L, seed = 1L)
  expect_lte(inf$ci_lower, inf$estimate + 1e-9)
  expect_gte(inf$ci_upper, inf$estimate - 1e-9)
})

test_that("Phase 21: si_inference multi-arm jackknife — SE finite and positive", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  inf <- si_inference(fit, method = "jackknife", seed = 1L)
  expect_true(is.finite(inf$se))
  expect_gt(inf$se, 0)
})

test_that("Phase 21: si_inference multi-arm jackknife_global — error with guidance", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_error(si_inference(fit, method = "jackknife_global"),
               "jackknife")
})

test_that("Phase 21: si_inference — rejects non-SI fit", {
  fit_mc <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  expect_error(si_inference(fit_mc), "method = 'si'")
})

# § F: edge cases ─────────────────────────────────────────────────────────────

test_that("Phase 21: SI multi-arm — aggregate ATT equals weighted arm mean", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  ws   <- vapply(fit$arm_fits, `[[`, numeric(1L), "weight")
  taus <- fit$arm_estimates
  expect_equal(sum(ws * taus) / sum(ws), fit$estimate, tolerance = 1e-10)
})

test_that("Phase 21: SI multi-arm — N_tr_d=1 arm works without error", {
  tp1 <- make_tensor_panel(N_co = 10L, N1 = 1L, N2 = 3L, seed = 99L)
  expect_no_error(scm_fit(y ~ d | id + time, data = tp1, method = "si"))
})

test_that("Phase 21: SI multi-arm — k > min(T_pre,N_co) raises error", {
  expect_error(
    scm_fit(y ~ d | id + time, data = tensor_panel, method = "si", k = 100L),
    "k"
  )
})

test_that("Phase 16b: SI staggered — 2 cohorts present", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expect_equal(nrow(fit$cohort_estimates), 2L)
})

test_that("Phase 16b: SI staggered — estimate within ±4 of true ATT (1.5)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expect_lt(abs(fit$estimate - 1.5), 4.0)
})

test_that("Phase 16b: SI staggered — Y_treat nrow matches panel T (24)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expect_equal(nrow(fit$Y_treat), 24L)
})

test_that("Phase 16b: SI staggered — control_group='never_treated' gives finite ATT", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si",
                 control_group = "never_treated")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 16b: SI staggered — multi-unit cohort gives finite ATT", {
  fit <- scm_fit(y ~ d | id + time, data = staggered_multi, method = "si")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 16b: SI staggered — k argument respected", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si", k = 2L)
  for (cf in fit$cohort_fits) {
    expect_lte(cf$k, 2L)
  }
})

test_that("Phase 16b: SI staggered — per-cohort k auto-adjusted for short pre-period", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  for (cf in fit$cohort_fits) {
    expect_gte(cf$k, 1L)
    expect_lte(cf$k, min(cf$T_pre, length(cf$idx_co)))
  }
})

test_that("Phase 16b: SI staggered — per-cohort weights matrix has correct dims", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  for (cf in fit$cohort_fits) {
    expect_equal(nrow(cf$weights), length(cf$idx_co))
    expect_equal(ncol(cf$weights), cf$n_treated)
  }
})

test_that("Phase 16b: SI staggered — tidy() returns cohort rows", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  td <- broom::tidy(fit)
  expect_s3_class(td, "data.frame")
  expect_true(nrow(td) == 2L)
  expect_true(grepl("^cohort_", td$term[1]))
})

# ── Phase 17: GSC/SI 推論 API(Bootstrap/Jackknife)────────────────────────────

# ── Phase 17a: gsc_inference — sharp ─────────────────────────────────────────

test_that("Phase 17a: gsc_inference sharp bootstrap — no error, class correct", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc")
  inf <- gsc_inference(fit, method = "bootstrap", n_boot = 99L, seed = 1L)
  expect_s3_class(inf, "coresynth_inference")
})

test_that("Phase 17a: gsc_inference sharp bootstrap — finite SE and p_value in [0,1]", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc")
  inf <- gsc_inference(fit, method = "bootstrap", n_boot = 99L, seed = 2L)
  expect_true(is.finite(inf$se))
  expect_gte(inf$p_value, 0)
  expect_lte(inf$p_value, 1)
})

test_that("Phase 17a: gsc_inference sharp bootstrap — CI contains estimate", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc")
  inf <- gsc_inference(fit, method = "bootstrap", n_boot = 99L, seed = 3L)
  expect_lte(inf$ci_lower, inf$estimate + 1e-8)
  expect_gte(inf$ci_upper, inf$estimate - 1e-8)
})

test_that("Phase 17a: gsc_inference sharp jackknife — finite SE and p_value in [0,1]", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc")
  inf <- gsc_inference(fit, method = "jackknife")
  expect_true(is.finite(inf$se))
  expect_gte(inf$p_value, 0)
  expect_lte(inf$p_value, 1)
})

test_that("Phase 17a: gsc_inference — rejects non-GSC fit", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_error(gsc_inference(fit), "method = 'gsc'")
})

# ── Phase 17b: gsc_inference — staggered ──────────────────────────────────────

test_that("Phase 17b: gsc_inference staggered bootstrap — no error, staggered=TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  inf <- gsc_inference(fit, method = "bootstrap", n_boot = 99L, seed = 10L)
  expect_s3_class(inf, "coresynth_inference")
  expect_true(inf$staggered)
})

test_that("Phase 17b: gsc_inference staggered bootstrap — finite SE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  inf <- gsc_inference(fit, method = "bootstrap", n_boot = 99L, seed = 11L)
  expect_true(is.finite(inf$se))
})

test_that("Phase 17b: gsc_inference staggered bootstrap — CI contains estimate", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  inf <- gsc_inference(fit, method = "bootstrap", n_boot = 99L, seed = 12L)
  expect_lte(inf$ci_lower, inf$estimate + 1e-8)
  expect_gte(inf$ci_upper, inf$estimate - 1e-8)
})

test_that("Phase 17b: gsc_inference staggered jackknife — finite SE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  inf <- gsc_inference(fit, method = "jackknife")
  expect_true(is.finite(inf$se))
  expect_true(inf$staggered)
})

# ── Phase 17c: si_inference — sharp ───────────────────────────────────────────

test_that("Phase 17c: si_inference sharp bootstrap — no error, class correct", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 99L, seed = 20L)
  expect_s3_class(inf, "coresynth_inference")
})

test_that("Phase 17c: si_inference sharp bootstrap — finite SE and p_value in [0,1]", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 99L, seed = 21L)
  expect_true(is.finite(inf$se))
  expect_gte(inf$p_value, 0)
  expect_lte(inf$p_value, 1)
})

test_that("Phase 17c: si_inference sharp bootstrap — CI contains estimate", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 99L, seed = 22L)
  expect_lte(inf$ci_lower, inf$estimate + 1e-8)
  expect_gte(inf$ci_upper, inf$estimate - 1e-8)
})

test_that("Phase 17c: si_inference sharp jackknife — finite SE and p_value in [0,1]", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  inf <- si_inference(fit, method = "jackknife")
  expect_true(is.finite(inf$se))
  expect_gte(inf$p_value, 0)
  expect_lte(inf$p_value, 1)
})

test_that("Phase 17c: si_inference — rejects non-SI fit", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_error(si_inference(fit), "method = 'si'")
})

# ── Phase 17d: si_inference — staggered ───────────────────────────────────────

test_that("Phase 17d: si_inference staggered bootstrap — no error, staggered=TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 99L, seed = 30L)
  expect_s3_class(inf, "coresynth_inference")
  expect_true(inf$staggered)
})

test_that("Phase 17d: si_inference staggered bootstrap — finite SE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 99L, seed = 31L)
  expect_true(is.finite(inf$se))
})

test_that("Phase 17d: si_inference staggered bootstrap — CI contains estimate", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 99L, seed = 32L)
  expect_lte(inf$ci_lower, inf$estimate + 1e-8)
  expect_gte(inf$ci_upper, inf$estimate - 1e-8)
})

test_that("Phase 17d: si_inference staggered jackknife — finite SE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  inf <- si_inference(fit, method = "jackknife")
  expect_true(is.finite(inf$se))
  expect_true(inf$staggered)
})

test_that("Phase 17d: gsc_boot staggered error mentions gsc_inference()", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  expect_error(gsc_boot(fit), "gsc_inference")
})

# ── Phase 18a: sdid_inference() staggered bootstrap / jackknife ───────────────

test_that("Phase 18a: sdid_inference staggered bootstrap returns sdid_inference", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "bootstrap", n_boot = 50L, seed = 1L)
  expect_s3_class(inf, "sdid_inference")
  expect_true(isTRUE(inf$staggered))
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
})

test_that("Phase 18a: sdid_inference staggered bootstrap CI straddles estimate", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "bootstrap", n_boot = 50L, seed = 1L)
  expect_true(inf$ci_lower <= inf$estimate && inf$estimate <= inf$ci_upper)
  expect_false(is.null(inf$boot_ests))
})

test_that("Phase 18a: sdid_inference staggered jackknife returns sdid_inference", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "jackknife")
  expect_s3_class(inf, "sdid_inference")
  expect_true(isTRUE(inf$staggered))
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
})

test_that("Phase 18a: sdid_inference staggered jackknife CI straddles estimate", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "jackknife")
  expect_true(inf$ci_lower <= inf$estimate && inf$estimate <= inf$ci_upper)
  expect_null(inf$boot_ests)
})

# ── Phase 18b: jackknife_global for SDID / GSC / SI staggered ────────────────

test_that("Phase 18b: sdid_inference jackknife_global staggered", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "jackknife_global")
  expect_s3_class(inf, "sdid_inference")
  expect_true(isTRUE(inf$staggered))
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(inf$ci_lower <= inf$estimate && inf$estimate <= inf$ci_upper)
})

test_that("Phase 18b: sdid_inference jackknife_global errors on sharp fit", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  expect_error(sdid_inference(fit, method = "jackknife_global"),
               regexp = "staggered fit")
})

test_that("Phase 18b: gsc_inference jackknife_global staggered", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "gsc")
  inf <- gsc_inference(fit, method = "jackknife_global")
  expect_s3_class(inf, "coresynth_inference")
  expect_true(isTRUE(inf$staggered))
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(inf$ci_lower <= inf$estimate && inf$estimate <= inf$ci_upper)
})

test_that("Phase 18b: si_inference jackknife_global staggered", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  inf <- si_inference(fit, method = "jackknife_global")
  expect_s3_class(inf, "coresynth_inference")
  expect_true(isTRUE(inf$staggered))
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(inf$ci_lower <= inf$estimate && inf$estimate <= inf$ci_upper)
})

# ── Phase 20: sdid_inference() staggered placebo ─────────────────────────────

test_that("Phase 20: sdid_inference staggered placebo returns sdid_inference", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  expect_s3_class(inf, "sdid_inference")
  expect_true(isTRUE(inf$staggered))
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  # Phase 34: placebo now reports the placebo-distribution SE and a normal CI
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(is.finite(inf$ci_lower))
  expect_true(is.finite(inf$ci_upper))
  expect_null(inf$boot_ests)
})

test_that("Phase 20: staggered placebo_effects length = never-treated controls", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  # never_co = intersect({u2..u12}, {u3..u12}) = {u3..u12} = 10 units
  never_co <- Reduce(intersect, lapply(fit$cohort_fits, `[[`, "idx_co"))
  expect_equal(length(inf$placebo_effects), length(never_co))
  expect_false(is.null(names(inf$placebo_effects)))
})

test_that("Phase 20: staggered placebo n_controls = never-treated count", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  never_co <- Reduce(intersect, lapply(fit$cohort_fits, `[[`, "idx_co"))
  expect_equal(inf$n_controls, length(never_co))
})

test_that("Phase 20: staggered placebo alternative one-sided p-values sum >= 1", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  p_g <- sdid_inference(fit, method = "placebo", alternative = "greater")$p_value
  p_l <- sdid_inference(fit, method = "placebo", alternative = "less")$p_value
  expect_gte(p_g + p_l, 1 - .Machine$double.eps)
  expect_lte(p_g, 1)
  expect_lte(p_l, 1)
})

test_that("Phase 20: staggered placebo with large effect yields small p-value", {
  d_strong <- staggered
  d_strong$y[d_strong$d == 1] <- d_strong$y[d_strong$d == 1] + 100
  fit <- scm_fit(y ~ d | id + time, data = d_strong, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  never_co <- Reduce(intersect, lapply(fit$cohort_fits, `[[`, "idx_co"))
  # ATT ≈ 100 >> noise; observed effect should rank first in placebo distribution
  expect_lte(inf$p_value, 2 / (length(never_co) + 1))
})

test_that("Phase 20: staggered placebo control_group=never_treated uses all controls", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid",
                 control_group = "never_treated")
  inf <- sdid_inference(fit, method = "placebo")
  # never_treated: intersection == full idx_co (same across all cohorts)
  expect_equal(inf$n_controls, length(fit$cohort_fits[[1]]$idx_co))
  expect_true(isTRUE(inf$staggered))
})

# ── Phase 23: SI multi-arm + staggered adoption ───────────────────────────────
# Two cohorts (g=5, g=7) × two arms (d=1, d=2) + N_co=4 controls, T=10.
# cohort g=5: 1 arm-1 unit + 1 arm-2 unit
# cohort g=7: 1 arm-1 unit + 1 arm-2 unit  →  4 (g,d) cells total
make_stag_tensor_panel <- function(N_co = 4L, TT = 10L, seed = 77L) {
  set.seed(seed)
  r <- 2L
  N <- N_co + 4L  # 4 treated (1 per cohort-arm cell)
  U      <- matrix(cumsum(rnorm(TT * r)), TT, r)
  V      <- matrix(rnorm(N * r, 0, 0.5), N, r)
  Lambda <- rbind(c(1.0, 1.0), c(1.5, 0.8), c(0.7, 1.3))
  # arm assignment: cols 1..N_co = arm 0, N_co+1 = arm1 g=5, N_co+2 = arm2 g=5,
  #                N_co+3 = arm1 g=7, N_co+4 = arm2 g=7
  arm_unit   <- c(rep(0L, N_co), 1L, 2L, 1L, 2L)
  cohort_of  <- c(rep(NA_integer_, N_co), 5L, 5L, 7L, 7L)
  ids  <- paste0("u", seq_len(N))
  rows <- expand.grid(time = seq_len(TT), id = ids, stringsAsFactors = FALSE)
  rows <- rows[order(rows$id, rows$time), ]
  i_idx <- match(rows$id, ids)
  rows$d <- mapply(function(i, t) {
    arm_i <- arm_unit[i]
    if (arm_i == 0L) return(0L)
    g <- cohort_of[i]
    if (!is.na(g) && t >= g) as.integer(arm_i) else 0L
  }, i_idx, rows$time)
  rows$y <- mapply(function(t, i) {
    sum(U[t, ] * Lambda[arm_unit[i] + 1L, ] * V[i, ]) + rnorm(1L, 0, 0.05)
  }, rows$time, i_idx)
  rows
}

stag_tensor_panel <- make_stag_tensor_panel()

# § A: fit structure ──────────────────────────────────────────────────────────

test_that("Phase 23: SI staggered-multi — no error", {
  expect_no_error(
    scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  )
})

test_that("Phase 23: SI staggered-multi — multi_arm=TRUE and staggered=TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  expect_true(isTRUE(fit$multi_arm))
  expect_true(isTRUE(fit$staggered))
})

test_that("Phase 23: SI staggered-multi — estimate is finite", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 23: SI staggered-multi — arm_levels is c(1L, 2L)", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  expect_equal(fit$arm_levels, c(1L, 2L))
})

test_that("Phase 23: SI staggered-multi — arm_estimates has 2 finite entries", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  expect_equal(length(fit$arm_estimates), 2L)
  expect_true(all(is.finite(fit$arm_estimates)))
})

test_that("Phase 23: SI staggered-multi — cohort_arm_estimates shape (4 rows)", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  ca  <- fit$cohort_arm_estimates
  expect_s3_class(ca, "data.frame")
  expect_equal(nrow(ca), 4L)
  for (col in c("cohort", "arm", "n_treated", "T_pre", "T_post", "estimate", "weight"))
    expect_true(col %in% names(ca))
})

test_that("Phase 23: SI staggered-multi — cohort_arm_estimates weights sum to 1", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  expect_equal(sum(fit$cohort_arm_estimates$weight), 1, tolerance = 1e-8)
})

test_that("Phase 23: SI staggered-multi — Y_all stored with correct dims", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  expect_false(is.null(fit$Y_all))
  expect_equal(nrow(fit$Y_all), 10L)                # TT
  expect_equal(ncol(fit$Y_all), 4L + 4L)            # N_co + N_treated
})

test_that("Phase 23: SI staggered-multi — control_group='never_treated' works", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si",
                 control_group = "never_treated")
  expect_true(isTRUE(fit$multi_arm))
  expect_true(isTRUE(fit$staggered))
  expect_true(is.finite(fit$estimate))
})

# § B: aggregate formula ──────────────────────────────────────────────────────

test_that("Phase 23: SI staggered-multi — ATT = sum(w*tau)/sum(w) over (g,d) cells", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  ws   <- vapply(fit$cohort_fits, `[[`, numeric(1L), "weight")
  taus <- vapply(fit$cohort_fits, `[[`, numeric(1L), "estimate")
  expect_equal(sum(ws * taus) / sum(ws), fit$estimate, tolerance = 1e-10)
})

test_that("Phase 23: SI staggered-multi — cohort_fits flat list with required fields", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  expect_type(fit$cohort_fits, "list")
  expect_gte(length(fit$cohort_fits), 1L)
  cf1 <- fit$cohort_fits[[1L]]
  for (fld in c("cohort", "arm", "n_treated", "T_pre", "T_post", "estimate",
                "weight", "idx_tr", "idx_co", "k", "weights", "Y_cf",
                "Y_treat", "gap"))
    expect_true(fld %in% names(cf1), info = paste("missing field:", fld))
})

test_that("Phase 23: SI staggered-multi — tidy() returns cohort_arm_estimate rows", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  td  <- tidy(fit)
  expect_s3_class(td, "data.frame")
  expect_equal(nrow(td), 4L)
  expect_true("type" %in% names(td))
  expect_true(all(td$type == "cohort_arm_estimate"))
})

# § C: inference — all three methods ─────────────────────────────────────────

test_that("Phase 23: si_inference staggered-multi bootstrap — class and staggered flag", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 49L, seed = 1L)
  expect_s3_class(inf, "coresynth_inference")
  expect_true(isTRUE(inf$staggered))
})

test_that("Phase 23: si_inference staggered-multi bootstrap — SE finite and CI valid", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  inf <- si_inference(fit, method = "bootstrap", n_boot = 49L, seed = 2L)
  expect_true(is.finite(inf$se))
  expect_lte(inf$ci_lower, inf$estimate)
  expect_gte(inf$ci_upper, inf$estimate)
  expect_gte(inf$p_value, 0); expect_lte(inf$p_value, 1)
})

test_that("Phase 23: si_inference staggered-multi jackknife — SE finite and CI valid", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  inf <- si_inference(fit, method = "jackknife")
  expect_true(is.finite(inf$se) && inf$se >= 0)
  expect_lte(inf$ci_lower, inf$estimate)
  expect_gte(inf$ci_upper, inf$estimate)
  expect_true(isTRUE(inf$staggered))
})

test_that("Phase 23: si_inference staggered-multi jackknife_global — SE finite and CI valid", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  inf <- si_inference(fit, method = "jackknife_global")
  expect_s3_class(inf, "coresynth_inference")
  expect_true(is.finite(inf$se) && inf$se >= 0)
  expect_lte(inf$ci_lower, inf$estimate)
  expect_gte(inf$ci_upper, inf$estimate)
  expect_true(isTRUE(inf$staggered))
})

# § D: regression — existing paths unaffected ─────────────────────────────────

test_that("Phase 23: sharp multi-arm still works (regression)", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_true(isTRUE(fit$multi_arm))
  expect_false(isTRUE(fit$staggered))
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 23: single-arm staggered still works (regression)", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  expect_true(isTRUE(fit$staggered))
  expect_false(isTRUE(fit$multi_arm))
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 23: sharp multi-arm jackknife_global still errors (regression)", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  expect_error(si_inference(fit, method = "jackknife_global"), "jackknife")
})

# ── Phase 24: applied-pattern regression tests ────────────────────────────────
# 1. SCM multi-predictor patterns (Abadie, Diamond & Hainmueller 2010 style)
# 2. tidy/glance.coresynth_inference and enriched glance.coresynth
# 3. augment.coresynth() staggered + multi-arm
# 4. Improved error messages (staggered + mspe_ratio_pval, staggered + predictors)
# 5. dead-code cleanup regression (Phase 14b/15b SCM staggered unchanged)

# § A: SCM multi-predictor patterns (Prop99 style) ───────────────────────────

test_that("Phase 24: SCM with three pred() at specific years works", {
  fit <- scm_fit(
    y ~ d | id + time, data = panel_cov, method = "scm",
    predictors = list(pred("y", 3), pred("y", 5), pred("y", 7))
  )
  expect_s3_class(fit, "coresynth")
  expect_equal(length(fit$v_weights), 3L)
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 24: SCM with mixed-aggregation pred() runs", {
  fit <- scm_fit(
    y ~ d | id + time, data = panel_cov, method = "scm",
    predictors = list(
      pred("y", 1:10),
      pred("cov1", 1:10, op = "median"),
      pred("cov2", 1:10, op = "mean")
    )
  )
  expect_s3_class(fit, "coresynth")
  expect_equal(nrow(fit$predictor_table), 3L)
})

test_that("Phase 24: SCM oos + multiple pred() returns finite ATT", {
  fit <- scm_fit(
    y ~ d | id + time, data = panel_cov, method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10),
                      pred("y", c(2, 5, 8))),
    v_selection = "oos"
  )
  expect_true(is.finite(fit$estimate))
  expect_equal(length(fit$v_weights), 3L)
})

# § B: tidy/glance.coresynth_inference ────────────────────────────────────────

test_that("Phase 24: tidy.coresynth_inference SDID sharp returns expected columns", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "jackknife")
  td  <- broom::tidy(inf)
  expect_s3_class(td, "data.frame")
  expect_equal(nrow(td), 1L)
  required <- c("term", "estimate", "std.error", "statistic", "p.value",
                "conf.low", "conf.high", "method", "alternative",
                "n_controls", "staggered")
  expect_true(all(required %in% names(td)))
  expect_equal(td$term, "ATT")
  expect_false(td$staggered)
})

test_that("Phase 24: tidy.coresynth_inference SDID staggered carries staggered=TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "jackknife")
  td  <- broom::tidy(inf)
  expect_true(td$staggered)
})

test_that("Phase 24: tidy.coresynth_inference GSC works", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "gsc", r = 2L)
  inf <- gsc_inference(fit, method = "jackknife")
  td  <- broom::tidy(inf)
  expect_equal(td$term, "ATT")
  expect_true(is.finite(td$estimate))
  expect_true(is.finite(td$std.error))
})

test_that("Phase 24: tidy.coresynth_inference SI works", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  inf <- si_inference(fit, method = "jackknife")
  td  <- broom::tidy(inf)
  expect_equal(td$term, "ATT")
  expect_true(is.finite(td$std.error))
})

test_that("Phase 24: glance.coresynth_inference returns one-row summary", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "bootstrap", n_boot = 31L, seed = 1L)
  g   <- broom::glance(inf)
  expect_equal(nrow(g), 1L)
  expect_true(all(c("method", "n_controls", "staggered", "estimate",
                    "std.error", "p.value", "conf.low", "conf.high",
                    "alternative", "n_boot_valid") %in% names(g)))
  expect_gte(g$n_boot_valid, 1L)
})

test_that("Phase 24: tidy.coresynth_inference carries placebo SE/CI (finite since Phase 34)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  td  <- broom::tidy(inf)
  expect_true(is.finite(td$std.error))
  expect_true(is.finite(td$conf.low))
  expect_true(is.finite(td$conf.high))
  expect_true(is.finite(td$p.value))
})

# § C: enriched glance.coresynth ─────────────────────────────────────────────

test_that("Phase 24: glance.coresynth returns new columns", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  g <- broom::glance(fit)
  expect_true(all(c("method", "estimate", "n_controls", "n_treated",
                    "T_pre", "T_post", "staggered", "multi_arm") %in% names(g)))
  expect_equal(g$T_pre, 10L)
  expect_equal(g$T_post, 10L)
  expect_equal(g$n_treated, 1L)
  expect_false(g$staggered)
  expect_false(g$multi_arm)
})

test_that("Phase 24: glance.coresynth staggered SDID sets staggered=TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  g <- broom::glance(fit)
  expect_true(g$staggered)
  expect_false(g$multi_arm)
})

test_that("Phase 24: glance.coresynth multi-arm SI sets multi_arm=TRUE", {
  fit <- scm_fit(y ~ d | id + time, data = tensor_panel, method = "si")
  g <- broom::glance(fit)
  expect_true(g$multi_arm)
  expect_false(g$staggered)
})

# § D: augment.coresynth staggered & multi-arm ────────────────────────────────

test_that("Phase 24: augment.coresynth() SDID staggered returns long df", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  aug <- broom::augment(fit)
  expect_s3_class(aug, "data.frame")
  expect_true(nrow(aug) > 0L)
  expect_true(all(c(".cohort", ".unit", ".type", ".time", ".observed",
                    ".fitted", ".period") %in% names(aug)))
  # SDID staggered reconstructs fitted via unit_weights + omega0
  expect_true(any(is.finite(aug$.fitted)))
})

test_that("Phase 24: augment.coresynth() SI staggered carries .cohort column", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  aug <- broom::augment(fit)
  expect_true(".cohort" %in% names(aug))
  expect_true(nrow(aug) > 0L)
  # SI cohort_fits carry Y_synth → .fitted should be finite
  expect_true(any(is.finite(aug$.fitted)))
})

test_that("Phase 24: augment.coresynth() multi-arm SI carries .arm column", {
  fit <- scm_fit(y ~ d | id + time, data = stag_tensor_panel, method = "si")
  aug <- broom::augment(fit)
  expect_true(".arm" %in% names(aug))
  expect_true(any(!is.na(aug$.arm)))
})

test_that("Phase 24: augment.coresynth() staggered include_donors adds control rows", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "si")
  aug <- broom::augment(fit, include_donors = TRUE)
  expect_true("control" %in% aug$.type)
  expect_true("treated" %in% aug$.type)
})

test_that("Phase 24: augment.coresynth() sharp unchanged (regression)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  aug <- broom::augment(fit)
  expect_false(".cohort" %in% names(aug))
  expect_true(all(c(".time", ".observed", ".fitted", ".resid") %in% names(aug)))
})

# § E: improved error messages ────────────────────────────────────────────────

test_that("Phase 24: mspe_ratio_pval() rejects staggered fits with helpful message", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_error(mspe_ratio_pval(fit), regexp = "staggered")
})

test_that("Phase 24: SCM staggered + predictors error message names alternatives", {
  err <- tryCatch(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm",
            predictors = list(pred("y", 1:3))),
    error = function(e) conditionMessage(e)
  )
  expect_true(grepl("covariates", err))
  expect_true(grepl("sdid|gsc|mc|tasc", err))
})

# § F: dead-code cleanup regression ───────────────────────────────────────────

test_that("Phase 24: SCM staggered unchanged after max_iter/tol cleanup", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_true(isTRUE(fit$staggered))
  expect_true(is.finite(fit$estimate))
  expect_equal(sum(fit$cohort_estimates$weight), 1, tolerance = 1e-8)
})

test_that("Phase 24: SCM staggered + covariates unchanged after cleanup", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "scm",
                 covariates = "cov1")
  expect_true(isTRUE(fit$staggered))
  expect_true(is.finite(fit$estimate))
  expect_length(fit$beta_hat, 1L)
})

test_that("Phase 24: sharp SCM tidy() unchanged (regression)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  td <- broom::tidy(fit)
  expect_true("unit_weight" %in% td$type)
})

test_that("Phase 24: sdid_inference object inherits coresynth_inference", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  expect_s3_class(inf, "sdid_inference")
  expect_s3_class(inf, "coresynth_inference")
})

# ── Conformal inference (Chernozhukov, Wuthrich & Zhu 2021) ───────────────────
test_that("conformal_inference works across supported sharp methods", {
  for (m in c("scm", "sdid", "gsc", "mc", "si")) {
    fit <- scm_fit(y ~ d | id + time, data = panel, method = m)
    ci  <- conformal_inference(fit, n_grid = 60L)
    expect_s3_class(ci, "conformal_inference")
    expect_s3_class(ci, "coresynth_inference")
    expect_equal(ci$method, "conformal")
    expect_true(ci$p_value >= 0 && ci$p_value <= 1)
    expect_true(is.finite(ci$estimate))
    expect_true(ci$n_controls == 9L)
  }
})

test_that("conformal p-value for true effect exceeds p-value for tau0 = 0", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  p_null   <- conformal_inference(fit, tau0 = 0, ci = FALSE)$p_value
  p_truth  <- conformal_inference(fit, tau0 = fit$estimate, ci = FALSE)$p_value
  # The estimate itself should not be rejected; tau0 = 0 (no effect) should be
  # less compatible than the point estimate.
  expect_true(p_truth >= p_null)
})

test_that("conformal CI contains the point estimate", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  ci  <- conformal_inference(fit, n_grid = 120L)
  if (!is.na(ci$ci_lower)) {
    expect_true(ci$ci_lower <= fit$estimate && fit$estimate <= ci$ci_upper)
  }
  expect_true(is.numeric(ci$grid) && length(ci$grid) == 120L)
})

test_that("conformal one-sided alternatives run and tidy()/glance() work", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "si")
  ci_g <- conformal_inference(fit, alternative = "greater", ci = FALSE)
  ci_l <- conformal_inference(fit, alternative = "less", ci = FALSE)
  expect_true(ci_g$p_value >= 0 && ci_g$p_value <= 1)
  expect_true(ci_l$p_value >= 0 && ci_l$p_value <= 1)
  td <- broom::tidy(conformal_inference(fit, n_grid = 40L))
  expect_equal(td$term, "ATT")
  expect_equal(td$method, "conformal")
  gl <- broom::glance(conformal_inference(fit, n_grid = 40L))
  expect_equal(gl$method, "conformal")
})

test_that("conformal_inference rejects staggered and tasc fits", {
  stag <- make_panel()
  stag$d <- 0L
  stag$d[stag$id == "u1" & stag$time > 10] <- 1L
  stag$d[stag$id == "u2" & stag$time > 15] <- 1L
  stag$y[stag$d == 1] <- stag$y[stag$d == 1] + 2.0
  fit_stag <- scm_fit(y ~ d | id + time, data = stag, method = "sdid")
  expect_error(conformal_inference(fit_stag), "sharp")

  fit_tasc <- scm_fit(y ~ d | id + time, data = panel, method = "tasc")
  expect_error(conformal_inference(fit_tasc), "tasc")
})

# ── Phase 26: Abadie 理論照合 — predictor scaling・OOS CV・ADH 2015 検証ツール ─

test_that("Phase 26: predictor scaling makes W invariant to predictor units", {
  pan_sc <- panel_cov
  pan_sc$cov1 <- pan_sc$cov1 * 1000  # same information, different units
  fit_base <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                      predictors = list(pred(c("cov1", "cov2"), 1:10)))
  fit_resc <- scm_fit(y ~ d | id + time, data = pan_sc, method = "scm",
                      predictors = list(pred(c("cov1", "cov2"), 1:10)))
  expect_equal(fit_base$unit_weights, fit_resc$unit_weights, tolerance = 1e-6)
  expect_equal(fit_base$estimate, fit_resc$estimate, tolerance = 1e-6)
})

test_that("Phase 26: scale_predictors = FALSE disables Synth-style scaling", {
  fit_off <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                     predictors = list(pred(c("cov1", "cov2"), 1:10)),
                     scale_predictors = FALSE)
  expect_s3_class(fit_off, "coresynth")
  expect_equal(sum(fit_off$unit_weights), 1, tolerance = 1e-4)
})

test_that("Phase 26: predictor_table reports original (unscaled) values", {
  fit <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                 predictors = list(pred(c("cov1", "cov2"), 1:10)))
  tr_cov1 <- mean(panel_cov$cov1[panel_cov$id == "u1" & panel_cov$time %in% 1:10])
  expect_equal(fit$predictor_table$treated[1], tr_cov1, tolerance = 1e-8)
})

test_that("Phase 26: build_predictor_matrices rejects non-finite predictors", {
  expect_error(
    scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
            predictors = list(pred("cov1", 100:110))),
    regexp = "non-finite"
  )
})

test_that("Phase 26: OOS V selection uses train/validation split (no leakage)", {
  fit_oos <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                     v_selection = "oos")
  # V applies to the last floor(T_pre/2) outcome rows (Abadie 2021 §3.2 step 4)
  expect_equal(names(fit_oos$v_weights), paste0("V", 6:10))
  expect_equal(fit_oos$v_rows, 6:10)
  expect_equal(sum(fit_oos$unit_weights), 1, tolerance = 1e-4)
  expect_true(all(is.finite(fit_oos$Y_synth)))
})

test_that("Phase 26: OOS + bfgs uses the same train/validation split", {
  fit <- suppressWarnings(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            v_selection = "oos", v_optim = "bfgs"))
  expect_equal(length(fit$v_weights), 5L)
  expect_true(all(is.finite(fit$Y_synth)))
})

test_that("Phase 26: OOS + lambda_pen builds the QP on the last-half window", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 v_selection = "oos", lambda_pen = 0.1)
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-4)
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 26: placebo_in_time shows no effect when backdated (ADH 2015)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  pit <- placebo_in_time(fit, t0_placebo = 5L)
  expect_equal(pit$t0_placebo, 5L)
  expect_equal(length(pit$gap), 10L)
  expect_true(is.finite(pit$placebo_att))
  # no treatment occurs in the pre-period, so the placebo ATT should be
  # small relative to the true post-period effect (2.0)
  expect_lt(abs(pit$placebo_att), 1.5)
})

test_that("Phase 26: placebo_in_time validates inputs", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(placebo_in_time(fit, t0_placebo = 1L), "t0_placebo")
  expect_error(placebo_in_time(fit, t0_placebo = 10L), "t0_placebo")
  stag_fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_error(placebo_in_time(stag_fit), "sharp")
})

test_that("Phase 26: loo_donors returns one row per contributing donor", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  loo <- loo_donors(fit)
  expect_gte(nrow(loo$results), 1L)
  expect_true(all(is.finite(loo$results$att_loo)))
  expect_true(all(loo$results$weight > 1e-6))
  expect_equal(loo$att_original, fit$estimate)
  expect_true(all(c("donor", "weight", "att_loo") %in% names(loo$results)))
})

test_that("Phase 26: loo_donors works for OOS and predictor fits", {
  fit_oos <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                     v_selection = "oos")
  loo_oos <- loo_donors(fit_oos)
  expect_true(all(is.finite(loo_oos$results$att_loo)))

  fit_cov <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                     predictors = list(pred(c("cov1", "cov2"), 1:10)))
  loo_cov <- loo_donors(fit_cov)
  expect_true(all(is.finite(loo_cov$results$att_loo)))

  stag_fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_error(loo_donors(stag_fit), "sharp")
})

# == Input validation: empty / degenerate / unbalanced panels =================

test_that("0-row data errors with a clear message for every method", {
  empty <- panel[0, ]
  for (m in c("scm", "sdid", "gsc", "mc", "tasc", "si")) {
    expect_error(
      scm_fit(y ~ d | id + time, data = empty, method = m),
      "0 rows", info = m
    )
  }
})

test_that("data with no treated units errors clearly", {
  no_tr <- panel[panel$id != "u1", ]
  expect_error(
    scm_fit(y ~ d | id + time, data = no_tr, method = "scm"),
    "No treated units"
  )
})

test_that("sharp data with no control units errors clearly", {
  all_tr <- panel
  all_tr$d <- as.integer(all_tr$time > 10)
  expect_error(
    scm_fit(y ~ d | id + time, data = all_tr, method = "scm"),
    "No control units"
  )
})

test_that("staggered panel with no never-treated units still fits via clean controls", {
  stag_all <- make_panel()
  stag_all$d <- as.integer(
    (stag_all$id %in% paste0("u", 1:3) & stag_all$time > 10) |
    (!stag_all$id %in% paste0("u", 1:3) & stag_all$time > 16)
  )
  expect_warning(
    fit <- scm_fit(y ~ d | id + time, data = stag_all, method = "scm"),
    "skipped"
  )
  expect_true(isTRUE(fit$staggered))
  expect_true(is.finite(fit$estimate))
})

test_that("NA in the treatment indicator errors instead of silently dropping the unit", {
  na_d <- panel
  na_d$d[57] <- NA
  expect_error(
    scm_fit(y ~ d | id + time, data = na_d, method = "scm"),
    "NA value"
  )
})

test_that("negative treatment values are rejected", {
  neg <- panel
  neg$d[1] <- -1L
  expect_error(
    scm_fit(y ~ d | id + time, data = neg, method = "scm"),
    "non-negative"
  )
})

test_that("NA unit identifiers are rejected", {
  na_id <- panel
  na_id$id[3] <- NA
  expect_error(
    scm_fit(y ~ d | id + time, data = na_id, method = "scm"),
    "identifiers contain NA"
  )
})

test_that("missing (id, time) cells error clearly for dense-matrix methods", {
  unbal <- panel[!(panel$id == "u5" & panel$time == 10), ]
  for (m in c("scm", "sdid", "gsc", "si")) {
    expect_error(
      scm_fit(y ~ d | id + time, data = unbal, method = m),
      "balanced panel with fully observed outcomes", info = m
    )
  }
})

test_that("NA outcome values error the same way as missing cells", {
  na_y <- panel
  na_y$y[na_y$id == "u5" & na_y$time == 10] <- NA
  expect_error(
    scm_fit(y ~ d | id + time, data = na_y, method = "scm"),
    "balanced panel with fully observed outcomes"
  )
})

test_that("MC masks missing cells instead of treating them as observed zeros", {
  unbal <- panel[!(panel$id == "u5" & panel$time == 10), ]
  fit_full  <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  fit_unbal <- scm_fit(y ~ d | id + time, data = unbal, method = "mc")
  expect_true(is.finite(fit_unbal$estimate))
  expect_lt(abs(fit_unbal$estimate - fit_full$estimate), 0.5)
})

test_that("TASC fits a panel with missing outcome cells", {
  unbal <- panel[!(panel$id == "u5" & panel$time == 10), ]
  fit <- scm_fit(y ~ d | id + time, data = unbal, method = "tasc")
  expect_true(is.finite(fit$estimate))
})

test_that("factor treatment and non-numeric outcome columns are rejected", {
  fac_d <- panel
  fac_d$d <- factor(fac_d$d)
  expect_error(
    scm_fit(y ~ d | id + time, data = fac_d, method = "scm"),
    "must be numeric, integer, or logical"
  )
  chr_y <- panel
  chr_y$y <- as.character(chr_y$y)
  expect_error(
    scm_fit(y ~ d | id + time, data = chr_y, method = "scm"),
    "must be numeric"
  )
})

test_that("non-integer treatment values are rejected", {
  frac_d <- panel
  frac_d$d <- frac_d$d * 0.5
  expect_error(
    scm_fit(y ~ d | id + time, data = frac_d, method = "scm"),
    "non-integer"
  )
})

test_that("SCM with zero pre-treatment periods errors clearly", {
  t0 <- panel
  t0$d <- as.integer(t0$id == "u1")
  expect_error(
    scm_fit(y ~ d | id + time, data = t0, method = "scm"),
    "No pre-treatment periods"
  )
})

test_that("scm_design rejects 0-row data", {
  expect_error(
    scm_design(panel[0, ], outcome = "y", unit = "id", time = "time", T0 = 10),
    "0 rows"
  )
})

# ── Phase 28: in-space placebo gap paths + plot (ADH 2010, Figures 4-8) ───────

test_that("scm_placebo_cpp returns gap paths consistent with effects", {
  fit  <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  plac <- scm_placebo_cpp(fit$Y_co_pre, fit$Y_co_post)
  expect_equal(dim(plac$gaps), c(20L, 9L))
  expect_true(all(is.finite(plac$gaps)))
  # Mean of the post-period gap rows must reproduce the placebo effects
  expect_equal(colMeans(plac$gaps[11:20, , drop = FALSE]), as.numeric(plac$effects),
               tolerance = 1e-10)
  # Mean squared pre-period gap rows must reproduce mspe_pre
  expect_equal(colMeans(plac$gaps[1:10, , drop = FALSE]^2), as.numeric(plac$mspe_pre),
               tolerance = 1e-10)
})

test_that("mspe_ratio_pval returns scm_placebo object with gap paths", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  expect_s3_class(inf, "scm_placebo")
  expect_equal(dim(inf$gaps), c(20L, 9L))
  expect_equal(colnames(inf$gaps), sort(paste0("u", 2:10)))
  expect_equal(length(inf$treated_gap), 20L)
  expect_equal(inf$treated_gap, as.numeric(fit$gap))
  expect_equal(length(inf$mspe_pre_placebo), 9L)
  expect_true(is.finite(inf$mspe_pre_treated))
  expect_equal(inf$T_pre, 10L)
  expect_equal(inf$times, fit$times)
  # placebo_effects must equal post-period column means of the gap paths
  expect_equal(as.numeric(inf$placebo_effects), unname(colMeans(inf$gaps[11:20, ])),
               tolerance = 1e-10)
})

test_that("mspe_ratio_pval use_covariates=TRUE also returns gap paths", {
  fit <- scm_fit(
    y ~ d | id + time, data = panel, method = "scm",
    predictors = list(pred("y", 1:5, "mean"), pred("y", 6:10, "mean"))
  )
  inf <- mspe_ratio_pval(fit, use_covariates = TRUE)
  expect_s3_class(inf, "scm_placebo")
  expect_equal(dim(inf$gaps), c(20L, 9L))
  expect_true(all(is.finite(inf$gaps)))
  expect_equal(as.numeric(inf$placebo_effects), unname(colMeans(inf$gaps[11:20, ])),
               tolerance = 1e-10)
})

test_that("plot.scm_placebo type='gaps' returns a ggplot and prunes by MSPE multiple", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)

  p_all <- plot(inf, type = "gaps")
  expect_s3_class(p_all, "ggplot")
  expect_equal(length(unique(p_all$layers[[1]]$data$unit)), 9L)

  # ADH 2010 Figures 5-7 style pruning: only well-fitted placebos remain
  mult    <- 2
  n_keep  <- sum(inf$mspe_pre_placebo <= mult * inf$mspe_pre_treated)
  p_prune <- plot(inf, type = "gaps", mspe_prune = mult)
  expect_s3_class(p_prune, "ggplot")
  expect_equal(length(unique(p_prune$layers[[1]]$data$unit)), n_keep)
})

test_that("plot.scm_placebo warns when all placebos are pruned", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  expect_warning(p <- plot(inf, type = "gaps", mspe_prune = 1e-12),
                 "pruned")
  expect_s3_class(p, "ggplot")
})

test_that("plot.scm_placebo type='ratios' returns a ggplot with all units", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  p <- plot(inf, type = "ratios")
  expect_s3_class(p, "ggplot")
  expect_equal(nrow(p$data), 10L) # treated + 9 donors
  expect_true("Treated" %in% p$data$unit)
})

test_that("plot.scm_placebo rejects invalid mspe_prune", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  expect_error(plot(inf, mspe_prune = -1), "positive")
  expect_error(plot(inf, mspe_prune = c(1, 2)), "single")
})

# ── Phase 29: plot style customization (colors/vline/hline/fill) ────────────

test_that("plot.coresynth trend: color override reaches the built plot data", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  p_default <- plot(fit, type = "trend")
  expect_equal(length(p_default$layers), 2L) # geom_line + geom_vline

  p_custom <- plot(fit, type = "trend", colors = c(treated = "black"))
  built <- ggplot2::ggplot_build(p_custom)
  expect_true("black" %in% built$data[[1]]$colour)
  expect_true("#d73027" %in% built$data[[1]]$colour) # Synthetic Control unchanged
})

test_that("plot.coresynth trend: unknown color name errors", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(plot(fit, type = "trend", colors = c(Bogus = "red")),
               "unrecognized name")
})

test_that("plot.coresynth trend: vline = FALSE/NULL suppresses the line", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  n_default <- length(plot(fit, type = "trend")$layers)
  expect_equal(length(plot(fit, type = "trend", vline = FALSE)$layers), n_default - 1L)
  expect_equal(length(plot(fit, type = "trend", vline = NULL)$layers), n_default - 1L)
})

test_that("plot.coresynth trend: vline must be a list", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(plot(fit, type = "trend", vline = "red"),
               "list of aesthetic overrides")
})

test_that("plot.coresynth gap: color/vline/hline overrides work", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  n_default <- length(plot(fit, type = "gap")$layers)
  expect_equal(length(plot(fit, type = "gap", vline = NULL, hline = NULL)$layers),
               n_default - 2L)

  p_color <- plot(fit, type = "gap", colors = "black")
  built   <- ggplot2::ggplot_build(p_color)
  expect_true(all(built$data[[1]]$colour == "black"))
})

test_that("plot.coresynth weights: fill override reaches the built plot data", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  p <- plot(fit, type = "weights", fill = "darkorange")
  built <- ggplot2::ggplot_build(p)
  expect_true(all(built$data[[1]]$fill == "darkorange"))
})

test_that("plot.coresynth weights: top_n limits bars to the largest weights", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  n_all <- nrow(ggplot2::ggplot_build(plot(fit, type = "weights"))$data[[1]])

  built <- ggplot2::ggplot_build(plot(fit, type = "weights", top_n = 2))
  expect_equal(nrow(built$data[[1]]), 2L)
  w_kept <- sort(fit$unit_weights[fit$unit_weights > 1e-4], decreasing = TRUE)
  expect_setequal(built$data[[1]]$y, unname(w_kept[1:2]))

  # top_n beyond the donor count keeps everything (default Inf unchanged)
  expect_equal(
    nrow(ggplot2::ggplot_build(plot(fit, type = "weights", top_n = 999))$data[[1]]),
    n_all
  )

  expect_error(plot(fit, type = "weights", top_n = 0), "top_n")
  expect_error(plot(fit, type = "weights", top_n = c(1, 2)), "top_n")
})

test_that("plot.coresynth pred_weights: one bar per predictor, fill reaches plot", {
  fit <- scm_fit(
    y ~ d | id + time, data = panel_cov, method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  p <- plot(fit, type = "pred_weights", fill = "darkorange")
  expect_s3_class(p, "ggplot")
  built <- ggplot2::ggplot_build(p)
  expect_equal(nrow(built$data[[1]]), length(fit$v_weights))
  expect_true(all(built$data[[1]]$fill == "darkorange"))
  # bar heights are the V weights (order-independent)
  expect_setequal(built$data[[1]]$y, unname(fit$v_weights))
})

test_that("plot.coresynth pred_weights: top_n keeps the largest V weights", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  n_all <- nrow(ggplot2::ggplot_build(plot(fit, type = "pred_weights"))$data[[1]])
  expect_equal(n_all, length(fit$v_weights))

  built <- ggplot2::ggplot_build(plot(fit, type = "pred_weights", top_n = 2))
  expect_equal(nrow(built$data[[1]]), 2L)
  v_kept <- sort(fit$v_weights, decreasing = TRUE)[1:2]
  expect_setequal(built$data[[1]]$y, unname(v_kept))

  expect_error(plot(fit, type = "pred_weights", top_n = 0), "top_n")
})

test_that("plot.coresynth pred_weights: errors when no V matrix exists", {
  fit_mc <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  expect_error(plot(fit_mc, type = "pred_weights"), "V")

  # staggered SCM stores v_weights = NULL
  fit_stag <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_error(plot(fit_stag, type = "pred_weights"), "V")
})

test_that("plot.scm_placebo type='gaps': color/vline/hline overrides work", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)

  n_default <- length(plot(inf, type = "gaps")$layers)
  expect_equal(length(plot(inf, type = "gaps", vline = FALSE, hline = FALSE)$layers),
               n_default - 2L)

  p_custom <- plot(inf, type = "gaps",
                    colors = c(placebo = "lightblue"))
  built <- ggplot2::ggplot_build(p_custom)
  expect_true("lightblue" %in% built$data[[1]]$colour)
  expect_true("#2166ac" %in% built$data[[2]]$colour) # Treated unchanged
})

test_that("plot.scm_placebo: unknown color name errors", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  expect_error(plot(inf, type = "gaps", colors = c(Bogus = "red")),
               "unrecognized name")
  expect_error(plot(inf, type = "ratios", colors = c(Bogus = "red")),
               "unrecognized name")
})

test_that("plot.coresynth trend: labels relabel the legend on both scales", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  p <- plot(fit, type = "trend", labels = c(treated = "California"))
  expect_equal(p$scales$get_scales("colour")$labels[["Treated"]], "California")
  expect_equal(p$scales$get_scales("linetype")$labels[["Treated"]], "California")
  # unmentioned series keeps its default label
  expect_equal(p$scales$get_scales("colour")$labels[["Synthetic Control"]],
               "Synthetic Control")
})

test_that("plot.coresynth trend: labels validation", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(plot(fit, type = "trend", labels = c(Bogus = "x")),
               "unrecognized name")
  expect_error(plot(fit, type = "trend", labels = "x"), "named vector")
})

test_that("plot.scm_placebo: labels relabel the legend and the ratios axis tick", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)

  p_gaps <- plot(inf, type = "gaps", labels = c(placebo = "Donors"))
  expect_equal(p_gaps$scales$get_scales("colour")$labels[["Placebo (donor pool)"]],
               "Donors")

  p_ratios <- plot(inf, type = "ratios", labels = c(treated = "California"))
  expect_equal(p_ratios$scales$get_scales("colour")$labels[["Treated"]], "California")
  built <- ggplot2::ggplot_build(p_ratios)
  expect_true("California" %in% built$layout$panel_params[[1]]$y$get_labels())

  expect_error(plot(inf, type = "gaps", labels = c(Bogus = "x")),
               "unrecognized name")
})

test_that("plot.coresynth trend: linetypes override the series line types", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")

  built_default <- ggplot2::ggplot_build(plot(fit, type = "trend"))$data[[1]]
  expect_setequal(built_default$linetype, c("solid", "dashed"))

  # both series solid: the point of the argument
  built <- ggplot2::ggplot_build(
    plot(fit, type = "trend", linetypes = c(synthetic = "solid")))$data[[1]]
  expect_true(all(built$linetype == "solid"))

  # unmentioned series keeps its default, and integer codes are accepted
  built_num <- ggplot2::ggplot_build(
    plot(fit, type = "trend", linetypes = c(treated = 3)))$data[[1]]
  expect_setequal(built_num$linetype, c("dotted", "dashed"))
})

test_that("plot.coresynth trend: linetypes compose with colors and labels", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  p <- plot(fit, type = "trend",
            colors    = c(treated = "black"),
            labels    = c(treated = "California"),
            linetypes = c(treated = "dashed", synthetic = "solid"))
  built <- ggplot2::ggplot_build(p)$data[[1]]
  expect_true("black" %in% built$colour)
  expect_setequal(built$linetype, c("dashed", "solid"))
  # legend stays merged: identical labels on the color and linetype scales
  expect_equal(p$scales$get_scales("linetype")$labels[["Treated"]], "California")
  expect_equal(p$scales$get_scales("colour")$labels[["Treated"]], "California")
})

test_that("plot.coresynth trend: linetypes validation", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(plot(fit, type = "trend", linetypes = c(Bogus = "solid")),
               "unrecognized name")
  expect_error(plot(fit, type = "trend", linetypes = "solid"), "named vector")
  expect_error(plot(fit, type = "trend", linetypes = c(treated = 9)),
               "integer code")
  # the donors key exists only while donor paths are drawn
  expect_error(plot(fit, type = "trend", linetypes = c(donors = "solid")),
               "unrecognized name")
  p <- plot(fit, type = "trend", show_donors = 2,
            linetypes = c(donors = "dotted"))
  expect_equal(ggplot2::ggplot_build(p)$data[[1]]$linetype[1], "dotted")
})

test_that("plot.coresynth gap: linetypes take a single line type", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_true(all(
    ggplot2::ggplot_build(plot(fit, type = "gap"))$data[[1]]$linetype == "solid"))

  built <- ggplot2::ggplot_build(
    plot(fit, type = "gap", linetypes = "dashed"))$data[[1]]
  expect_true(all(built$linetype == "dashed"))

  built_num <- ggplot2::ggplot_build(
    plot(fit, type = "gap", linetypes = 3))$data[[1]]
  expect_true(all(built_num$linetype == "dotted"))

  expect_error(plot(fit, type = "gap", linetypes = 9), "integer code")
})

test_that("plot.scm_placebo gaps: linetypes override the series line types", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)

  built_default <- ggplot2::ggplot_build(plot(inf, type = "gaps"))$data
  expect_true(all(built_default[[1]]$linetype == "solid"))
  expect_true(all(built_default[[2]]$linetype == "solid"))

  p <- plot(inf, type = "gaps", linetypes = c(placebo = "dotted"),
            labels = c(placebo = "Donors"))
  built <- ggplot2::ggplot_build(p)$data
  expect_true(all(built[[1]]$linetype == "dotted"))  # placebo paths
  expect_true(all(built[[2]]$linetype == "solid"))   # treated path unchanged
  # merged legend: both scales relabel the same series
  expect_equal(p$scales$get_scales("linetype")$labels[["Placebo (donor pool)"]],
               "Donors")

  expect_error(plot(inf, type = "gaps", linetypes = c(Bogus = "solid")),
               "unrecognized name")
})

# ── plot_data(): tidy extraction of the plot() data ──────────────────────────

test_that("plot_data(type='trend') mirrors the treated/synthetic accessors", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  df  <- plot_data(fit, type = "trend")
  expect_s3_class(df, "data.frame")
  expect_identical(names(df), c("time", "value", "series"))
  expect_setequal(unique(df$series), c("Treated", "Synthetic Control"))
  expect_equal(df$value[df$series == "Treated"],
               as.numeric(treated_outcomes(fit, na.rm = TRUE)))
  expect_equal(df$value[df$series == "Synthetic Control"],
               as.numeric(synthetic_outcomes(fit, na.rm = TRUE)))
})

test_that("plot_data(type='gap') equals treated minus synthetic", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  df  <- plot_data(fit, type = "gap")
  expect_identical(names(df), c("time", "gap"))
  expect_equal(df$gap,
               as.numeric(treated_outcomes(fit, na.rm = TRUE) -
                            synthetic_outcomes(fit, na.rm = TRUE)))
})

test_that("plot_data trend align shifts synthetic by the pre-period gap", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  raw <- plot_data(fit, type = "gap")
  ali <- plot_data(fit, type = "gap", align = TRUE)
  # SDID: lambda-aligned post-period mean gap equals the point estimate
  post <- seq(fit$T_pre + 1L, length(ali$gap))
  expect_equal(mean(ali$gap[post]), unname(fit$estimate), tolerance = 1e-6)
  expect_false(isTRUE(all.equal(raw$gap, ali$gap)))
})

test_that("plot_data trend show_donors adds Donors rows and a unit column", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  df  <- plot_data(fit, type = "trend", show_donors = 3)
  expect_true("unit" %in% names(df))
  expect_setequal(unique(df$series),
                  c("Treated", "Synthetic Control", "Donors"))
  # exactly three distinct donor units
  expect_equal(length(unique(df$unit[df$series == "Donors"])), 3L)
})

test_that("plot_data(type='weights') returns every donor (no 1e-4 pruning)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  df  <- plot_data(fit, type = "weights")
  expect_identical(names(df), c("unit", "weight"))
  # complete: one row per donor, unlike the plot which drops tiny weights
  expect_equal(nrow(df), length(fit$unit_weights))
  expect_equal(sort(df$weight), sort(as.numeric(fit$unit_weights)))
})

test_that("plot_data(type='weights') top_n keeps the largest weights", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  df  <- plot_data(fit, type = "weights", top_n = 2)
  expect_equal(nrow(df), 2L)
  expect_equal(df$weight, sort(as.numeric(fit$unit_weights),
                               decreasing = TRUE)[1:2])
  expect_error(plot_data(fit, type = "weights", top_n = 0), "top_n")
})

test_that("plot_data(type='weights') for SDID adds an omega/lambda panel", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  df  <- plot_data(fit, type = "weights")
  expect_true("panel" %in% names(df))
  expect_setequal(unique(df$panel), c("omega", "lambda"))
  expect_equal(df$weight[df$panel == "lambda"], as.numeric(fit$time_weights))
})

test_that("plot_data(type='pred_weights') returns one row per predictor", {
  fit <- scm_fit(
    y ~ d | id + time, data = panel_cov, method = "scm",
    predictors = list(pred(c("cov1", "cov2"), 1:10))
  )
  df <- plot_data(fit, type = "pred_weights")
  expect_identical(names(df), c("predictor", "weight"))
  expect_equal(df$predictor, names(fit$v_weights))
  expect_equal(df$weight, as.numeric(fit$v_weights))
})

test_that("plot_data errors on fits without the requested component", {
  fit_mc <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  expect_error(plot_data(fit_mc, type = "pred_weights"), "V")
  expect_error(plot_data(fit_mc, type = "weights"), "No unit weights")

  fit_stag <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_error(plot_data(fit_stag, type = "trend"), "per cohort")
})

test_that("plot_data has no method for unsupported objects", {
  expect_error(plot_data(1:10), "no method")
})

test_that("plot_data.scm_placebo returns tidy gaps and ratios frames", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)

  g <- plot_data(inf, type = "gaps")
  expect_identical(names(g), c("time", "gap", "unit", "series"))
  expect_setequal(unique(g$series), c("Treated", "Placebo (donor pool)"))
  expect_equal(g$gap[g$series == "Treated"], as.numeric(inf$treated_gap))

  r <- plot_data(inf, type = "ratios")
  expect_identical(names(r), c("unit", "ratio", "series"))
  expect_equal(r$ratio, as.numeric(inf$mspe_ratios_all))
  expect_equal(r$series[1], "Treated")
})

test_that("plot_data.scm_placebo gaps honors mspe_prune", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  n_all  <- length(unique(plot_data(inf, type = "gaps")$unit))
  n_keep <- length(unique(plot_data(inf, type = "gaps", mspe_prune = 1)$unit))
  expect_lte(n_keep, n_all)
  expect_error(plot_data(inf, type = "gaps", mspe_prune = -1), "positive")
})

# ── Phase 33: Partially pooled staggered SCM (Ben-Michael et al. 2022) ────────

test_that("Phase 33: nu/fixedeff on a sharp fit error informatively", {
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm", nu = 0.5),
    regexp = "staggered"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm", fixedeff = TRUE),
    regexp = "staggered"
  )
})

test_that("Phase 33: invalid nu is rejected", {
  expect_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm", nu = 1.5),
    regexp = "nu must be"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm", nu = "bogus"),
    regexp = "nu must be"
  )
})

test_that("Phase 33: nu is incompatible with the legacy-path knobs", {
  expect_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm",
            nu = 0.5, lambda_pen = 0.1),
    regexp = "cannot be combined"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm",
            nu = 0.5, v_selection = "oos"),
    regexp = "cannot be combined"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = staggered, method = "scm",
            nu = 0.5, donor_mspe_threshold = 20),
    regexp = "cannot be combined"
  )
})

test_that("Phase 33: partially pooled SCM returns a valid staggered fit", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm", nu = 0.5)
  expect_s3_class(fit, "coresynth_staggered")
  expect_true(is.finite(fit$estimate))
  expect_lt(abs(fit$estimate - 1.5), 4.0)
  expect_equal(sum(fit$cohort_estimates$weight), 1, tolerance = 1e-8)
  for (cf in fit$cohort_fits) {
    expect_equal(sum(cf$unit_weights), 1, tolerance = 1e-5)
    expect_true(all(cf$unit_weights >= -1e-5))
  }
})

test_that("Phase 33: pooling diagnostics are stored and consistent", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm", nu = 0.5)
  p <- fit$pooling
  expect_equal(p$nu, 0.5)
  expect_true(p$fixedeff == FALSE)
  expect_true(all(is.finite(c(p$q_sep, p$q_pool,
                              p$q_sep_separate, p$q_pool_separate))))
  # pooling can only improve the pooled fit relative to separate SCM
  expect_lte(p$q_pool, p$q_pool_separate + 1e-10)
  # and can only worsen the separate fit
  expect_gte(p$q_sep, p$q_sep_separate - 1e-10)
})

test_that("Phase 33: nu = 0 equals the separate (uniform-V) solution", {
  fit0 <- scm_fit(y ~ d | id + time, data = staggered, method = "scm", nu = 0)
  p <- fit0$pooling
  expect_equal(p$q_sep, p$q_sep_separate, tolerance = 1e-12)
  expect_equal(p$q_pool, p$q_pool_separate, tolerance = 1e-12)
})

test_that("Phase 33: pooled fit improves monotonically in nu", {
  fit0 <- scm_fit(y ~ d | id + time, data = staggered, method = "scm", nu = 0)
  fit1 <- scm_fit(y ~ d | id + time, data = staggered, method = "scm", nu = 1)
  expect_lte(fit1$pooling$q_pool, fit0$pooling$q_pool + 1e-10)
  expect_gte(fit1$pooling$q_sep,  fit0$pooling$q_sep - 1e-10)
})

test_that("Phase 33: nu = 'auto' picks the heuristic in [0, 1]", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm",
                 nu = "auto")
  expect_true(fit$pooling$nu >= 0 && fit$pooling$nu <= 1)
  expect_equal(fit$pooling$nu, fit$pooling$nu_heuristic)
  expect_true(is.finite(fit$estimate))
})

test_that("Phase 33: fixedeff works on both staggered paths", {
  # legacy path
  fit_l <- scm_fit(y ~ d | id + time, data = staggered, method = "scm",
                   fixedeff = TRUE)
  expect_true(is.finite(fit_l$estimate))
  expect_true(fit_l$fixedeff)
  # pooled path
  fit_p <- scm_fit(y ~ d | id + time, data = staggered, method = "scm",
                   nu = 0.5, fixedeff = TRUE)
  expect_true(is.finite(fit_p$estimate))
  expect_true(fit_p$pooling$fixedeff)
  # Y_synth is reported on the raw outcome scale: pre-period levels track
  # the treated series (alpha restores the intercept)
  for (cf in fit_p$cohort_fits) {
    pre <- seq_len(cf$T_pre)
    expect_lt(abs(mean(cf$Y_treat[pre]) - mean(cf$Y_synth[pre])), 1.0)
  }
})

test_that("Phase 33: fixedeff recovers ATT under large unit intercepts", {
  # DGP where donors share the treated units' trend but sit at shifted
  # levels: plain SCM cannot match the level, the intercept shift can.
  set.seed(11)
  N <- 12; TT <- 24
  f <- cumsum(rnorm(TT, 0, 0.5))
  shift <- c(0, 0.5, seq(4, 14, length.out = N - 2))
  rows <- expand.grid(time = seq_len(TT), id = paste0("u", seq_len(N)),
                      stringsAsFactors = FALSE)
  rows$y <- rep(f, N) + rep(shift, each = TT) + rnorm(nrow(rows), 0, 0.2)
  rows$d <- as.integer(
    (rows$id == "u1" & rows$time >= 9) | (rows$id == "u2" & rows$time >= 15)
  )
  rows$y[rows$d == 1] <- rows$y[rows$d == 1] + 1.5

  fit_fe <- scm_fit(y ~ d | id + time, data = rows, method = "scm",
                    nu = 0.5, fixedeff = TRUE)
  expect_lt(abs(fit_fe$estimate - 1.5), 0.5)
})

test_that("Phase 33: staggered cohort_fits expose idx_tr/idx_co/Y_treat_mat", {
  for (nu in list(NULL, 0.5)) {
    fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm",
                   nu = nu)
    for (cf in fit$cohort_fits) {
      expect_true(is.integer(cf$idx_tr) || is.numeric(cf$idx_tr))
      expect_true(length(cf$idx_co) >= 2L)
      expect_true(is.matrix(cf$Y_treat_mat))
      expect_equal(ncol(cf$Y_treat_mat), cf$n_treated)
      expect_equal(nrow(cf$Y_treat_mat), 24L)
    }
  }
})

# ── Phase 33: scm_inference() wild bootstrap ─────────────────────────────────

test_that("Phase 33: scm_inference errors on sharp and non-SCM fits", {
  fit_sharp <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(scm_inference(fit_sharp), regexp = "staggered")
  fit_sdid <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  expect_error(scm_inference(fit_sdid), regexp = "method = 'scm'")
})

test_that("Phase 33: scm_inference returns a valid coresynth_inference", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  inf <- suppressWarnings(scm_inference(fit, n_boot = 200, seed = 1))
  expect_s3_class(inf, "coresynth_inference")
  expect_equal(inf$estimate, fit$estimate)
  expect_true(is.finite(inf$se) && inf$se > 0)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  expect_lt(inf$ci_lower, inf$ci_upper)
  expect_equal(length(inf$boot_ests), 200L)
  expect_true(inf$staggered)
  expect_equal(inf$method, "wild_bootstrap")
  expect_equal(inf$n_treated, 2L)
})

test_that("Phase 33: scm_inference warns with fewer than 5 treated units", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  expect_warning(scm_inference(fit, n_boot = 50, seed = 1),
                 regexp = "fewer than 5 treated units")
})

test_that("Phase 33: scm_inference is reproducible with seed and works with tidy/glance", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm",
                 nu = 0.5, fixedeff = TRUE)
  inf1 <- suppressWarnings(scm_inference(fit, n_boot = 100, seed = 42))
  inf2 <- suppressWarnings(scm_inference(fit, n_boot = 100, seed = 42))
  expect_equal(inf1$boot_ests, inf2$boot_ests)

  td <- broom::tidy(inf1)
  expect_equal(nrow(td), 1L)
  expect_equal(td$method, "wild_bootstrap")
  expect_true(td$staggered)
  gl <- broom::glance(inf1)
  expect_equal(gl$n_boot_valid, 100L)
})

test_that("Phase 33: scm_inference alternative directions behave sensibly", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "scm")
  # true ATT = 1.5 > 0: 'greater' should not have a larger p than 'less'
  p_g <- suppressWarnings(scm_inference(fit, n_boot = 200, seed = 3,
                                        alternative = "greater"))$p_value
  p_l <- suppressWarnings(scm_inference(fit, n_boot = 200, seed = 3,
                                        alternative = "less"))$p_value
  expect_lte(p_g, p_l)
})

test_that("Phase 33: partially pooled SCM works with covariates partial-out", {
  fit <- scm_fit(y ~ d | id + time, data = stag_cov, method = "scm",
                 nu = 0.5, covariates = "cov1")
  expect_true(is.finite(fit$estimate))
  expect_length(fit$beta_hat, 1L)
})

# ── Phase 34: plot align / show_donors / SDID weight panels / placebo SE ─────

test_that("Phase 34: plot(align=TRUE) makes the SDID post-period mean gap equal the estimate", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  d   <- ggplot2::ggplot_build(plot(fit, type = "gap", align = TRUE))$data[[1]]
  expect_equal(mean(d$y[d$x > fit$T_pre]), fit$estimate, tolerance = 1e-8)
  # the raw gap plot lacks this property (SDID has a free intercept)
  d0 <- ggplot2::ggplot_build(plot(fit, type = "gap"))$data[[1]]
  expect_false(isTRUE(all.equal(mean(d0$y[d0$x > fit$T_pre]), fit$estimate,
                                tolerance = 1e-8)))
})

test_that("Phase 34: plot(align=TRUE) shifts only the synthetic series, by the lambda-weighted pre-gap", {
  fit  <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  pre  <- seq_len(fit$T_pre)
  off  <- sum(fit$time_weights * (fit$Y_treat[pre] - fit$Y_synth[pre]))
  base <- ggplot2::ggplot_build(plot(fit, type = "trend"))$data[[1]]
  al   <- ggplot2::ggplot_build(plot(fit, type = "trend", align = TRUE))$data[[1]]
  shift <- al$y - base$y
  expect_true(all(abs(shift) < 1e-10 | abs(shift - off) < 1e-10))
  expect_true(any(abs(shift - off) < 1e-10))  # synthetic series moved
  expect_true(any(abs(shift) < 1e-10))        # treated series did not
})

test_that("Phase 34: plot(align=TRUE) centers the pre-period gap at zero for non-SDID fits", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  d   <- ggplot2::ggplot_build(plot(fit, type = "gap", align = TRUE))$data[[1]]
  expect_equal(mean(d$y[d$x <= fit$T_pre]), 0, tolerance = 1e-10)
})

test_that("Phase 34: plot align argument is validated", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(plot(fit, type = "trend", align = "yes"), "align")
  expect_error(plot(fit, type = "gap", align = 1), "align")
})

test_that("Phase 34: show_donors overlays top-weight donor paths", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  p0  <- plot(fit, type = "trend")
  p2  <- plot(fit, type = "trend", show_donors = 2)
  expect_equal(length(p2$layers), length(p0$layers) + 1L)
  don <- ggplot2::ggplot_build(p2)$data[[1]]  # donor layer is drawn first
  expect_equal(length(unique(don$group)), 2L)
  don_all <- ggplot2::ggplot_build(
    plot(fit, type = "trend", show_donors = Inf))$data[[1]]
  expect_equal(length(unique(don_all$group)), length(fit$unit_weights))
})

test_that("Phase 34: Donors series accepts color/label overrides only when shown", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  p <- plot(fit, type = "trend", show_donors = 3,
            colors = c(donors = "black"), labels = c(donors = "Donor pool"))
  expect_s3_class(p, "ggplot")
  expect_error(plot(fit, type = "trend", colors = c(donors = "black")),
               "unrecognized")
})

test_that("Phase 34: show_donors validation and unsupported methods", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(plot(fit, type = "trend", show_donors = -1), "show_donors")
  expect_error(plot(fit, type = "trend", show_donors = c(1, 2)), "show_donors")
  fit_mc <- scm_fit(y ~ d | id + time, data = panel, method = "mc")
  expect_error(plot(fit_mc, type = "trend", show_donors = 2), "show_donors")
})

test_that("Phase 34: SDID weights plot shows unit and time weight panels", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  p <- plot(fit, type = "weights")
  expect_s3_class(p, "ggplot")
  built <- ggplot2::ggplot_build(p)
  expect_equal(length(unique(built$data[[1]]$PANEL)), 2L)
  n_units <- sum(fit$unit_weights > 1e-4)
  n_times <- sum(fit$time_weights > 1e-4)
  expect_equal(nrow(built$data[[1]]), n_units + n_times)
  # top_n prunes only the unit panel
  built2 <- ggplot2::ggplot_build(plot(fit, type = "weights", top_n = 2))
  expect_equal(nrow(built2$data[[1]]), 2L + n_times)
})

test_that("Phase 34: non-SDID weights plot keeps the single unit panel", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  built <- ggplot2::ggplot_build(plot(fit, type = "weights"))
  expect_equal(length(unique(built$data[[1]]$PANEL)), 1L)
})

test_that("Phase 34: sdid_inference placebo returns SE and normal CI (Clarke Alg. 4)", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  eff <- inf$placebo_effects
  expect_equal(inf$se, sqrt(mean((eff - mean(eff))^2)), tolerance = 1e-12)
  expect_true(inf$se > 0)
  expect_equal(inf$ci_lower, inf$estimate - qnorm(0.975) * inf$se,
               tolerance = 1e-12)
  expect_equal(inf$ci_upper, inf$estimate + qnorm(0.975) * inf$se,
               tolerance = 1e-12)
  inf80 <- sdid_inference(fit, method = "placebo", level = 0.80)
  expect_equal(inf80$ci_upper, inf80$estimate + qnorm(0.90) * inf80$se,
               tolerance = 1e-12)
  td <- broom::tidy(inf)
  expect_true(is.finite(td$std.error))
  expect_true(is.finite(td$conf.low) && is.finite(td$conf.high))
})

test_that("Phase 34: sdid_inference staggered placebo also returns SE and CI", {
  fit <- scm_fit(y ~ d | id + time, data = staggered, method = "sdid")
  inf <- sdid_inference(fit, method = "placebo")
  eff <- inf$placebo_effects
  expect_equal(inf$se, sqrt(mean((eff - mean(eff))^2)), tolerance = 1e-12)
  expect_true(is.finite(inf$ci_lower) && is.finite(inf$ci_upper))
  expect_lt(inf$ci_lower, inf$ci_upper)
})

# ── Phase 35: vline positioning (vline_offset / vline xintercept) ─────────────

# xintercept values of every geom_vline layer in a built plot
.vline_x <- function(p) {
  d  <- ggplot2::ggplot_build(p)$data
  xs <- lapply(d, function(l) l[["xintercept"]])
  as.numeric(unlist(xs[!vapply(xs, is.null, logical(1))]))
}

test_that("plot.coresynth: vline_offset moves the treatment line", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  # panel: times 1:20, T_pre = 10, first post-treatment period = 11
  expect_equal(.vline_x(plot(fit, type = "trend")), 11)
  expect_equal(.vline_x(plot(fit, type = "trend", vline_offset = -1)), 10)
  # fractional offsets interpolate between adjacent observed times
  expect_equal(.vline_x(plot(fit, type = "trend", vline_offset = -0.5)), 10.5)
  expect_equal(.vline_x(plot(fit, type = "gap", vline_offset = -1)), 10)
})

test_that("plot.coresynth: vline_offset is validated and range-checked", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_error(plot(fit, type = "trend", vline_offset = "a"), "single finite")
  expect_error(plot(fit, type = "trend", vline_offset = c(-1, 0)), "single finite")
  expect_error(plot(fit, type = "trend", vline_offset = Inf), "single finite")
  # out of range: warn and drop the line (layer count shrinks by one)
  n_default <- length(plot(fit, type = "trend")$layers)
  expect_warning(p <- plot(fit, type = "trend", vline_offset = 99),
                 "outside the observed time range")
  expect_equal(length(p$layers), n_default - 1L)
})

test_that("plot.coresynth: vline xintercept pins the line to an absolute time", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_equal(.vline_x(plot(fit, type = "trend",
                             vline = list(xintercept = 15.5))), 15.5)
  # aesthetics in the same list still apply
  p <- plot(fit, type = "gap", vline = list(xintercept = 5, color = "red"))
  expect_equal(.vline_x(p), 5)
  # several positions draw several lines
  expect_equal(.vline_x(plot(fit, type = "trend",
                             vline = list(xintercept = c(5, 11)))), c(5, 11))
  expect_error(plot(fit, type = "trend", vline_offset = -1,
                    vline = list(xintercept = 5)), "not both")
  expect_error(plot(fit, type = "trend", vline = list(xintercept = NA)),
               "non-missing")
})

test_that("plot.coresynth: vline_offset works on a Date time axis", {
  dpanel <- panel
  dpanel$time <- as.Date("2000-01-01") + 30 * (dpanel$time - 1L)
  fit <- scm_fit(y ~ d | id + time, data = dpanel, method = "scm")

  d0 <- as.numeric(as.Date("2000-01-01") + 30 * 10) # 11th period
  expect_equal(.vline_x(plot(fit, type = "trend")), d0)
  expect_equal(.vline_x(plot(fit, type = "trend", vline_offset = -1)), d0 - 30)
  expect_equal(.vline_x(plot(fit, type = "trend", vline_offset = -0.5)), d0 - 15)
  # absolute positions accept a date string on a Date axis
  expect_equal(
    .vline_x(plot(fit, type = "trend",
                  vline = list(xintercept = "2000-03-01"))),
    as.numeric(as.Date("2000-03-01"))
  )
})

test_that("plot.scm_placebo gaps: vline_offset and xintercept work", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  inf <- mspe_ratio_pval(fit)
  expect_equal(.vline_x(plot(inf, type = "gaps")), 11)
  expect_equal(.vline_x(plot(inf, type = "gaps", vline_offset = -1)), 10)
  expect_equal(.vline_x(plot(inf, type = "gaps",
                             vline = list(xintercept = 4))), 4)
  expect_error(plot(inf, type = "gaps", vline_offset = 2.5,
                    vline = list(xintercept = 4)), "not both")
})

# ── Phase 36: multi-start outer V optimisation & v_window ─────────────────────

test_that("Phase 36: predictor-path default resolves to multistart", {
  fit <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                 predictors = list(pred(c("cov1", "cov2"), 1:10)))
  expect_identical(fit$v_optim_effective, "multistart")
})

test_that("Phase 36: outcomes-only default keeps single-start coord descent", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_identical(fit$v_optim_effective, "coord_descent")
  fit_cd <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                    v_optim = "coord_descent")
  expect_equal(fit$unit_weights, fit_cd$unit_weights, tolerance = 1e-12)
})

test_that("Phase 36: multistart is never worse than coord_descent in pre-loss", {
  spec <- list(pred(c("cov1", "cov2"), 1:10))
  fit_ms <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                    predictors = spec, v_optim = "multistart")
  fit_cd <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                    predictors = spec, v_optim = "coord_descent")
  expect_lte(fit_ms$loss, fit_cd$loss + 1e-8)
})

test_that("Phase 36: multistart is deterministic and leaves the R RNG alone", {
  spec <- list(pred(c("cov1", "cov2"), 1:10), pred("y", c(3, 6, 9)))
  set.seed(123)
  seed_before <- get(".Random.seed", envir = globalenv())
  fit1 <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                  predictors = spec, v_optim = "multistart")
  expect_identical(get(".Random.seed", envir = globalenv()), seed_before)
  fit2 <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                  predictors = spec, v_optim = "multistart")
  expect_identical(fit1$unit_weights, fit2$unit_weights)
  expect_identical(fit1$v_weights, fit2$v_weights)
})

test_that("Phase 36: multistart errors on outcomes-only fits", {
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            v_optim = "multistart"),
    regexp = "predictor"
  )
})

test_that("Phase 36: mspe_ratio_pval mirrors the multistart fit", {
  fit <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                 predictors = list(pred(c("cov1", "cov2"), 1:10)))
  inf <- mspe_ratio_pval(fit)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
  expect_equal(length(inf$mspe_ratios_all), 10L)
})

test_that("Phase 36: v_window equal to the full pre-period reproduces default", {
  fit_def <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  fit_win <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                     v_window = 1:10)
  expect_equal(fit_def$unit_weights, fit_win$unit_weights, tolerance = 1e-12)
  expect_equal(fit_win$v_window, fit_def$times[1:10])
  expect_equal(fit_win$z_rows, 1:10)
})

test_that("Phase 36: v_window subset restricts the outer evaluation", {
  fit_win <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                     v_window = 4:10)
  expect_equal(fit_win$z_rows, 4:10)
  expect_equal(sum(fit_win$unit_weights), 1, tolerance = 1e-4)
  expect_true(is.finite(fit_win$estimate))
  # loss is still reported on the full pre-treatment window
  gap_pre <- (fit_win$Y_treat - fit_win$Y_synth)[seq_len(fit_win$T_pre)]
  expect_equal(fit_win$loss, sqrt(sum(gap_pre^2)), tolerance = 1e-8)
})

test_that("Phase 36: v_window works with predictors and mspe_ratio_pval", {
  fit <- scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
                 predictors = list(pred(c("cov1", "cov2"), 1:10)),
                 v_window = 5:10)
  expect_equal(fit$z_rows, 5:10)
  inf <- mspe_ratio_pval(fit)
  expect_true(inf$p_value >= 0 && inf$p_value <= 1)
})

test_that("Phase 36: v_window validation errors", {
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "sdid", v_window = 1:5),
    regexp = "method"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm", v_window = 5:12),
    regexp = "pre-treatment"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm", v_window = 5),
    regexp = "at least 2"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            v_window = 1:8, v_selection = "oos"),
    regexp = "oos"
  )
})

test_that("Phase 36: v_window errors on staggered fits", {
  set.seed(99)
  sp <- expand.grid(id = paste0("u", 1:8), time = 1:12)
  sp$d <- as.integer((sp$id == "u1" & sp$time >= 7) |
                     (sp$id == "u2" & sp$time >= 9))
  sp$y <- rnorm(nrow(sp)) + as.numeric(sub("u", "", sp$id)) + 0.3 * sp$time
  expect_error(
    scm_fit(y ~ d | id + time, data = sp, method = "scm", v_window = 1:6),
    regexp = "sharp"
  )
})

# ── Phase 38: outcomes-only cheap bordered-KKT inner-QP face solve ─────────────
# The outcomes-only path solves each active-set face with a cheap bordered-KKT
# direct solve (cheap_face = TRUE) rather than the scale-robust null-space
# projection used for user predictors. It is the same nested V/W estimator, so
# the two solvers must agree to round-off; the flag only changes speed.

# Build outcome matrices with a low-rank factor structure; optionally make two
# donor pairs near-collinear to stress the inner face solver.
.mk_oo <- function(seed, N_co, T_pre, sd, degen = FALSE) {
  set.seed(seed)
  r  <- 3L
  L  <- matrix(rnorm((N_co + 1L) * r), N_co + 1L, r)
  Ft <- matrix(rnorm(T_pre * r), T_pre, r)
  Y  <- Ft %*% t(L) +
    matrix(rnorm(T_pre * (N_co + 1L), sd = sd), T_pre, N_co + 1L)
  if (degen) {
    Y[, 3] <- Y[, 2] + rnorm(T_pre, sd = 1e-4)
    Y[, 5] <- Y[, 4] + rnorm(T_pre, sd = 1e-4)
  }
  list(X1 = Y[, 1], X0 = Y[, -1, drop = FALSE])
}

test_that("Phase 38: cheap_face matches the null-space solver (outcomes-only)", {
  grid <- expand.grid(seed = 1:5, sd = c(0.1, 0.3, 1.0),
                      degen = c(FALSE, TRUE))
  for (i in seq_len(nrow(grid))) {
    g <- grid[i, ]
    d <- .mk_oo(g$seed, N_co = 25L, T_pre = 18L, sd = g$sd, degen = g$degen)
    cheap  <- scm_weights_cpp(d$X0, d$X1, d$X0, d$X1, cheap_face = TRUE)
    robust <- scm_weights_cpp(d$X0, d$X1, d$X0, d$X1, cheap_face = FALSE)
    # Same estimator, different inner face solve: agree to round-off.
    expect_equal(as.numeric(cheap$W), as.numeric(robust$W), tolerance = 1e-8)
    expect_equal(cheap$loss, robust$loss, tolerance = 1e-10)
    # Simplex feasibility of the fast-path solution.
    expect_true(all(cheap$W >= -1e-10))
    expect_equal(sum(cheap$W), 1, tolerance = 1e-8)
  }
})

test_that("Phase 38: scm_fit outcomes-only equals the robust null-space solver", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_equal(fit$v_optim_effective, "coord_descent")
  Yp <- fit$Y_co_pre
  yt <- fit$Y_treat[seq_len(fit$T_pre)]
  # The fit uses cheap_face internally; the null-space solver is the reference.
  ref <- scm_weights_cpp(Yp, yt, Yp, yt, cheap_face = FALSE)
  expect_equal(unname(fit$unit_weights), as.numeric(ref$W), tolerance = 1e-8)
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-8)
})

test_that("Phase 38: near-collinear donor pool stays feasible (outcomes-only)", {
  d  <- .mk_oo(123L, N_co = 20L, T_pre = 15L, sd = 0.05, degen = TRUE)
  df <- data.frame(
    id   = rep(seq_len(ncol(d$X0) + 1L), each = nrow(d$X0) + 4L),
    time = rep(seq_len(nrow(d$X0) + 4L), times = ncol(d$X0) + 1L)
  )
  Yfull <- cbind(d$X1, d$X0)
  Ypost <- Yfull[rep(nrow(Yfull), 4L), ] + 0.1  # 4 post periods
  Yall  <- rbind(Yfull, Ypost)
  df$y  <- as.vector(Yall)
  df$d  <- as.integer(df$id == 1L & df$time > nrow(d$X0))
  fit <- scm_fit(y ~ d | id + time, data = df, method = "scm")
  expect_true(all(fit$unit_weights >= -1e-10))
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-8)
  expect_true(is.finite(fit$estimate))
})

# ── outcomes-only routing of per-period outcome predictor specs ───────────────
# A predictors list naming solely the outcome at each single pre-treatment
# period defines X0/X1 equal to the pre-treatment outcome rows, so it must
# return the identical outcomes-only fit rather than entering the scaled
# multistart predictor path.

test_that("per-period outcome spec returns the identical outcomes-only fit", {
  spec <- lapply(1:10, function(t) pred("y", t))
  ref  <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_message(
    fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                   predictors = spec),
    regexp = "outcomes-only"
  )
  expect_equal(fit$v_optim_effective, "coord_descent")
  expect_null(fit$X0_mat)
  expect_null(fit$predictor_table)
  expect_identical(fit$unit_weights, ref$unit_weights)
  expect_identical(fit$v_weights, ref$v_weights)
  expect_identical(fit$estimate, ref$estimate)
})

test_that("outcomes-only routing is order- and op-invariant", {
  spec <- lapply(rev(1:10), function(t) pred("y", t, op = "median"))
  ref  <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  expect_message(
    fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                   predictors = spec),
    regexp = "outcomes-only"
  )
  expect_identical(fit$unit_weights, ref$unit_weights)
  expect_identical(fit$estimate, ref$estimate)
})

test_that("partial-coverage or windowed outcome specs keep the predictor path", {
  fit_sub <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                     predictors = lapply(1:9, function(t) pred("y", t)))
  expect_equal(fit_sub$v_optim_effective, "multistart")
  expect_false(is.null(fit_sub$X0_mat))

  spec_win <- c(lapply(1:8, function(t) pred("y", t)), list(pred("y", 9:10)))
  fit_win  <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                      predictors = spec_win)
  expect_equal(fit_win$v_optim_effective, "multistart")
  expect_false(is.null(fit_win$X0_mat))
})

test_that("explicit v_optim = 'multistart' overrides the outcomes-only routing", {
  spec <- lapply(1:10, function(t) pred("y", t))
  fit  <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                  predictors = spec, v_optim = "multistart")
  expect_equal(fit$v_optim_effective, "multistart")
  expect_false(is.null(fit$X0_mat))
  expect_equal(sum(fit$unit_weights), 1, tolerance = 1e-8)
})

# Coverage is judged against the pre-treatment periods `data` carries, so the
# same predictors list lands on either path depending on how far back `data`
# reaches. That switch changes the estimand, so it has to be announced.

test_that("a partial-coverage outcome-lag spec reports the missed routing", {
  spec <- lapply(3:10, function(t) pred("y", t))
  expect_message(
    fit_full <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                        predictors = spec),
    regexp = "covering 8 of the 10 pre-treatment periods"
  )
  expect_equal(fit_full$v_optim_effective, "multistart")
  expect_false(is.null(fit_full$X0_mat))

  # Same spec, `data` trimmed to the periods it covers: now it is the
  # outcomes-only fit, on the same call.
  panel_trim <- panel[panel$time >= 3, ]
  expect_message(
    fit_trim <- scm_fit(y ~ d | id + time, data = panel_trim, method = "scm",
                        predictors = spec),
    regexp = "outcomes-only"
  )
  expect_equal(fit_trim$v_optim_effective, "coord_descent")
  expect_null(fit_trim$X0_mat)
})

test_that("the missed-routing message names v_window only when it is unset", {
  spec <- lapply(3:10, function(t) pred("y", t))
  expect_message(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            predictors = spec),
    regexp = "v_window"
  )
  msgs <- capture_messages(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            predictors = spec, v_window = 3:10)
  )
  expect_true(any(grepl("covering 8 of the 10", msgs)))
  expect_false(any(grepl("v_window", msgs)))
})

test_that("specs that are not single-period outcome lags stay silent", {
  msgs <- capture_messages(
    scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
            predictors = list(pred(c("cov1", "cov2"), 1:10)))
  )
  expect_false(any(grepl("pre-treatment periods in 'data'", msgs)))

  msgs_ms <- capture_messages(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            predictors = lapply(3:10, function(t) pred("y", t)),
            v_optim = "multistart")
  )
  expect_false(any(grepl("pre-treatment periods in 'data'", msgs_ms)))
})

# pred() windows that do not aggregate what they read as ──────────────────────

test_that("pred() windows reaching past the treatment date warn", {
  expect_warning(
    scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
            predictors = list(pred("cov1", 1:10), pred("cov2", 8:14))),
    regexp = "post-treatment"
  )
})

test_that("pred() windows only partly present in the panel warn", {
  expect_warning(
    scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
            predictors = list(pred("cov1", -3:10))),
    regexp = "does not contain"
  )
})

test_that("a pred() window with no overlap errors instead of warning", {
  expect_silent(suppressMessages(expect_error(
    scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
            predictors = list(pred("cov1", 100:110))),
    regexp = "non-finite"
  )))
})

test_that("windows that lie inside the pre-treatment period do not warn", {
  expect_no_warning(
    scm_fit(y ~ d | id + time, data = panel_cov, method = "scm",
            predictors = list(pred(c("cov1", "cov2"), 1:10),
                              pred("y", c(3, 6, 9))))
  )
})

# ── Inner-QP solver: exactness and the Wolfe (Caratheodory-sparse) option ─────

# KKT residual of W for the inner QP min ||X1 - X0 w||_V^2 over the simplex,
# normalised by the gradient scale.
scm_kkt_residual <- function(X0, X1, V, W) {
  g    <- as.numeric(crossprod(X0, V * (as.numeric(X0 %*% W) - X1)))
  supp <- W > 1e-10
  mu   <- -mean(g[supp])
  viol <- max(abs(g[supp] + mu))
  if (!all(supp)) viol <- max(viol, 0, -min(g[!supp] + mu))
  viol / (1 + max(abs(g)))
}

test_that("scm_inner_weights_cpp returns a KKT-exact optimum", {
  set.seed(4242L)
  for (cfg in list(c(k = 1, N = 40), c(k = 3, N = 60), c(k = 12, N = 25),
                   c(k = 30, N = 20))) {
    k <- cfg[["k"]]; N <- cfg[["N"]]
    X0 <- matrix(rnorm(k * N), k, N) + 5
    X1 <- rowMeans(X0) + rnorm(k, 0, 0.2)
    V  <- runif(k); V <- V / sum(V)
    W  <- as.numeric(scm_inner_weights_cpp(X0, X1, V))
    expect_equal(sum(W), 1, tolerance = 1e-8)
    expect_true(all(W >= -1e-10))
    expect_lt(scm_kkt_residual(X0, X1, V, W), 1e-8)
  }
})

test_that("Wolfe solver matches the active set in objective but is sparse", {
  set.seed(99L)
  for (cfg in list(c(k = 1, N = 50), c(k = 2, N = 80), c(k = 5, N = 120))) {
    k <- cfg[["k"]]; N <- cfg[["N"]]
    X0 <- matrix(rnorm(k * N), k, N) + 5
    X1 <- rowMeans(X0) + rnorm(k, 0, 0.2)
    V  <- rep(1 / k, k)
    wa <- as.numeric(scm_inner_weights_cpp(X0, X1, V, wolfe = FALSE))
    ww <- as.numeric(scm_inner_weights_cpp(X0, X1, V, wolfe = TRUE))
    obj <- function(w) sum(V * (X1 - as.numeric(X0 %*% w))^2)
    expect_equal(obj(ww), obj(wa), tolerance = 1e-8)
    expect_lt(scm_kkt_residual(X0, X1, V, ww), 1e-8)
    # Caratheodory: an optimum supported on at most k+1 donors always exists
    expect_lte(sum(ww > 1e-8), k + 1L)
  }
})

test_that("qp_solver = 'wolfe' leaves a non-degenerate outcomes-only fit alone", {
  fit_a <- scm_fit(y ~ d | id + time, data = panel, method = "scm")
  fit_w <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                   qp_solver = "wolfe")
  expect_equal(fit_w$unit_weights, fit_a$unit_weights, tolerance = 1e-6)
  expect_equal(fit_w$estimate, fit_a$estimate, tolerance = 1e-6)
  expect_equal(fit_w$qp_solver, "wolfe")
  expect_equal(fit_a$qp_solver, "active_set")
})

test_that("qp_solver = 'wolfe' rejects the paths it is not wired through", {
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            qp_solver = "wolfe", v_selection = "oos"),
    "not available with"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "scm",
            qp_solver = "wolfe", lambda_pen = 0.1),
    "not available with"
  )
  expect_error(
    scm_fit(y ~ d | id + time, data = panel, method = "sdid",
            qp_solver = "wolfe"),
    "method = \"scm\" only"
  )
})

test_that("mspe_ratio_pval mirrors the fit's qp_solver", {
  fit <- scm_fit(y ~ d | id + time, data = panel, method = "scm",
                 qp_solver = "wolfe")
  res <- mspe_ratio_pval(fit)
  expect_s3_class(res, "scm_placebo")
  expect_true(is.finite(res$p_value))
  expect_gte(res$p_value, 0)
  expect_lte(res$p_value, 1)
})

Try the coresynth package in your browser

Any scripts or data that you put into this service are public.

coresynth documentation built on Aug. 28, 2026, 1:06 a.m.