test_that("pqa validates input correctly", {
  # Test 1: non-table input should error
  expect_error(pqa(matrix(c(1,2,3,4,5,6), nrow=2, ncol=3)),
               "x must be an object of class `table'")

  # Test 2: one-dimensional table should error
  expect_error(pqa(table(c(1,1,2,2,3))),
               "Quetelet-Pearson Analysis is currently implemented only for two-way tables")

  # Test 3: three-dimensional table should error
  x <- table(c(1,1,2,2), c(1,2,1,2), c(1,1,2,2))
  expect_error(pqa(x),
               "Quetelet-Pearson Analysis is currently implemented only for two-way tables")
})

test_that("pqa processes valid table input correctly", {
  # Test 1: simple 2x2 table
  testdata <- as.table(matrix(c(10, 20, 15, 25), nrow=2, ncol=2))
  dimnames(testdata) <- list(Gender = c("M", "F"), Preference = c("A", "B"))
  result <- pqa(testdata)

  expect_s3_class(result, "pqa")
  expect_named(result, c("abs", "rel", "q", "pq", "apex", "chisq"))
  expect_equal(result$abs, testdata)

  # Test 2: table from vectors
  x <- c(1,1,1,2,2,2)
  y <- c(1,2,1,2,1,2)
  testdata <- table(x, y)
  result <- pqa(testdata)

  expect_s3_class(result, "pqa")
  expect_equal(result$abs, testdata)

  # Test 3: table from factors
  gender <- factor(c("M", "F", "M", "F", "M", "F"))
  preference <- factor(c("A", "B", "A", "B", "A", "B"))
  testdata <- table(gender, preference)
  result <- pqa(testdata)

  expect_s3_class(result, "pqa")
  expect_equal(result$abs, testdata)
})

test_that("pqa output structure and classes are correct", {
  testdata <- as.table(matrix(c(10, 20, 15, 25), nrow=2, ncol=2))
  dimnames(testdata) <- list(Gender = c("M", "F"), Preference = c("A", "B"))
  result <- pqa(testdata)

  # Test output is a list with pqa class
  expect_type(result, "list")
  expect_s3_class(result, "pqa")

  # Test all subtables have correct classes
  expect_s3_class(result$abs, "pqa.subtable.absfreq")
  expect_s3_class(result$rel, "pqa.subtable.relfreq")
  expect_s3_class(result$q, "pqa.subtable.quetelet_values")
  expect_s3_class(result$pq, "pqa.subtable.pq_decomposition")
  expect_s3_class(result$apex, "pqa.subtable.apex_tab")
  expect_s3_class(result$chisq, "pqa.chisq")

  # Test dimensions are preserved
  expect_equal(dim(result$abs), dim(testdata))
  expect_equal(dim(result$rel), dim(testdata))
  expect_equal(dim(result$q), dim(testdata))
  expect_equal(dim(result$pq), dim(testdata))
  expect_equal(dim(result$apex), dim(testdata))

  # Test chi-square structure
  expect_named(result$chisq, c("stat", "df", "flag", "pval", "fct_names"))
  expect_type(result$chisq$stat, "double")
  expect_type(result$chisq$df, "integer")
  expect_type(result$chisq$flag, "double")
  expect_type(result$chisq$pval, "double")
  expect_length(result$chisq$fct_names, 2)
})

test_that("pqa chi-square test calculations are correct", {
  # Test 1: Simple 2x2 table with known chi-square value
  testdata <- as.table(matrix(c(10, 20, 15, 25), nrow=2, ncol=2))
  dimnames(testdata) <- list(Gender = c("M", "F"), Preference = c("A", "B"))
  result <- pqa(testdata)

  # Calculate expected values manually
  tot_obs <- sum(testdata)
  expected_df <- (nrow(testdata) - 1) * (ncol(testdata) - 1)

  # Test degrees of freedom
  expect_equal(result$chisq$df, expected_df)
  expect_equal(result$chisq$df, 1)

  # Test chi-square statistic calculation: stat = tot_obs * phisq
  phisq <- sum(result$pq)
  expect_equal(result$chisq$stat, tot_obs * phisq)

  # Test p-value is valid probability
  expect_true(result$chisq$pval >= 0 && result$chisq$pval <= 1)
  expect_equal(result$chisq$pval, 1 - pchisq(result$chisq$stat, result$chisq$df))

  # Test factor names are captured
  expect_equal(result$chisq$fct_names, c("Gender", "Preference"))

  # Test 2: 3x3 table
  testdata2 <- as.table(matrix(c(5, 10, 15, 8, 12, 18, 7, 14, 21), nrow=3, ncol=3))
  dimnames(testdata2) <- list(Row = c("A", "B", "C"), Col = c("X", "Y", "Z"))
  result2 <- pqa(testdata2)

  expect_equal(result2$chisq$df, (3-1) * (3-1))
  expect_equal(result2$chisq$df, 4)
  expect_equal(result2$chisq$fct_names, c("Row", "Col"))

  # Test 3: 4x2 table
  testdata3 <- as.table(matrix(c(20, 30, 25, 35, 15, 20, 18, 22), nrow=4, ncol=2))
  result3 <- pqa(testdata3)

  expect_equal(result3$chisq$df, (4-1) * (2-1))
  expect_equal(result3$chisq$df, 3)
})

