///|
pub(all) enum InterpolationMode {
  Nearest
  Linear
  Quadratic
} derive(Eq, Debug, ToJson)

///|
pub(all) struct InterpolationResult {
  value : Double
  left : Int
  right : Int
  fraction : Double
} derive(Debug, ToJson)

///|
pub fn locate_periodic(grid : Grid1D, x : Double) -> InterpolationResult {
  let y = grid.wrap(x)
  let scaled = y / grid.dx - 0.5
  let left0 = @math.floor(scaled).to_int()
  let fraction = scaled - left0.to_double()
  let left = modulo_index(left0, grid.cells)
  let right = modulo_index(left + 1, grid.cells)
  { value: scaled, left, right, fraction }
}

///|
pub fn interpolate_pair(
  left : Double,
  right : Double,
  fraction : Double,
  mode : InterpolationMode,
) -> Double {
  match mode {
    Nearest => if fraction < 0.5 { left } else { right }
    Linear => lerp(left, right, fraction)
    Quadratic => {
      let smooth = fraction * fraction * (3.0 - 2.0 * fraction)
      lerp(left, right, smooth)
    }
  }
}

///|
pub fn sample_periodic(
  grid : Grid1D,
  values : ArrayView[Double],
  x : Double,
  mode : InterpolationMode,
) -> Double {
  if values.length() == 0 {
    0.0
  } else {
    let location = locate_periodic(grid, x)
    interpolate_pair(
      values[location.left],
      values[location.right],
      location.fraction,
      mode,
    )
  }
}

///|
pub fn sample_periodic_linear(
  grid : Grid1D,
  values : ArrayView[Double],
  x : Double,
) -> Double {
  sample_periodic(grid, values, x, Linear)
}

///|
pub fn nearest_index(grid : Grid1D, x : Double) -> Int {
  let location = locate_periodic(grid, x)
  if location.fraction < 0.5 {
    location.left
  } else {
    location.right
  }
}

///|
pub fn resample_periodic(
  source_grid : Grid1D,
  source : ArrayView[Double],
  target_grid : Grid1D,
  mode : InterpolationMode,
) -> Array[Double] {
  let output = zeros(target_grid.cells)
  for i in 0.. Double {
  if values.length() == 0 {
    0.0
  } else {
    sum_values(values) * grid.dx
  }
}

///|
pub fn periodic_average(
  grid : Grid1D,
  values : ArrayView[Double],
  radius : Int,
) -> Array[Double] {
  let output = zeros(grid.cells)
  if values.length() == grid.cells {
    for i in 0.. Array[Double] {
  let output = zeros(grid.cells)
  let mut total = 0.0
  for i in 0.. Double {
  if values.length() != grid.cells {
    0.0
  } else {
    let left = values[modulo_index(index - 1, grid.cells)]
    let right = values[modulo_index(index + 1, grid.cells)]
    (right - left) / (2.0 * grid.dx)
  }
}

///|
pub fn forward_difference(
  grid : Grid1D,
  values : ArrayView[Double],
  index : Int,
) -> Double {
  if values.length() != grid.cells {
    0.0
  } else {
    (
      values[modulo_index(index + 1, grid.cells)] -
      values[modulo_index(index, grid.cells)]
    ) /
    grid.dx
  }
}

///|
pub fn backward_difference(
  grid : Grid1D,
  values : ArrayView[Double],
  index : Int,
) -> Double {
  if values.length() != grid.cells {
    0.0
  } else {
    (
      values[modulo_index(index, grid.cells)] -
      values[modulo_index(index - 1, grid.cells)]
    ) /
    grid.dx
  }
}

///|
pub fn interpolation_error(
  grid : Grid1D,
  values : ArrayView[Double],
  points : ArrayView[Double],
  expected : ArrayView[Double],
) -> Double {
  if points.length() != expected.length() || points.length() == 0 {
    0.0
  } else {
    let mut total = 0.0
    for i in 0..