///|
/// All rates use the same time unit. No runtime measurements are collected.
pub(all) struct Mmc {
  arrival : Double
  service : Double
  servers : Int
} derive(Debug, ToJson)

///|
/// Means are unconditional. Waiting excludes service; response includes it.
pub(all) struct QueueMetrics {
  utilization : Double
  delay_probability : Double
  mean_queue : Double
  mean_system : Double
  mean_wait : Double
  mean_response : Double
  throughput : Double
} derive(Debug, ToJson)

///|
fn Mmc::validate(self : Mmc) -> Unit raise QueueError {
  rate(self.arrival, "arrival")
  rate(self.service, "service")
  positive(self.service, "service")
  servers(self.servers)
  if self.arrival >= self.service * self.servers.to_double() {
    raise Unstable("arrival must be below total service capacity")
  }
}

///|
pub fn Mmc::analyze(self : Mmc) -> QueueMetrics raise QueueError {
  self.validate()
  let offered = self.arrival / self.service
  let utilization = offered / self.servers.to_double()
  let delay_probability = erlang_c(offered, self.servers)
  let spare = self.service * self.servers.to_double() - self.arrival
  let mean_wait = checked_ratio(delay_probability, spare, "mean wait")
  let mean_response = mean_wait + 1.0 / self.service
  let mean_queue = self.arrival * mean_wait
  let mean_system = mean_queue + offered
  for value in [mean_response, mean_queue, mean_system] {
    finite(value, "queue metric")
  }
  {
    utilization,
    delay_probability,
    mean_queue,
    mean_system,
    mean_wait,
    mean_response,
    throughput: self.arrival,
  }
}

///|
/// P(Wq > t), with the point mass at zero kept separate: survival(0)=C.
pub fn Mmc::wait_survival(self : Mmc, time : Double) -> Double raise QueueError {
  nonnegative(time, "waiting threshold")
  let metrics = self.analyze()
  let spare = self.service * self.servers.to_double() - self.arrival
  metrics.delay_probability * @math.exp(-spare * time)
}

///|
/// Smallest t>=0 with P(Wq<=t)>=p. p=1 has no finite quantile when C>0.
pub fn Mmc::wait_quantile(self : Mmc, p : Double) -> Double raise QueueError {
  probability(p, "quantile probability")
  let metrics = self.analyze()
  let c = metrics.delay_probability
  if c == 0.0 || p <= 1.0 - c {
    return 0.0
  }
  if p == 1.0 {
    raise Invalid("unit waiting quantile is infinite")
  }
  let spare = self.service * self.servers.to_double() - self.arrival
  checked_ratio(@math.ln(c / (1.0 - p)), spare, "waiting quantile")
}

///|
/// Conditional exponential waiting-time mean, undefined with zero arrivals.
pub fn Mmc::conditional_wait(self : Mmc) -> Double raise QueueError {
  self.validate()
  if self.arrival == 0.0 {
    raise Invalid("no delayed customers at zero arrival rate")
  }
  checked_ratio(
    1.0,
    self.service * self.servers.to_double() - self.arrival,
    "conditional waiting mean",
  )
}

///|
fn Mmc::head(self : Mmc) -> Array[Double] raise QueueError {
  self.validate()
  let a = self.arrival / self.service
  if a == 0.0 {
    let values = Array::make(self.servers + 1, 0.0)
    values[0] = 1.0
    return values
  }
  let weights : Array[Double] = [0.0]
  let mut log_weight = 0.0
  let mut maximum = 0.0
  for k in 1..<=self.servers {
    log_weight += @math.ln(a) - @math.ln(k.to_double())
    weights.push(log_weight)
    maximum = maximum.max(log_weight)
  }
  let scaled = weights.map(x => @math.exp(x - maximum))
  let mut normalizer = scaled[self.servers] /
    (1.0 - a / self.servers.to_double())
  for k in 0.. x / normalizer)
}

///|
/// Stationary P(N=n). Underflow of astronomically unlikely states is zero.
pub fn Mmc::state_probability(
  self : Mmc,
  state : Int,
) -> Double raise QueueError {
  if state < 0 || state > 1000000000 {
    raise Limit("state index must be 0..1e9")
  }
  if state <= self.servers {
    return self.head()[state]
  }
  let metrics = self.analyze()
  if metrics.utilization == 0.0 {
    return 0.0
  }
  metrics.delay_probability *
  (1.0 - metrics.utilization) *
  @math.exp((state - self.servers).to_double() * @math.ln(metrics.utilization))
}

///|
/// Tail P(N>=n); sums head probabilities without 1-CDF cancellation.
pub fn Mmc::state_tail(self : Mmc, state : Int) -> Double raise QueueError {
  if state < 0 || state > 1000000000 {
    raise Limit("tail state must be 0..1e9")
  }
  let metrics = self.analyze()
  if state == 0 {
    return 1.0
  }
  if metrics.utilization == 0.0 {
    return 0.0
  }
  if state >= self.servers {
    return metrics.delay_probability *
      @math.exp(
        (state - self.servers).to_double() * @math.ln(metrics.utilization),
      )
  }
  let head = self.head()
  let mut tail = metrics.delay_probability
  for k in state..