test_that("seair_solve conserves total mass", {
  m <- seq(0, 1, length.out = 41)
  f <- stats::dbeta(m, 2, 2)
  init <- seair_init(m, f, I_seed = 1e-3)
  pars <- seair_params(beta = 1.5, sigma = 0.3, kappa = 0.2,
                       gamma_A = 0.1, gamma_I = 0.13,
                       alpha = 0.6, delta = 0.4)
  sol <- seair_solve(init, pars, times = seq(0, 50, by = 5))
  agg <- seair_aggregate(sol)
  totals <- rowSums(agg[, c("S", "E", "A", "I", "R")])
  expect_true(all(abs(totals - 1) < 1e-6))
})

test_that("seair_solve agrees with baseline SIR in the degenerate limit", {
  ## When alpha = delta = 1 and kappa is small, the asymptomatic stage
  ## dominates and the model is close to SIR-like dynamics; here we
  ## only check that trajectories stay non-negative and bounded.
  m <- seq(0, 1, length.out = 41)
  f <- stats::dbeta(m, 5, 5)
  init <- seair_init(m, f, I_seed = 1e-4)
  pars <- seair_params(beta = 1.0, sigma = 1.0, kappa = 1e-3,
                       gamma_A = 0.15, gamma_I = 0.15,
                       alpha = 1, delta = 1)
  sol <- seair_solve(init, pars, times = seq(0, 60, by = 2))
  expect_true(all(sol$S >= -1e-8))
  expect_true(all(sol$I >= -1e-8))
  expect_true(all(sol$R >= -1e-8))
  expect_true(max(sol$I) <= 1 + 1e-6)
})