# ---- Imputacion ----
set.seed(1)
A <- matrix(rnorm(40 * 2), 40, 2) %*% matrix(rnorm(2 * 6), 2, 6)
A <- sweep(A, 2, c(10, 20, 5, 0, 3, 8), "+")
colnames(A) <- paste0("v", 1:6)

test_that("impute_svd recupera una matriz de rango 2 exacto", {
  Xna <- A
  idx <- cbind(c(2, 5, 9, 17, 30), c(1, 3, 6, 2, 4))
  Xna[idx] <- NA
  imp <- impute_svd(Xna, rank = 2, transform = "center", tol = 1e-10, maxit = 5000)
  expect_true(imp$converged)
  expect_equal(unname(imp$X[idx]), unname(A[idx]), tolerance = 1e-4)
  expect_equal(nrow(imp$values), 5L)
})

test_that("impute_svd no modifica los valores observados", {
  Xna <- as.matrix(USArrests); Xna[c(1, 7), 2] <- NA
  imp <- impute_svd(Xna)
  obs <- !is.na(Xna)
  expect_equal(imp$X[obs], Xna[obs])
  expect_false(anyNA(imp$X))
})

test_that("bigarx imputa con na_action = 'impute' y lo informa", {
  Xna <- as.matrix(USArrests); Xna[c(1, 7), 2] <- NA
  expect_error(bigarx(Xna), "impute")
  f <- bigarx(Xna, na_action = "impute")
  expect_s3_class(f$imputation, "bigarx_imputation")
  expect_output(print(f), "imputados")
})

test_that("impute_svd valida sus argumentos", {
  Xna <- as.matrix(USArrests); Xna[, 2] <- NA
  expect_error(impute_svd(Xna), "sin ning")
  expect_error(impute_svd(as.matrix(USArrests), rank = 10), "rank")
})

# ---- Biplot canonico ----
X <- as.matrix(iris[, 1:4]); g <- iris$Species
cb <- bigarx_canonical(X, g)

test_that("los valores propios coinciden con los de HE^-1 de manova()", {
  m <- summary(manova(X ~ g), test = "Wilks")
  expect_equal(unname(cb$canonical$valor_propio), unname(m$Eigenvalues[1, 1:2]), tolerance = 1e-8)
  expect_equal(cb$manova$wilks, unname(m$stats[1, "Wilks"]), tolerance = 1e-10)
})

test_that("los individuos son las funciones discriminantes de MASS::lda", {
  skip_if_not_installed("MASS")
  ld <- predict(MASS::lda(X, g))$x
  expect_equal(abs(unname(cb$individuals[, 1:2])), abs(unname(ld)), tolerance = 1e-8)
})

test_that("la varianza dentro de los grupos es 1 en cada eje canonico", {
  Ind <- cb$individuals[, 1:2]
  res <- Ind - cb$rows[as.integer(g), 1:2]
  W <- crossprod(res) / (nrow(X) - nlevels(g))
  expect_equal(unname(W), diag(2), tolerance = 1e-10)
})

test_that("las medias proyectadas son el promedio de los individuos", {
  expect_equal(unname(rowsum(cb$individuals, g) / 50), unname(cb$rows), tolerance = 1e-10)
})

test_that("P Q' reproduce las medias de grupo centradas", {
  M <- rowsum(scale(X, scale = FALSE), g) / 50
  expect_equal(unname(cb$rows %*% t(cb$cols)), unname(M), tolerance = 1e-10)
})

test_that("radios de confianza", {
  expect_equal(unname(cb$radius), rep(sqrt(qchisq(0.95, 2) / 50), 3))
  cbt <- bigarx_canonical(X, g, region = "t")
  expect_equal(unname(cbt$radius), rep(qt(0.975, 147) / sqrt(50), 3))
  cbb <- bigarx_canonical(X, g, region = "bonferroni")
  expect_true(all(cbb$radius > cbt$radius))
})

test_that("calidades de medias suman 100 en el espacio completo", {
  expect_equal(unname(rowSums(cb$quality_rows)), rep(100, 3))
})

test_that("estandarizar no cambia el analisis canonico", {
  cs <- bigarx_canonical(X, g, transform = "standardize")
  expect_equal(cs$rows, cb$rows, tolerance = 1e-10)
  expect_equal(cs$canonical, cb$canonical, tolerance = 1e-10)
})

test_that("predict proyecta nuevos individuos", {
  expect_equal(unname(predict(cb, X)), unname(cb$individuals), tolerance = 1e-10)
})

test_that("errores y salidas del biplot canonico", {
  expect_error(bigarx_canonical(X[1:100, ], droplevels(g[1:100])), "3 grupos")
  expect_error(bigarx_canonical(X, g[-1]), "una entrada")
  expect_output(print(cb), "Wilks")
  expect_output(print(summary(cb)), "correlaci")
  pdf(NULL); on.exit(dev.off())
  expect_invisible(plot(cb))
  td <- bigarx_tidy(cb)
  expect_true(all(c("media", "individuo", "columna") %in% td$tipo))
  skip_if_not_installed("ggplot2")
  expect_s3_class(ggplot2::autoplot(cb), "ggplot")
})

