# =============================================================================
# Tests for DACE package
# =============================================================================

# Test dace_manifold() function
# -----------------------------------------------------------------------------

test_that("dace_manifold creates valid manifold object", {
  # Create simple test data
  set.seed(123)
  features <- matrix(rnorm(50 * 10), nrow = 50, ncol = 10)
  rownames(features) <- paste0("obs_", 1:50)
  colnames(features) <- paste0("feat_", 1:10)
  
  # Create manifold
  manifold <- dace_manifold(features, dims = 3)
  
  # Check object structure
  expect_s3_class(manifold, "dace_manifold")
  expect_true(is.list(manifold))
  
  # Check required components
  expect_true("coordinates" %in% names(manifold))
  expect_true("model" %in% names(manifold))
  expect_true("variance_explained" %in% names(manifold))
  expect_true("dims" %in% names(manifold))
  
  # Check dimensions
  expect_equal(nrow(manifold$coordinates), 50)
  expect_equal(ncol(manifold$coordinates), 3)
  expect_equal(manifold$dims, 3)
})

test_that("dace_manifold handles scaling options correctly", {
  set.seed(456)
  features <- matrix(rnorm(30 * 5), nrow = 30, ncol = 5)
  
  # With scaling
  m_scaled <- dace_manifold(features, scale = TRUE, dims = 3)
  expect_true(m_scaled$scale)
  
  # Without scaling
  m_unscaled <- dace_manifold(features, scale = FALSE, dims = 3)
  expect_false(m_unscaled$scale)
})

test_that("dace_manifold validates input correctly", {
  # Invalid: not a matrix/data.frame
  expect_error(
    dace_manifold(c(1, 2, 3)),
    "must be a matrix or data.frame"
  )
  
  # Invalid: too few rows for PCA
  expect_error(
    dace_manifold(matrix(1:5, nrow = 1)),
    "must have at least 2 rows"
  )
  
  # Invalid: dims exceeds features
  features <- matrix(rnorm(20 * 3), nrow = 20, ncol = 3)
  expect_error(
    dace_manifold(features, dims = 5),
    "cannot exceed number of features"
  )
  
  # Invalid: non-finite values
  bad_features <- matrix(rnorm(20 * 5), nrow = 20, ncol = 5)
  bad_features[1, 1] <- NA
  expect_error(
    dace_manifold(bad_features),
    "non-finite values"
  )
})

test_that("dace_manifold assigns rownames if missing", {
  features <- matrix(rnorm(20 * 5), nrow = 20, ncol = 5)
  # No rownames provided
  
  manifold <- dace_manifold(features, dims = 3)
  expect_true(!is.null(rownames(manifold$coordinates)))
  expect_true(all(grepl("obs_", rownames(manifold$coordinates))))
})

test_that("dace_manifold variance explained is valid", {
  set.seed(789)
  features <- matrix(rnorm(40 * 8), nrow = 40, ncol = 8)
  
  manifold <- dace_manifold(features, dims = 3)
  
  expect_true(manifold$variance_explained >= 0)
  expect_true(manifold$variance_explained <= 1)
  expect_true(is.numeric(manifold$variance_explained))
})

# Test anchor_manifold() function
# -----------------------------------------------------------------------------

test_that("anchor_manifold creates valid anchored object", {
  set.seed(234)
  features <- matrix(rnorm(50 * 10), nrow = 50, ncol = 10)
  
  # Add signal for controls
  features[1:5, 1:5] <- features[1:5, 1:5] + 3    # POS
  features[6:10, 1:5] <- features[6:10, 1:5] - 3  # NEG
  
  rownames(features) <- paste0("obs_", 1:50)
  
  manifold <- dace_manifold(features, dims = 3)
  
  anchored <- anchor_manifold(
    manifold,
    pos_controls = paste0("obs_", 1:5),
    neg_controls = paste0("obs_", 6:10),
    scale = "unit",
    verbose = FALSE
  )
  
  # Check object structure
  expect_s3_class(anchored, "dace_manifold_anchored")
  expect_true(is.list(anchored))
  
  # Check required components
  expect_true("anchored_coordinates" %in% names(anchored))
  expect_true("rotation_matrix" %in% names(anchored))
  expect_true("scale_factor_x" %in% names(anchored))
  expect_true("pos_controls" %in% names(anchored))
  expect_true("neg_controls" %in% names(anchored))
})

