X <- as.matrix(USArrests)
Z <- scale(X)
full <- function(method, ...) bigarx(X, method = method, ndim = 4, ...)

test_that("GH y JK reproducen exactamente la matriz con todos los ejes (Gabriel, 1971)", {
  for (m in c("gh", "jk", "sqrt")) {
    f <- full(m)
    expect_equal(unname(f$rows %*% t(f$cols)), unname(Z[, ]), tolerance = 1e-10,
                 ignore_attr = TRUE)
  }
  f <- bigarx(X, alpha = 0.3, ndim = 4)
  expect_equal(unname(f$rows %*% t(f$cols)), unname(Z[, ]), tolerance = 1e-10,
               ignore_attr = TRUE)
})

test_that("HJ: filas del JK y columnas del GH (Galindo, 1986)", {
  hj <- full("hj"); gh <- full("gh"); jk <- full("jk")
  expect_equal(hj$rows, jk$rows)
  expect_equal(hj$cols, gh$cols)
})

test_that("JK: las filas son las puntuaciones de componentes principales", {
  p <- prcomp(X, scale. = TRUE)
  jk <- full("jk")
  expect_equal(abs(unname(jk$rows)), abs(unname(p$x)), tolerance = 1e-10)
})

test_that("GH: el coseno entre flechas es la correlaci\u00F3n (rango completo)", {
  gh <- full("gh")
  H <- gh$cols
  cosang <- (H %*% t(H)) / sqrt(outer(rowSums(H^2), rowSums(H^2)))
  expect_equal(unname(cosang), unname(cor(X)), tolerance = 1e-10)
})

test_that("GH con pca_scaling: longitud de flecha = desviaci\u00F3n t\u00EDpica", {
  gh <- bigarx(X, method = "gh", transform = "center", ndim = 4, pca_scaling = TRUE)
  expect_equal(unname(sqrt(rowSums(gh$cols^2))), unname(apply(X, 2, sd)), tolerance = 1e-10)
})

test_that("las calidades de cada elemento suman 100 sobre todos los ejes", {
  f <- full("hj")
  expect_equal(unname(rowSums(f$quality_rows)), rep(100, nrow(X)))
  expect_equal(unname(rowSums(f$quality_cols)), rep(100, ncol(X)))
  expect_equal(unname(colSums(f$contrib_rows)), rep(100, 4))
})

test_that("bondad de ajuste seg\u00FAn Gabriel (1971) y Galindo (1986)", {
  d <- svd(Z)$d
  r <- length(d)
  hj <- full("hj"); gh <- full("gh"); jk <- full("jk")
  p4 <- 100 * sum(d[1:2]^4) / sum(d^4)
  p2 <- 100 * sum(d[1:2]^2) / sum(d^2)
  expect_equal(unname(hj$gof), c(p2, p4, p4))
  expect_equal(unname(gh$gof), c(p2, 100 * 2 / r, p4))
  expect_equal(unname(jk$gof), c(p2, p4, 100 * 2 / r))
})

test_that("la inercia suma 100", {
  f <- bigarx(X)
  expect_equal(sum(f$inertia$inercia), 100)
  expect_equal(f$inertia$acumulada[nrow(f$inertia)], 100)
})

test_that("predict() con los datos activos devuelve sus coordenadas", {
  for (m in c("hj", "gh", "jk")) {
    f <- full(m)
    expect_equal(unname(predict(f, X)), unname(f$rows), tolerance = 1e-10)
    expect_equal(unname(predict(f, X, type = "cols")), unname(f$cols), tolerance = 1e-10)
  }
  fd <- bigarx(X, transform = "double", ndim = 3)
  expect_equal(unname(predict(fd, X)), unname(fd$rows), tolerance = 1e-10)
})

test_that("filas suplementarias se a\u00F1aden sin cambiar los ejes", {
  f1 <- bigarx(X[-(1:5), ])
  f2 <- bigarx(X[-(1:5), ], sup_rows = X[1:5, ])
  expect_equal(f1$rows, f2$rows)
  expect_equal(dim(f2$sup_rows), c(5L, 3L))
})

test_that("as_bigarx() desde prcomp coincide con bigarx()", {
  p <- prcomp(X, scale. = TRUE)
  expect_equal(unname(as_bigarx(p, ndim = 4)$rows), unname(full("hj")$rows), tolerance = 1e-10)
})

test_that("mensajes de error claros", {
  expect_error(bigarx(iris), "no num\u00E9ricas")
  Xna <- X; Xna[1, 1] <- NA
  expect_error(bigarx(Xna), "NA")
  expect_error(bigarx(X, alpha = 2), "entre 0 y 1")
  expect_error(bigarx(X, groups = 1:3), "una entrada por fila")
  expect_error(summary(bigarx(X), axes = c(1, 5)), "ejes")
})

test_that("summary, print, plot, tidy y autoplot funcionan", {
  f <- bigarx(iris[, 1:4], groups = iris$Species)
  expect_s3_class(summary(f), "summary.bigarx")
  expect_output(print(f), "HJ-Biplot")
  expect_output(print(summary(f)), "Inercia")
  pdf(NULL); on.exit(dev.off())
  expect_invisible(plot(f))
  td <- bigarx_tidy(f)
  expect_equal(nrow(td), 150 + 4)
  skip_if_not_installed("ggplot2")
  expect_s3_class(ggplot2::autoplot(f), "ggplot")
})

test_that("plot() acepta titulo y etiquetas propias", {
  pdf(NULL); on.exit(dev.off())
  expect_invisible(plot(bigarx(USArrests), main = "Mi biplot", xlab = "A", ylab = "B", cex = 0.8))
})

test_that("min_quality acepta umbrales separados para filas y columnas", {
  skip_if_not_installed("ggplot2")
  f <- bigarx(USArrests, method = "jk")
  g1 <- ggplot2::autoplot(f, min_quality = 60)
  g2 <- ggplot2::autoplot(f, min_quality = c(60, 0))
  ncols <- function(g) sum(ggplot2::layer_data(g, 4)$xend != 0 | TRUE)
  expect_lt(nrow(ggplot2::layer_data(g1, 4)), nrow(ggplot2::layer_data(g2, 4)))
  expect_equal(nrow(ggplot2::layer_data(g2, 4)), 4L)
})