///|
/// A finite birth-death stationary law. Birth[k] moves k to k+1;
/// death[k] moves k+1 to k. All death rates must be strictly positive.
/// Zero birth rates make upper states transient; their mass is exactly zero.
pub fn stationary_birth_death(
  birth : Array[Double],
  death : Array[Double],
) -> Array[Double] raise QueueError {
  if birth.length() != death.length() || birth.length() > 10000 {
    raise Limit("birth/death shapes must match, at most 10000 transitions")
  }
  let logs : Array[Double] = [0.0]
  let mut active = true
  let mut value = 0.0
  let mut largest = 0.0
  for i in 0.. @math.exp(x - largest))
  let mut total = 0.0
  let mut compensation = 0.0
  for w in weights {
    let corrected = w - compensation
    let next = total + corrected
    compensation = next - total - corrected
    total = next
  }
  positive(total, "stationary normalizer")
  weights.map(w => w / total)
}

///|
/// Distribution summaries, with variance from central deviations (not E[N2]-E[N]^2).
pub(all) struct StateSummary {
  mean : Double
  variance : Double
  full_probability : Double
  empty_probability : Double
} derive(Debug, ToJson)

///|
fn validate_law(law : Array[Double]) -> Unit raise QueueError {
  if law.is_empty() || law.length() > 10001 {
    raise Limit("law needs 1..10001 states")
  }
  let mut sum = 0.0
  for p in law {
    probability(p, "state mass")
    sum += p
  }
  if (sum - 1.0).abs() > 1.0e-9 {
    raise Invalid("state masses must sum to one")
  }
}

///|
pub fn summarize_states(law : Array[Double]) -> StateSummary raise QueueError {
  validate_law(law)
  let mut mean = 0.0
  for i in 0..= p; p=0 returns 0.
pub fn state_quantile(law : Array[Double], p : Double) -> Int raise QueueError {
  validate_law(law)
  probability(p, "state quantile")
  let mut cumulative = 0.0
  for i in 0..= p {
      return i
    }
  }
  // Rounding at p=1; trailing zero-mass states must not change the answer.
  for offset in 0.. 0.0 {
      return i
    }
  }
  raise Numerical("empty probability support")
}

///|
/// Largest absolute detailed-balance residual, scaled by each edge's rate.
pub fn balance_residual(
  law : Array[Double],
  birth : Array[Double],
  death : Array[Double],
) -> Double raise QueueError {
  validate_law(law)
  if birth.length() + 1 != law.length() || death.length() != birth.length() {
    raise Invalid("detailed-balance shape")
  }
  let mut result = 0.0
  for i in 0.. 0.0 {
      result = result.max(
        (law[i] * (birth[i] / scale) - law[i + 1] * (death[i] / scale)).abs(),
      )
    }
  }
  result
}