test_that("anchor_manifold rotation matrix is orthonormal", {
  set.seed(345)
  features <- matrix(rnorm(50 * 10), nrow = 50, ncol = 10)
  features[1:5, 1:5] <- features[1:5, 1:5] + 2
  features[6:10, 1:5] <- features[6:10, 1:5] - 2
  rownames(features) <- paste0("obs_", 1:50)
  
  manifold <- dace_manifold(features, dims = 3)
  anchored <- anchor_manifold(
    manifold,
    pos_controls = paste0("obs_", 1:5),
    neg_controls = paste0("obs_", 6:10),
    verbose = FALSE
  )
  
  R <- anchored$rotation_matrix
  RtR <- t(R) %*% R
  
  # Check orthonormality: R^T R = I
  expect_equal(as.numeric(RtR), as.numeric(diag(3)), tolerance = 1e-6)
  
  # Check determinant = +1 (proper rotation)
  expect_equal(det(R), 1, tolerance = 1e-10)
})

test_that("anchor_manifold ensures POS mean X > NEG mean X", {
  set.seed(456)
  features <- matrix(rnorm(50 * 10), nrow = 50, ncol = 10)
  features[1:5, 1:5] <- features[1:5, 1:5] + 3
  features[6:10, 1:5] <- features[6:10, 1:5] - 3
  rownames(features) <- paste0("obs_", 1:50)
  
  manifold <- dace_manifold(features, dims = 3)
  anchored <- anchor_manifold(
    manifold,
    pos_controls = paste0("obs_", 1:5),
    neg_controls = paste0("obs_", 6:10),
    verbose = FALSE
  )
  
  pos_x <- mean(anchored$anchored_coordinates[paste0("obs_", 1:5), "AnchoredX"])
  neg_x <- mean(anchored$anchored_coordinates[paste0("obs_", 6:10), "AnchoredX"])
  
  expect_true(pos_x > neg_x)
})

test_that("anchor_manifold unit scaling works correctly", {
  set.seed(567)
  features <- matrix(rnorm(50 * 10), nrow = 50, ncol = 10)
  features[1:5, 1:5] <- features[1:5, 1:5] + 3
  features[6:10, 1:5] <- features[6:10, 1:5] - 3
  rownames(features) <- paste0("obs_", 1:50)
  
  manifold <- dace_manifold(features, dims = 3)
  
  # Unit scaling
  anchored_unit <- anchor_manifold(
    manifold,
    pos_controls = paste0("obs_", 1:5),
    neg_controls = paste0("obs_", 6:10),
    scale = "unit",
    verbose = FALSE
  )
  
  pos_x <- mean(anchored_unit$anchored_coordinates[paste0("obs_", 1:5), "AnchoredX"])
  neg_x <- mean(anchored_unit$anchored_coordinates[paste0("obs_", 6:10), "AnchoredX"])
  
  # With unit scaling, POS-NEG separation should be approximately 2
  expect_equal(pos_x - neg_x, 2, tolerance = 0.01)
  
  # No scaling
  anchored_none <- anchor_manifold(
    manifold,
    pos_controls = paste0("obs_", 1:5),
    neg_controls = paste0("obs_", 6:10),
    scale = "none",
    verbose = FALSE
  )
  
  expect_equal(anchored_none$scale_factor_x, 1)
})

