test_that("double_ridge produces valid numerical equivalence outputs", {
  set.seed(42)
  n_p <- 100
  d <- 4
  X_p <- matrix(rnorm(n_p * d), n_p, d)
  Y <- X_p %*% c(1, 2, -1, 0.5) + rnorm(n_p)
  target_mean <- colMeans(matrix(rnorm(50 * d, mean = 0.3), 50, d))

  fit <- double_ridge(Y, X_p, target_mean, lambda = 0.1, delta = 0.05)

  expect_true(is.numeric(fit$estimate))
  expect_false(is.na(fit$estimate))
  expect_false(is.nan(fit$estimate))
  expect_equal(length(fit$beta_aug), d)
  expect_equal(length(fit$weights), n_p)

  # Proposition 3.1 & 3.2 equivalence check: estimate == sum(target_mean * beta_aug)
  expect_equal(fit$estimate, sum(target_mean * fit$beta_aug), tolerance = 1e-6)

  # Check exact OLS collapse when delta = 0 and lambda = 0
  fit_ols <- double_ridge(Y, X_p, target_mean, lambda = 0, delta = 0)
  expect_equal(fit_ols$beta_aug, fit_ols$beta_ols, tolerance = 1e-5)
})