test_that("pqa chi-square test consistency with R's chisq.test", {
  # Test that our chi-square statistic matches R's built-in chisq.test
  testdata <- as.table(matrix(c(10, 20, 15, 25), nrow=2, ncol=2))
  dimnames(testdata) <- list(Gender = c("M", "F"), Preference = c("A", "B"))

  result <- pqa(testdata)
  r_chisq <- suppressWarnings(chisq.test(testdata, correct = FALSE))

  # Compare chi-square statistics
  expect_equal(result$chisq$stat, r_chisq$statistic[[1]], tolerance = 1e-10)

  # Compare degrees of freedom
  expect_equal(result$chisq$df, r_chisq$parameter[[1]])

  # Compare p-values
  expect_equal(result$chisq$pval, r_chisq$p.value, tolerance = 1e-10)
})

test_that("pqa chi-square validity_flag works correctly", {
  # Test 1: Valid test (all expected frequencies >= 5)
  testdata1 <- as.table(matrix(c(50, 60, 55, 65), nrow=2, ncol=2))
  result1 <- pqa(testdata1)
  expect_equal(result1$chisq$flag, 0)
  expect_false(is.na(result1$chisq$pval))

  # Test 2: Unreliable test (expected frequency >= 1 but < 5)
  testdata2 <- as.table(matrix(c(2, 3, 4, 5), nrow=2, ncol=2))
  result2 <- pqa(testdata2)
  expect_equal(result2$chisq$flag, 1)
  expect_false(is.na(result2$chisq$pval))

  # Test 3: Cannot compute (expected frequency < 1)
  testdata3 <- as.table(matrix(c(1, 0, 0, 1), nrow=2, ncol=2))
  result3 <- pqa(testdata3)
  expect_equal(result3$chisq$flag, 2)
  expect_true(is.na(result3$chisq$pval))
})

test_that("pqa handles perfectly independent tables without NaN apex", {
  testdata <- as.table(matrix(c(5, 5, 5, 5), nrow = 2, ncol = 2))
  dimnames(testdata) <- list(Row = c("A", "B"), Col = c("X", "Y"))
  result <- pqa(testdata)

  expect_equal(result$chisq$stat, 0)
  expect_true(all(is.finite(as.numeric(result$apex))))
  expect_true(all(as.numeric(result$apex) == 0))
  expect_no_error(capture.output(summary(result)))
})

test_that("pqa test against (orig) usa_voting_prefs", {
  testdata <- as.table(matrix(nrow=4, ncol=3, byrow = TRUE,
              data = c(2388, 2034, 4423, 2286,  938,
                2696, 3885, 1126, 3712, 3258,  695, 3049))
  )
  dimnames(testdata) <- list(
    Income=c("I", "II", "III", "IV"), Leaning=c("R", "U", "D")
  )

  q <- pqa(usa_voting_prefs)

  # Test 1: absolute frequencies (function should do nothing to the input)
  expect_equal(q[[1]], usa_voting_prefs)

  # Test 2: relative frequencies
  expect_equal(round(as.numeric(q[[2]]),4), c(0.0783, 0.0750, 0.1274, 0.1069,
    0.0667, 0.0308, 0.0369, 0.0228, 0.1451, 0.0884, 0.1217, 0.1000))

  # Test 3: test quetelet values
  expect_equal(round(as.numeric(q[[3]])*100,4), c(-30.3395, -0.3667,
    14.9146, 20.0547, 46.2861, 0.7932, -17.8850, -36.8588, 9.8467,
    0.0383, -6.5219, -4.3460))

  # Test 4: pearson-Quetelet decomposition
  expect_equal(round(sum(q[[4]]), 4), 0.0345)
  expect_equal(round(as.numeric(q[[4]]), 4), c(-0.0238,
   -0.0003, 0.0190, 0.0214, 0.0309, 0.0002, -0.0066, -0.0084, 0.0143,
    0.0000, -0.0079, -0.0043))

  # Test 5: apex
  expect_equal(round(as.numeric(q[[5]]*100),4), c(-68.7894,
    -0.7958, 55.0152, 62.0364, 89.3883, 0.7064, -19.1208, -24.3223,
    41.3511, 0.0979, -22.9858, -12.5812))

  # Test 6: chisq statistic
  expect_equal(q[[6]]$df, 3*2)
  expect_equal(round(q[[6]]$stat,3), 1053.224)
})