///|
/// Summary statistics for a scalar field.
pub(all) struct FieldStatistics {
  count : Int
  sum : Double
  mean : Double
  minimum : Double
  maximum : Double
  variance : Double
  l2_norm : Double
} derive(Debug)

///|
/// A row-major scalar field with clamped interpolation support.
pub(all) struct Field2D {
  size : Size
  data : Array[Double]
} derive(Debug)

///|
/// Allocate a scalar field.
pub fn Field2D::new(size~ : Size, initial? : Double = 0.0) -> Field2D {
  { size, data: Array::make(size.width.max(0) * size.height.max(0), initial) }
}

///|
/// Return whether integer coordinates are inside the field.
pub fn Field2D::contains(self : Field2D, x : Int, y : Int) -> Bool {
  x >= 0 && y >= 0 && x < self.size.width && y < self.size.height
}

///|
/// Convert a coordinate to a row-major index.
pub fn Field2D::index(self : Field2D, x : Int, y : Int) -> Int {
  y * self.size.width + x
}

///|
/// Read a value, returning zero outside the domain.
pub fn Field2D::get(self : Field2D, x : Int, y : Int) -> Double {
  if self.contains(x, y) {
    self.data[self.index(x, y)]
  } else {
    0.0
  }
}

///|
/// Write a value when the coordinate is inside the domain.
pub fn Field2D::set(
  self : Field2D,
  x~ : Int,
  y~ : Int,
  value~ : Double,
) -> Unit {
  if self.contains(x, y) {
    self.data[self.index(x, y)] = value
  }
}

///|
/// Fill all cells with one value.
pub fn Field2D::fill(self : Field2D, value~ : Double) -> Unit {
  self.data.fill(value)
}

///|
/// Copy the underlying data in row-major order.
pub fn Field2D::values(self : Field2D) -> Array[Double] {
  self.data.copy()
}

///|
/// Apply a pointwise transformation and return a new field.
pub fn Field2D::map(
  self : Field2D,
  transform : (Int, Int, Double) -> Double,
) -> Field2D {
  let result = Field2D::new(size=self.size)
  for y in 0.. Double {
  if self.size.width == 0 || self.size.height == 0 {
    0.0
  } else {
    let xx = clamp_double(x, low=0.0, high=(self.size.width - 1).to_double())
    let yy = clamp_double(y, low=0.0, high=(self.size.height - 1).to_double())
    let x0 = xx.floor().to_int()
    let y0 = yy.floor().to_int()
    let x1 = (x0 + 1).min(self.size.width - 1)
    let y1 = (y0 + 1).min(self.size.height - 1)
    let tx = xx - x0.to_double()
    let ty = yy - y0.to_double()
    let a = lerp(self.get(x0, y0), self.get(x1, y0), tx)
    let b = lerp(self.get(x0, y1), self.get(x1, y1), tx)
    lerp(a, b, ty)
  }
}

///|
/// Compute unmasked descriptive statistics.
pub fn Field2D::statistics(self : Field2D) -> FieldStatistics {
  let accumulator = StatisticsAccumulator::new()
  for value in self.data {
    accumulator.push(value)
  }
  accumulator.finish()
}

///|
/// Compute statistics while optionally excluding solid mask cells.
pub fn Field2D::statistics_masked(
  self : Field2D,
  mask~ : DomainMask,
  include_solid? : Bool = false,
) -> FieldStatistics {
  let accumulator = StatisticsAccumulator::new()
  for y in 0.. VectorField2D {
  let cells = size.width.max(0) * size.height.max(0)
  { size, ux: Array::make(cells, 0.0), uy: Array::make(cells, 0.0) }
}

///|
/// Set a velocity at one cell.
pub fn VectorField2D::set(
  self : VectorField2D,
  x~ : Int,
  y~ : Int,
  ux~ : Double,
  uy~ : Double,
) -> Unit {
  if x >= 0 && y >= 0 && x < self.size.width && y < self.size.height {
    let index = y * self.size.width + x
    self.ux[index] = ux
    self.uy[index] = uy
  }
}

///|
/// Fill a velocity field with a uniform vector.
pub fn VectorField2D::fill(
  self : VectorField2D,
  ux~ : Double,
  uy~ : Double,
) -> Unit {
  self.ux.fill(ux)
  self.uy.fill(uy)
}

///|
/// Read a velocity with zero extension outside the grid.
pub fn VectorField2D::get(self : VectorField2D, x : Int, y : Int) -> Point {
  if x >= 0 && y >= 0 && x < self.size.width && y < self.size.height {
    let index = y * self.size.width + x
    Point::new(x=self.ux[index], y=self.uy[index])
  } else {
    Point::new(x=0.0, y=0.0)
  }
}

///|
/// Sample both velocity components bilinearly.
pub fn VectorField2D::sample(
  self : VectorField2D,
  x~ : Double,
  y~ : Double,
) -> Point {
  let x_field = Field2D::{ size: self.size, data: self.ux }
  let y_field = Field2D::{ size: self.size, data: self.uy }
  Point::new(x=x_field.sample(x~, y~), y=y_field.sample(x~, y~))
}