# Seber (1984, p. 460) y Amaro, Vicente-Villardon y Galindo (2004), tablas I a VII
cesped <- data.frame(
  trat = c("GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GC", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI", "GI"),
  temp = c("T14", "T14", "T14", "T14", "T14", "T14", "T18", "T18", "T18", "T18", "T18", "T18", "T22", "T22", "T22", "T22", "T22", "T22", "T26", "T26", "T26", "T26", "T26", "T26", "T14", "T14", "T14", "T14", "T14", "T14", "T18", "T18", "T18", "T18", "T18", "T18", "T22", "T22", "T22", "T22", "T22", "T22", "T26", "T26", "T26", "T26", "T26", "T26"),
  X1 = c(2.2, 3, 3.3, 2.2, 2, 3.5, 21.8, 11, 16.4, 13.1, 15.4, 14.5, 13.6, 6.2, 16.7, 12.2, 8.7, 12.3, 3, 5.3, 3.1, 4.8, 3.4, 7.4, 2.3, 3, 2.3, 2.5, 2.4, 2.7, 10.1, 7.6, 19.7, 4.3, 5.2, 3.9, 10, 12.3, 4.9, 9.6, 6.5, 13.6, 4.2, 2.2, 2.8, 1.3, 4.2, 3), X2 = c(23.5, 27, 24.5, 20.5, 19, 23.5, 41.5, 32.5, 46.5, 31, 41.5, 46, 29.5, 23.5, 58.5, 40.5, 37, 41.5, 24, 26.5, 24.5, 34, 22.5, 32, 23.5, 21, 22, 22.5, 21.5, 25, 43.5, 27, 32.5, 28.5, 33.5, 24.5, 21, 49, 28.5, 27, 29, 30.5, 25.5, 23.5, 19.5, 21.5, 28.5, 25), X3 = c(1.7, 2.3, 3.2, 1.5, 2, 2.9, 23, 15.4, 22.8, 21.5, 20.8, 20.3, 30.8, 14.6, 36, 23.9, 20.3, 27.7, 10.2, 15.6, 14.7, 20.5, 14.3, 23.2, 2, 2.7, 1.8, 2.4, 1.1, 2.6, 14.2, 14.7, 21.4, 9.7, 12.2, 8.2, 23.6, 28.1, 13.3, 24.6, 19.3, 31.5, 13.3, 8.5, 11.8, 7.8, 15.1, 11.8))
X <- log(cesped[, 3:5])
ajusta <- function(e) bigarx_manova2(X, cesped$trat, cesped$temp, efecto = e)

test_that("contrastes MANOVA de dos factores (tabla I)", {
  pr <- ajusta("total")$diseno$pruebas
  w <- pr$estadistico[pr$contraste == "Wilks"]
  expect_equal(w, c(0.765, 0.0126, 0.8017), tolerance = 1e-3)
  ref <- summary(manova(as.matrix(X) ~ trat * temp, data = cesped), test = "Wilks")$stats
  expect_equal(w, unname(ref[1:3, 2]), tolerance = 1e-10)
})

test_that("biplot total, columnas e interaccion reproducen las tablas II a VII", {
  inercia_t2 <- function(cb) { th <- cb$canonical$valor_propio; 100 * th^2 / sum(th^2) }
  tot <- ajusta("total"); col <- ajusta("columnas"); int <- ajusta("interaccion")
  expect_equal(inercia_t2(tot), c(98.3906, 1.6083, 0.0011), tolerance = 1e-4)
  expect_equal(inercia_t2(col)[1:2], c(98.8244, 1.1756), tolerance = 1e-4)
  expect_equal(inercia_t2(int)[1:2], c(79.1750, 20.8084), tolerance = 1e-4)
  expect_equal(unname(tot$quality_rows[c("GC-T14", "GC-T18"), 1]), c(99.6858, 8.0247), tolerance = 1e-4)
  expect_equal(unname(tot$quality_cols_principal[, 1]), c(77.0349, 77.5299, 98.4058), tolerance = 1e-4)
  expect_equal(unname(col$quality_rows[, 1]), c(99.6234, 13.2106, 92.0380, 72.4175), tolerance = 1e-4)
  expect_equal(unname(col$quality_cols_principal[, 1]), c(81.6214, 83.4358, 98.8201), tolerance = 1e-4)
  expect_equal(unname(int$quality_rows[c("GC-T14", "GC-T18"), 1]), c(3.1963, 94.3761), tolerance = 1e-4)
})

test_that("un efecto con 1 g.l. se proyecta sobre el biplot total", {
  f <- ajusta("filas")
  expect_true(f$diseno$proyectado)
  expect_equal(nrow(f$rows), 2)
})

test_that("T2 de Hotelling por pares coincide con el contraste exacto", {
  cb <- bigarx_canonical(iris[, 1:4], iris$Species)
  hp <- hotelling_pares(cb)
  W <- cb$.W; d <- colMeans(cb$.Z[iris$Species == "setosa", ]) - colMeans(cb$.Z[iris$Species == "versicolor", ])
  T2 <- 25 * drop(t(d) %*% solve(W, d))
  expect_equal(hp$T2[1], T2, tolerance = 1e-8)
  expect_equal(hp$gl2[1], 147 - 4 + 1)
})

test_that("datos composicionales: clr con inversa generalizada y alr", {
  set.seed(1)
  g <- rep(c("a", "b", "c"), each = 30)
  P <- t(sapply(seq_along(g), function(i) { x <- rgamma(4, shape = c(2, 3, 4, 5) + 2 * (g[i] == "b") * c(1, 0, 0, 0)); x / sum(x) }))
  colnames(P) <- c("p1", "p2", "p3", "p4")
  expect_true(es_composicional(P))
  expect_error(bigarx_canonical(transformar_composicional(P, "clr"), g), "singular")
  cl <- bigarx_canonical(transformar_composicional(P, "clr"), g, singular = "ginv")
  al <- bigarx_canonical(transformar_composicional(P, "alr", "p4"), g)
  expect_true(cl$singular)
  expect_equal(cl$.krz$p_eff, 3)
  # clr con inversa generalizada y alr dan el mismo MANOVA (invariancia)
  expect_equal(cl$canonical$valor_propio, al$canonical$valor_propio, tolerance = 1e-8)
})

test_that("el rango nunca supera K - 1 (eje espurio por redondeo)", {
  set.seed(2026)
  idx <- unlist(lapply(split(1:150, iris$Species), function(i) sample(i, 50, TRUE)))
  for (k in 1:200) {
    idx <- unlist(lapply(split(1:150, iris$Species), function(i) sample(i, 50, TRUE)))
    expect_lte(bigarx_canonical(iris[idx, 1:4], iris$Species[idx], fix_signs = FALSE)$rank, 2)
  }
})

test_that("invertir_ejes solo cambia el signo", {
  cb <- bigarx_canonical(iris[, 1:4], iris$Species)
  ci <- invertir_ejes(cb, 2)
  expect_equal(ci$rows[, 2], -cb$rows[, 2]); expect_equal(ci$rows[, 1], cb$rows[, 1])
  expect_equal(canonical_regions(ci)$semieje_mayor, canonical_regions(cb)$semieje_mayor)
})