test_that("anchor_manifold validates control IDs", {
  set.seed(678)
  features <- matrix(rnorm(30 * 5), nrow = 30, ncol = 5)
  rownames(features) <- paste0("obs_", 1:30)
  
  manifold <- dace_manifold(features, dims = 3)
  
  # Missing positive controls
  expect_error(
    anchor_manifold(
      manifold,
      pos_controls = c("missing_1", "missing_2"),
      neg_controls = paste0("obs_", 6:10),
      verbose = FALSE
    ),
    "Positive controls not found"
  )
  
  # Missing negative controls
  expect_error(
    anchor_manifold(
      manifold,
      pos_controls = paste0("obs_", 1:5),
      neg_controls = c("missing_1", "missing_2"),
      verbose = FALSE
    ),
    "Negative controls not found"
  )
})

test_that("anchor_manifold requires 3D coordinates", {
  set.seed(789)
  features <- matrix(rnorm(30 * 5), nrow = 30, ncol = 5)
  rownames(features) <- paste0("obs_", 1:30)
  
  # Create 2D manifold (should fail anchoring)
  manifold_2d <- dace_manifold(features, dims = 2)
  
  expect_error(
    anchor_manifold(
      manifold_2d,
      pos_controls = paste0("obs_", 1:5),
      neg_controls = paste0("obs_", 6:10),
      verbose = FALSE
    ),
    "assumes 3D coordinates"
  )
})

test_that("anchor_manifold accepts plain coordinate matrix", {
  set.seed(890)
  
  # Create coordinate matrix directly
  coords <- matrix(rnorm(30 * 3), nrow = 30, ncol = 3)
  rownames(coords) <- paste0("obs_", 1:30)
  
  # Add signal
  coords[1:5, 1] <- coords[1:5, 1] + 3
  coords[6:10, 1] <- coords[6:10, 1] - 3
  
  # Should work with plain matrix
  anchored <- anchor_manifold(
    coords,
    pos_controls = paste0("obs_", 1:5),
    neg_controls = paste0("obs_", 6:10),
    verbose = FALSE
  )
  
  expect_s3_class(anchored, "dace_manifold_anchored")
})

# Test S3 methods
# -----------------------------------------------------------------------------

test_that("print methods work without errors", {
  set.seed(111)
  features <- matrix(rnorm(30 * 5), nrow = 30, ncol = 5)
  features[1:5, 1:3] <- features[1:5, 1:3] + 2
  features[6:10, 1:3] <- features[6:10, 1:3] - 2
  rownames(features) <- paste0("obs_", 1:30)
  
  manifold <- dace_manifold(features, dims = 3)
  
  expect_output(print(manifold), "DACE Manifold")
  expect_output(summary(manifold), "Top 5 loadings")
  
  anchored <- anchor_manifold(
    manifold,
    pos_controls = paste0("obs_", 1:5),
    neg_controls = paste0("obs_", 6:10),
    verbose = FALSE
  )
  
  expect_output(print(anchored), "DACE Anchored Manifold")
  expect_output(summary(anchored), "Effect axis")
})

# Test edge cases
# -----------------------------------------------------------------------------

test_that("anchor_manifold handles identical controls gracefully", {
  set.seed(222)
  features <- matrix(rnorm(30 * 5), nrow = 30, ncol = 5)
  rownames(features) <- paste0("obs_", 1:30)
  
  manifold <- dace_manifold(features, dims = 3)
  
  # Make POS and NEG identical (should fail)
  manifold$coordinates[1:5, ] <- manifold$coordinates[6:10, ]
  
  expect_error(
    anchor_manifold(
      manifold,
      pos_controls = paste0("obs_", 1:5),
      neg_controls = paste0("obs_", 6:10),
      verbose = FALSE
    ),
    "numerically"
  )
})

test_that("anchor_manifold handles single control per group", {
  set.seed(333)
  features <- matrix(rnorm(30 * 5), nrow = 30, ncol = 5)
  features[1, 1:3] <- features[1, 1:3] + 3
  features[2, 1:3] <- features[2, 1:3] - 3
  rownames(features) <- paste0("obs_", 1:30)
  
  manifold <- dace_manifold(features, dims = 3)
  
  # Should work with single controls
  anchored <- anchor_manifold(
    manifold,
    pos_controls = "obs_1",
    neg_controls = "obs_2",
    verbose = FALSE
  )
  
  expect_s3_class(anchored, "dace_manifold_anchored")
})