///|
/// Numerical health signals used by reports and quality gates.
pub(all) struct NumericalHealth {
  finite : Bool
  positive_density : Bool
  mass_residual : Double
  max_mach : Double
  min_density : Double
  max_density : Double
  pass : Bool
} derive(Debug)

///|
/// Check a distribution vector for finite and nonnegative populations.
pub fn distribution_is_physical(distribution : ArrayView[Double]) -> Bool {
  let mut physical = true
  for value in distribution {
    physical = physical && !value.is_nan() && !value.is_inf() && value >= 0.0
  }
  physical
}

///|
/// Maximum absolute deviation from equilibrium moments.
pub fn equilibrium_error(
  distribution : ArrayView[Double],
  rho~ : Double,
  ux~ : Double,
  uy~ : Double,
) -> Double {
  let expected = equilibrium_distribution(rho~, ux~, uy~)
  linf_error(expected, distribution)
}

///|
/// Check density/velocity fields for finite values.
pub fn fields_are_finite(density : Field2D, velocity : VectorField2D) -> Bool {
  let mut result = true
  for value in density.data {
    result = result && !value.is_nan() && !value.is_inf()
  }
  for value in velocity.ux {
    result = result && !value.is_nan() && !value.is_inf()
  }
  for value in velocity.uy {
    result = result && !value.is_nan() && !value.is_inf()
  }
  result
}

///|
/// Compute health metrics for a reusable simulation.
pub fn Simulation::health(self : Simulation) -> NumericalHealth {
  let stability = self.stability()
  let velocity = self.velocity_field()
  let finite = !stability.has_nan &&
    fields_are_finite(self.density_field(), velocity)
  let positive_density = stability.min_rho > self.options.density_floor
  let max_mach = mach_number(ux=max_speed(velocity), uy=0.0)
  let expected_mass = self.size.width.to_double() * self.size.height.to_double()
  let mass_residual = abs_double(stability.mass - expected_mass)
  {
    finite,
    positive_density,
    mass_residual,
    max_mach,
    min_density: stability.min_rho,
    max_density: stability.max_rho,
    pass: finite && positive_density && max_mach < self.options.max_mach,
  }
}

///|
/// Compute a finite residual between two scalar fields.
pub fn field_residual(left : Field2D, right : Field2D) -> Double {
  relative_l2_error(left.data, right.data)
}

///|
/// Check that a field has no negative values below a tolerance.
pub fn field_is_nonnegative(
  field : Field2D,
  tolerance? : Double = 0.000000001,
) -> Bool {
  let mut result = true
  for value in field.data {
    result = result && value >= -tolerance
  }
  result
}

///|
/// Return the fraction of finite values in a scalar field.
pub fn finite_fraction(field : Field2D) -> Double {
  let mut count = 0
  for value in field.data {
    if !value.is_nan() && !value.is_inf() {
      count += 1
    }
  }
  safe_divide(count.to_double(), field.data.length().to_double(), fallback=1.0)
}

///|
/// Return the maximum absolute speed in a simulation.
pub fn Simulation::max_mach(self : Simulation) -> Double {
  mach_number(ux=max_speed(self.velocity_field()), uy=0.0)
}

///|
/// Return the normalized density variation.
pub fn Simulation::density_variation(self : Simulation) -> Double {
  let statistics = self.density_field().statistics()
  safe_divide(
    statistics.maximum - statistics.minimum,
    statistics.mean,
    fallback=0.0,
  )
}