test_that("Laplace mechanism is unbiased", {
  # The expected value of noisy output should equal true output
  set.seed(42)
  true_value <- 100
  sensitivity <- 10
  epsilon <- 0.5
  n_samples <- 10000

  noisy_values <- replicate(n_samples, {
    tidydp:::add_laplace_noise(true_value, sensitivity, epsilon)
  })

  mean_noisy <- mean(noisy_values)

  # With 10000 samples, mean should be very close to true value
  # Allow for 3 standard errors
  se <- sd(noisy_values) / sqrt(n_samples)
  expect_true(abs(mean_noisy - true_value) < 3 * se)
})

test_that("Laplace mechanism has correct scale", {
  # The scale of Laplace noise should be sensitivity / epsilon
  set.seed(42)
  true_value <- 0  # Center at 0 for easier analysis
  sensitivity <- 10
  epsilon <- 0.5
  expected_scale <- sensitivity / epsilon
  n_samples <- 10000

  noisy_values <- replicate(n_samples, {
    tidydp:::add_laplace_noise(true_value, sensitivity, epsilon)
  })

  # Laplace distribution has variance = 2 * scale^2
  expected_var <- 2 * expected_scale^2
  observed_var <- var(noisy_values)

  # Allow 10% tolerance due to sampling variation
  expect_true(abs(observed_var - expected_var) / expected_var < 0.1)
})

test_that("Gaussian mechanism is unbiased", {
  # The expected value of noisy output should equal true output
  set.seed(42)
  true_value <- 100
  sensitivity <- 10
  epsilon <- 0.5
  delta <- 1e-5
  n_samples <- 10000

  noisy_values <- replicate(n_samples, {
    tidydp:::add_gaussian_noise(true_value, sensitivity, epsilon, delta)
  })

  mean_noisy <- mean(noisy_values)

  # With 10000 samples, mean should be very close to true value
  se <- sd(noisy_values) / sqrt(n_samples)
  expect_true(abs(mean_noisy - true_value) < 3 * se)
})

test_that("Gaussian mechanism has correct variance", {
  # The variance should match the theoretical formula
  set.seed(42)
  true_value <- 0  # Center at 0 for easier analysis
  sensitivity <- 10
  epsilon <- 0.5
  delta <- 1e-5
  n_samples <- 10000

  expected_sigma <- sensitivity * sqrt(2 * log(1.25 / delta)) / epsilon
  expected_var <- expected_sigma^2

  noisy_values <- replicate(n_samples, {
    tidydp:::add_gaussian_noise(true_value, sensitivity, epsilon, delta)
  })

  observed_var <- var(noisy_values)

  # Allow 10% tolerance due to sampling variation
  expect_true(abs(observed_var - expected_var) / expected_var < 0.1)
})

test_that("Laplace mechanism respects epsilon parameter", {
  # Smaller epsilon should result in more noise (larger variance)
  set.seed(42)
  true_value <- 100
  sensitivity <- 10
  epsilon_small <- 0.1
  epsilon_large <- 1.0
  n_samples <- 5000

  noisy_small <- replicate(n_samples, {
    tidydp:::add_laplace_noise(true_value, sensitivity, epsilon_small)
  })

  noisy_large <- replicate(n_samples, {
    tidydp:::add_laplace_noise(true_value, sensitivity, epsilon_large)
  })

  var_small <- var(noisy_small)
  var_large <- var(noisy_large)

  # Variance should be larger for smaller epsilon
  expect_true(var_small > var_large)

  # Check ratio matches theory: var(small) / var(large) = (eps_large / eps_small)^2
  expected_ratio <- (epsilon_large / epsilon_small)^2
  observed_ratio <- var_small / var_large
  expect_true(abs(observed_ratio - expected_ratio) / expected_ratio < 0.2)
})

test_that("Gaussian mechanism respects epsilon parameter", {
  # Smaller epsilon should result in more noise
  set.seed(42)
  true_value <- 100
  sensitivity <- 10
  epsilon_small <- 0.1
  epsilon_large <- 1.0
  delta <- 1e-5
  n_samples <- 5000

  noisy_small <- replicate(n_samples, {
    tidydp:::add_gaussian_noise(true_value, sensitivity, epsilon_small, delta)
  })

  noisy_large <- replicate(n_samples, {
    tidydp:::add_gaussian_noise(true_value, sensitivity, epsilon_large, delta)
  })

  var_small <- var(noisy_small)
  var_large <- var(noisy_large)

  # Variance should be larger for smaller epsilon
  expect_true(var_small > var_large)

  # Check ratio matches theory
  expected_ratio <- (epsilon_large / epsilon_small)^2
  observed_ratio <- var_small / var_large
  expect_true(abs(observed_ratio - expected_ratio) / expected_ratio < 0.2)
})