test_that("la version publicada (ellipse_amaro) reproduce Amaro et al. (2008), tablas 4 a 7", {
  ce <- bigarx_canonical(iris[, 1:4], iris$Species, region = "ellipse_amaro")
  rg <- canonical_regions(ce)
  # Tabla 4 (G1 y G3) y tabla 6 (G2: coeficiente 0.028; la tabla 4 trae la errata 0.2830)
  expect_equal(rg$var1, c(0.0215, 0.0213, 0.0214), tolerance = 0.02)
  expect_equal(rg$var2, c(0.1643, 0.028, 0.1034), tolerance = 0.02)
  expect_equal(abs(rg$cov), c(0.0042, 0.0034, 0.0075), tolerance = 0.05)
  # Tabla 7: semiejes, angulos y areas. Para G3 la tabla 7 trae 0.67, pero su
  # propia varianza (tabla 4: 0.1034) y su area (0.873 = pi * a * 0.35) dan 0.79.
  expect_equal(rg$semieje_mayor, c(0.99, 0.42, 0.79), tolerance = 0.03)
  expect_equal(rg$semieje_menor, c(0.35, 0.34, 0.35), tolerance = 0.03)
  expect_equal(rg$area, c(1.114, 0.418, 0.873), tolerance = 0.03)
  expect_equal(sort(abs(rg$angulo)), sort(c(1.7, 23.1, 5.21)), tolerance = 0.06)
})

test_that("las elipses se dibujan y aparecen en summary", {
  ce <- bigarx_canonical(iris[, 1:4], iris$Species, region = "ellipse")
  expect_output(print(summary(ce)), "Elipses")
  pdf(NULL); on.exit(dev.off())
  expect_invisible(plot(ce))
  skip_if_not_installed("ggplot2")
  expect_s3_class(ggplot2::autoplot(ce), "ggplot")
  cc <- canonical_regions(bigarx_canonical(iris[, 1:4], iris$Species))
  expect_equal(cc$semieje_mayor, cc$semieje_menor)
})


test_that("elipses de Krzanowski coherentes en escala W (revision del 6-10-2026)", {
  ce <- bigarx_canonical(iris[, 1:4], iris$Species, region = "ellipse")
  rg <- canonical_regions(ce)
  expect_equal(attr(rg, "escala"), "escala W (varianza dentro de los grupos = 1)")
  # centros = x$rows (escala W)
  expect_equal(rg$centro1, unname(ce$rows[, 1]))
  # theta = 1 + lambda/(K-1), lambda = autovalores de SCH W^-1 = d^2
  expect_equal(ce$.krz$theta, unname(1 + ce$d^2 / 2))
  expect_equal(ce$.krz$P, ce$rows)
  # formula de la ecuacion 3 calculada a mano para setosa, eje 1, todo en escala W
  P <- ce$rows; th <- ce$.krz$theta; g <- 2 / 147
  v11 <- 1 / 50 + (0.5 * P[1, 1]^2 + th[1] * (th[2] + g * th[1]) / (g * (th[2] - th[1])^2) * P[1, 2]^2) / 147
  expect_equal(rg$var1[1], unname(v11), tolerance = 1e-12)
  expect_equal(rg$var1, c(0.2174, 0.0375, 0.1368), tolerance = 2e-3)
  expect_equal(rg$var2, c(0.2939, 0.0376, 0.1791), tolerance = 2e-3)
  expect_equal(rg$semieje_mayor, c(1.3288, 0.5129, 1.0482), tolerance = 1e-3)
  expect_true(all(rg$definida_positiva))
  expect_true(all(rg$error_contorno < 1e-10))
  # los puntos dibujados cumplen la ecuacion del contorno
  S <- .krz_cov(ce, 1:2); q <- qchisq(0.95, 2)
  df <- .regions_df(ce, 1:2, grupo = rownames(ce$rows))
  for (i in 1:3) {
    z <- t(as.matrix(df[df$grupo == rownames(ce$rows)[i], c("x", "y")])) - ce$rows[i, 1:2]
    expect_lt(max(abs(colSums(z * solve(S[[i]], z)) - q)), 1e-10)
  }
  # el estimador alternativo de theta y la inversion de ejes
  cm <- bigarx_canonical(iris[, 1:4], iris$Species, region = "ellipse", theta_est = "ms")
  expect_equal(cm$.krz$theta, unname(ce$d^2 / 2))
  ci <- invertir_ejes(ce, 2)
  expect_equal(canonical_regions(ci)$semieje_mayor, rg$semieje_mayor)
  expect_equal(canonical_regions(ci)$cov, -rg$cov)
})