#Partially generated by Deepseek + checked/fixed by human
test_that("deriv_pif_p works correctly", {
  # Simple case with known values
  p     <- c(0.5, 0.5)
  p_cft <- c(0.4, 0.6)
  rr    <- c(2, 1)

  # Calculate expected values manually
  mu_obs   <- sum(p * rr) # 0.5*2 + 0.5*1 = 1.5
  mu_cft   <- sum(p_cft * rr) # 0.4*2 + 0.6*1 = 1.4
  expected <- (mu_cft / (mu_obs^2)) * (rr - 1) # (1.4/2.25) * (1, 0) ≈ (0.622, 0)

  # Test with pre-calculated mu values
  expect_equal(
    deriv_pif_p(p, p_cft, rr),
    expected
  )

  expect_equal(
    deriv_pif_p(p, p_cft, rr),
    expected
  )

  expect_equal(
    deriv_pif_p(p, p_cft, rr),
    expected
  )

  expect_equal(
    deriv_pif_p(p, p_cft, rr = rr),
    expected
  )

  # Test with NULL mu values (should calculate internally)
  expect_equal(
    deriv_pif_p(p, p_cft, rr),
    expected
  )

  # Edge case: all rr = 1 (should give 0 derivative)
  expect_equal(
    deriv_pif_p(p, p_cft, c(1, 1)),
    c(0, 0)
  )

  # Edge case: p = p_cft (should give 0 derivative when rr=1)
  expect_equal(
    deriv_pif_p(p, p, c(2, 1)),
    (sum(p*rr) / (sum(p*rr)^2)) * (c(2,1) - 1)
  )
})

test_that("deriv_pif_beta works correctly", {
  # Simple case with known values
  p        <- c(0.5, 0.5)
  p_cft    <- c(0.4, 0.6)
  rr       <- c(2, 1)
  rr_deriv <- c(0.5, 0.2) # arbitrary derivative values

  # Calculate expected values manually
  mu_obs    <- sum(p * rr) # 1.5
  mu_cft    <- sum(p_cft * rr) # 1.4
  numerator <- (mu_cft * p - mu_obs * p_cft)
  expected  <- (numerator / (mu_obs^2)) * rr_deriv

  # Test with pre-calculated mu values
  expect_equal(
    deriv_pif_beta(p, p_cft, rr, rr_deriv),
    expected
  )

  # Test with NULL mu values (should calculate internally)
  expect_equal(
    deriv_pif_beta(p, p_cft, rr, rr_deriv),
    expected
  )

  # Edge case: all rr_deriv = 0 (should give 0 derivative)
  expect_equal(
    deriv_pif_beta(p, p_cft, rr, c(0, 0)),
    c(0, 0)
  )

  # Edge case: p = p_cft
  result <- deriv_pif_beta(p, p, rr, rr_deriv)
  manual <- (sum(p*rr)*p - sum(p*rr)*p) / (sum(p*rr)^2) * rr_deriv
  expect_equal(result, manual)
})

test_that("Input validation works", {
    p        <- c(0.5, 0.5)
    p_cft    <- c(0.4, 0.6)
    rr       <- c(2, 1)
    rr_deriv <- c(0.5, 0.2)

    # Test for length mismatches
    expect_error(deriv_pif_p(p, p_cft, c(1, 2, 3)), "different lengths") # rr too long
    expect_error(deriv_pif_beta(p, p_cft, rr, c(1)), "different lengths") # rr_deriv too short

    # Test for invalid probabilities
    expect_error(deriv_pif_p(c(-0.1, 0.2), p_cft, rr),"Probabilities|probability")
    expect_error(deriv_pif_p(c(0.1, 1.2), p_cft, rr),"Probabilities|probability")
    expect_error(deriv_pif_p(p, c(-0.1, 0.2), rr),"Probabilities|probability")
    expect_error(deriv_pif_beta(c(-0.1, 0.2), p_cft, rr, rr_deriv), "Probabilities|probability")
    expect_error(deriv_pif_beta(p, c(-0.1, 0.2), rr, rr_deriv), "Probabilities|probability")
    expect_error(deriv_pif_beta(p, c(0.5, 0.6), rr, rr_deriv),"Probabilities|probability") # sums > 1

    # Test for invalid relative risks
    expect_error(deriv_pif_p(p, p_cft, c(0, 1)),"<= 0") # rr <= 0
    expect_error(deriv_pif_beta(p, p_cft, c(-1, 2), rr_deriv),"<= 0")

  })

  test_that("NA handling works as expected", {

    # Should error if critical values are NA
    expect_error(deriv_pif_p(c(0.5, NA), c(0.5, 0.5), c(1, 1)))
    expect_error(deriv_pif_p(c(0.5, 1.0), c(NA, 0.5), c(1, 1)))
    expect_error(deriv_pif_p(c(0.5, 0.5), c(0.5, 0.5), c(1, NA)))

    expect_error(deriv_pif_beta(c(0.5, 1.0), c(NA, 0.5), c(1, 1)))
    expect_error(deriv_pif_beta(c(NA, 1.0), c(0.5, 0.5), c(1, 1)))
    expect_error(deriv_pif_beta(c(0.5, 0.5), c(0.5, 0.5), c(1, 1), c(1, NA)))
  })