test_that("Sensitivity calculations are correct", {
  # Test count sensitivity
  expect_equal(tidydp:::sensitivity_count(), 1)

  # Test sum sensitivity (substitution model)
  expect_equal(tidydp:::sensitivity_sum(0, 10), 10)
  expect_equal(tidydp:::sensitivity_sum(-5, 5), 10)
  expect_equal(tidydp:::sensitivity_sum(-10, 0), 10)
  expect_equal(tidydp:::sensitivity_sum(-3, 7), 10)

  # Test mean sensitivity
  expect_equal(tidydp:::sensitivity_mean(0, 10, 100), 0.1)
  expect_equal(tidydp:::sensitivity_mean(-5, 5, 50), 0.2)
  expect_equal(tidydp:::sensitivity_mean(0, 1, 10), 0.1)
})

test_that("Laplace mechanism handles vector inputs correctly", {
  set.seed(42)
  values <- c(10, 20, 30, 40, 50)
  sensitivity <- 5
  epsilon <- 0.5

  noisy_values <- tidydp:::add_laplace_noise(values, sensitivity, epsilon)

  # Should return same length
  expect_equal(length(noisy_values), length(values))

  # All values should be different from original (with extremely high probability)
  expect_false(all(noisy_values == values))

  # Each value should be reasonably close to original
  # (within 10 * scale with high probability)
  scale <- sensitivity / epsilon
  expect_true(all(abs(noisy_values - values) < 10 * scale))
})

test_that("Gaussian mechanism handles vector inputs correctly", {
  set.seed(42)
  values <- c(10, 20, 30, 40, 50)
  sensitivity <- 5
  epsilon <- 0.5
  delta <- 1e-5

  noisy_values <- tidydp:::add_gaussian_noise(values, sensitivity, epsilon, delta)

  # Should return same length
  expect_equal(length(noisy_values), length(values))

  # All values should be different from original
  expect_false(all(noisy_values == values))

  # Each value should be reasonably close to original
  sigma <- sensitivity * sqrt(2 * log(1.25 / delta)) / epsilon
  expect_true(all(abs(noisy_values - values) < 6 * sigma))
})

test_that("Mechanisms reject invalid parameters", {
  # Test Laplace mechanism
  expect_error(tidydp:::add_laplace_noise(100, 10, -0.5), "epsilon must be positive")
  expect_error(tidydp:::add_laplace_noise(100, 10, 0), "epsilon must be positive")
  expect_error(tidydp:::add_laplace_noise(100, -10, 0.5), "sensitivity must be positive")
  expect_error(tidydp:::add_laplace_noise(100, 0, 0.5), "sensitivity must be positive")

  # Test Gaussian mechanism
  expect_error(tidydp:::add_gaussian_noise(100, 10, -0.5, 1e-5), "epsilon must be positive")
  expect_error(tidydp:::add_gaussian_noise(100, 10, 0.5, -1e-5), "delta must be between 0 and 1")
  expect_error(tidydp:::add_gaussian_noise(100, 10, 0.5, 1.5), "delta must be between 0 and 1")
  expect_error(tidydp:::add_gaussian_noise(100, -10, 0.5, 1e-5), "sensitivity must be positive")
})

test_that("Laplace distribution shape is correct", {
  # Test that the noise follows Laplace distribution shape
  set.seed(42)
  true_value <- 0
  sensitivity <- 10
  epsilon <- 1.0
  scale <- sensitivity / epsilon
  n_samples <- 10000

  noisy_values <- replicate(n_samples, {
    tidydp:::add_laplace_noise(true_value, sensitivity, epsilon)
  })

  noise <- noisy_values - true_value

  # Test median is close to 0 (Laplace is symmetric)
  expect_true(abs(median(noise)) < 0.5)

  # Test mean absolute deviation
  # For Laplace(0, b), E[|X|] = b
  mad <- mean(abs(noise))
  expect_true(abs(mad - scale) / scale < 0.1)
})

test_that("Composition increases noise correctly", {
  # Running mechanism twice should result in sqrt(2) times the standard deviation
  # (for Gaussian) or sqrt(2) times scale (for variance)
  set.seed(42)
  true_value <- 100
  sensitivity <- 10
  epsilon <- 1.0
  delta <- 1e-5
  n_samples <- 5000

  # Single query
  single_noisy <- replicate(n_samples, {
    tidydp:::add_gaussian_noise(true_value, sensitivity, epsilon, delta)
  })

  # Two independent queries (basic composition: epsilon doubles, so sigma halves per query)
  # But we're adding two independent noises
  double_noisy <- replicate(n_samples, {
    noise1 <- tidydp:::add_gaussian_noise(true_value, sensitivity, epsilon/2, delta/2)
    noise2 <- tidydp:::add_gaussian_noise(noise1, sensitivity, epsilon/2, delta/2)
    noise2
  })

  # For same total epsilon, adding twice should roughly double variance
  # (This is basic composition)
  var_single <- var(single_noisy)
  var_double <- var(double_noisy)

  # var_double should be approximately 4 * var_single with basic composition
  # (since each query uses epsilon/2, so 4x variance, and we add them)
  # Actually, let me recalculate: if we split epsilon in half for two queries,
  # each has variance (2 * sensitivity / epsilon)^2 = 4 * (sensitivity/epsilon)^2
  # And adding two: 4 * var + 4 * var = 8 * var
  # Hmm, this is getting complicated. Let me just verify double is larger.

  expect_true(var_double > var_single)
})