///|
/// M/G/1 FCFS means. The second moment is E[S^2], not variance or stddev.
pub fn mg1(
  arrival : Double,
  mean_service : Double,
  second_moment : Double,
) -> QueueMetrics raise QueueError {
  rate(arrival, "arrival")
  positive(mean_service, "mean service")
  nonnegative(second_moment, "service second moment")
  let square = mean_service * mean_service
  finite(square, "squared service mean")
  if second_moment < square {
    raise Invalid("E[S^2] must be at least E[S]^2")
  }
  let utilization = arrival * mean_service
  finite(utilization, "utilization")
  if utilization >= 1.0 {
    raise Unstable("M/G/1 utilization must be below one")
  }
  let mean_wait = checked_ratio(
    arrival * second_moment,
    2.0 * (1.0 - utilization),
    "Pollaczek-Khinchine waiting mean",
  )
  let mean_response = mean_wait + mean_service
  let mean_queue = arrival * mean_wait
  let mean_system = arrival * mean_response
  for x in [mean_response, mean_queue, mean_system] {
    finite(x, "M/G/1 metric")
  }
  {
    utilization,
    delay_probability: utilization,
    mean_queue,
    mean_system,
    mean_wait,
    mean_response,
    throughput: arrival,
  }
}

///|
/// An explicitly approximate GI/G/c mean model using the variability factor
/// (Ca^2+Cs^2)/2 times the exact M/M/c wait. Not a distributional model.
pub(all) struct ApproximateQueue {
  mean_wait : Double
  mean_response : Double
  mean_queue : Double
  utilization : Double
  variability_factor : Double
} derive(Debug, ToJson)

///|
pub fn gigc_approximation(
  model : Mmc,
  arrival_scv : Double,
  service_scv : Double,
) -> ApproximateQueue raise QueueError {
  nonnegative(arrival_scv, "arrival SCV")
  nonnegative(service_scv, "service SCV")
  if arrival_scv > 1000000.0 || service_scv > 1000000.0 {
    raise Limit("SCV must not exceed 1e6")
  }
  let baseline = model.analyze()
  let variability_factor = (arrival_scv + service_scv) / 2.0
  let mean_wait = variability_factor * baseline.mean_wait
  let mean_response = mean_wait + 1.0 / model.service
  let mean_queue = model.arrival * mean_wait
  for x in [mean_wait, mean_response, mean_queue] {
    finite(x, "GI/G/c approximation")
  }
  {
    mean_wait,
    mean_response,
    mean_queue,
    utilization: baseline.utilization,
    variability_factor,
  }
}