test_that("tmpinv works for AP/TM on symmetric 20x20 matrix", {
  skip_if_not_installed("rclsp")
  
  set.seed(123456789)
  m <- 20L
  p <- 20L
  
  # --- create symmetric positive matrix ---
  X_true <- abs(matrix(rnorm(m * p), nrow = m, ncol = p))
  X_true <- 0.5 * (X_true + t(X_true))
  
  # 10% known values
  n_cells <- m * p
  idx <- sample.int(n_cells, size = max(1L, floor(0.1 * n_cells)),
                    replace = FALSE)
  M     <- diag(n_cells)[idx, , drop = FALSE]
  b_row <- rowSums(X_true)
  b_col <- colSums(X_true)
  b_val <- matrix(as.numeric(X_true)[idx], ncol = 1L)
  
  # run estimator
  result <- tmpinv(
    M = M, b_row = b_row, b_col = b_col, b_val = b_val,
    bounds = c(0, NA),
    symmetric = TRUE,
    r = 1L,
    alpha = 1.0
  )
  
  # --- tests ---
  expect_true(is.list(result))
  expect_true(is.matrix(result$x))
  expect_equal(dim(result$x), c(m, p))
  
  # model must exist
  expect_true(is.list(result$model))
  
  # stability metrics must be finite
  expect_true(is.finite(result$model$kappaC))
  expect_true(is.finite(result$model$kappaB))
  expect_true(is.finite(result$model$kappaA))
  
  # compute NRMSE
  expect_true(is.finite(result$model$nrmse))
  
  # confidence bands
  expect_true(all(is.finite(result$model$x_lower)))
  expect_true(all(is.finite(result$model$x_upper)))
  
  # bootstrap t-test (NRMSE)
  ttest_res <- rclsp::ttest(result$model,
    sample_size = 30L,
    seed = 123456789L,
    distribution = rnorm,
    partial = TRUE
  )
  
  expect_true(is.list(ttest_res))
  expect_true(length(ttest_res) > 0)
})

test_that("tmpinv works for reduced AP/TM with zero diagonal (40x40)", {
  skip_if_not_installed("rclsp")
  
  set.seed(123456789)
  m <- 40L
  p <- 40L
  
  # --- create zero-diagonal matrix ---
  X_true <- abs(matrix(rnorm(m * p), nrow = m, ncol = p))
  diag(X_true) <- 0
  
  # 20% known values
  n_cells <- m * p
  idx <- sample.int(n_cells, size = max(1L, floor(0.2 * n_cells)),
                    replace = FALSE)
  M     <- diag(n_cells)[idx, , drop = FALSE]
  b_row <- rowSums(X_true)
  b_col <- colSums(X_true)
  b_val <- matrix(as.numeric(X_true)[idx], ncol = 1L)
  
  # run estimator
  result <- tmpinv(
    M = M, b_row = b_row, b_col = b_col, b_val = b_val,
    zero_diagonal = TRUE,
    reduced = c(20L, 20L),
    bounds = c(0, NA),
    r = 1L,
    alpha = 1.0
  )
  
  # --- tests ---
  expect_true(is.list(result))
  expect_true(is.matrix(result$x))
  expect_equal(dim(result$x), c(m, p))
  
  # reduced model returns list of submodels
  expect_true(is.list(result$model))
  expect_gt(length(result$model), 0L)
  
  # each block model must be valid
  for (CLSP in result$model) {
    expect_true(is.finite(CLSP$kappaC))
    expect_true(is.finite(CLSP$kappaB))
    expect_true(is.finite(CLSP$kappaA))
    expect_true(is.finite(CLSP$nrmse))
    expect_true(all(is.finite(CLSP$x_lower)))
    expect_true(all(is.finite(CLSP$x_upper)))
  }
  
  # bootstrap t-test (NRMSE) (per block)
  for (CLSP in result$model) {
    ttest_res <- rclsp::ttest(CLSP, sample_size = 30L, seed = 123456789L,
                            distribution = rnorm, partial = TRUE)
    expect_true(is.list(ttest_res))
    expect_true(length(ttest_res) > 0)
  }
})