///|
/// A one-dimensional relaxation parameter sweep.
pub(all) struct ParameterSweep {
  omegas : Array[Double]
  steps : Int
  boundary : BoundaryMode
} derive(Debug)

///|
/// Build a sweep over supplied relaxation values.
pub fn ParameterSweep::new(
  omegas~ : Array[Double],
  steps~ : Int,
  boundary? : BoundaryMode = BounceBack,
) -> ParameterSweep {
  { omegas: omegas.copy(), steps: steps.max(0), boundary }
}

///|
/// One parameter sweep result.
pub(all) struct SweepResult {
  omega : Double
  steps : Int
  mass : Double
  max_speed : Double
  stable : Bool
  health_pass : Bool
} derive(Debug)

///|
/// Execute every point in a relaxation sweep.
pub fn run_parameter_sweep(
  size~ : Size,
  sweep~ : ParameterSweep,
) -> Array[SweepResult] {
  let result = Array::new()
  for omega in sweep.omegas {
    let simulation = Simulation::new(
      size~,
      options=SimulationOptions::new(
        model=CollisionModel::mrt(omega~),
        boundary=sweep.boundary,
      ),
    )
    simulation.run(sweep.steps)
    let stability = simulation.stability()
    let health = simulation.health()
    result.push({
      omega,
      steps: sweep.steps,
      mass: simulation.mass(),
      max_speed: stability.max_speed,
      stable: stability.recommended,
      health_pass: health.pass,
    })
  }
  result
}

///|
/// Return the stable sweep result with the smallest maximum speed.
pub fn best_sweep_result(results : ArrayView[SweepResult]) -> SweepResult? {
  let mut best : SweepResult? = None
  for result in results {
    if result.stable {
      match best {
        None => best = Some(result)
        Some(current) =>
          if result.max_speed < current.max_speed {
            best = Some(result)
          }
      }
    }
  }
  best
}

///|
/// Serialize sweep results as CSV.
pub fn sweep_results_csv(results : ArrayView[SweepResult]) -> String {
  let builder = StringBuilder(size_hint=96 + results.length() * 80)
  builder.write_string("omega,steps,mass,max_speed,stable,health_pass\n")
  for result in results {
    builder.write_string(
      "\{result.omega},\{result.steps},\{result.mass},\{result.max_speed},\{result.stable},\{result.health_pass}\n",
    )
  }
  builder.to_string()
}

///|
/// Create a regular relaxation grid, inclusive of its endpoints.
pub fn omega_grid(
  start~ : Double,
  end~ : Double,
  count~ : Int,
) -> Array[Double] {
  let result = Array::new()
  let n = count.max(0)
  if n == 1 {
    result.push(start)
  } else if n > 1 {
    for i in 0.. Double {
  let mut stable = 0
  for result in results {
    if result.stable {
      stable += 1
    }
  }
  safe_divide(stable.to_double(), results.length().to_double(), fallback=0.0)
}