///|
pub fn bracket_index(values : ArrayView[Double], query : Double) -> Int? {
  if values.length() < 2 ||
    query < values[0] ||
    query > values[values.length() - 1] {
    None
  } else {
    for i in 1.. Double {
  if right_x == left_x {
    (left_y + right_y) / 2.0
  } else {
    let fraction = ((query - left_x) / (right_x - left_x)).clamp(
      min=0.0,
      max=1.0,
    )
    lerp(left_y, right_y, fraction)
  }
}

///|
pub fn linear_interpolate(
  coordinates : ArrayView[Double],
  values : ArrayView[Double],
  query : Double,
) -> Double? {
  let count = Int::min(coordinates.length(), values.length())
  if count == 0 {
    None
  } else if count == 1 {
    Some(values[0])
  } else {
    match bracket_index(coordinates[:count], query) {
      Some(index) =>
        Some(
          linear_interpolate_clamped(
            coordinates[index],
            values[index],
            coordinates[index + 1],
            values[index + 1],
            query,
          ),
        )
      None =>
        if query < coordinates[0] {
          Some(values[0])
        } else {
          Some(values[count - 1])
        }
    }
  }
}

///|
pub fn nearest_index(values : ArrayView[Double], query : Double) -> Int? {
  if values.length() == 0 {
    None
  } else {
    let mut best = 0
    let mut distance = (values[0] - query).abs()
    for i in 1.. Double? {
  let count = Int::min(coordinates.length(), values.length())
  match nearest_index(coordinates[:count], query) {
    Some(index) => Some(values[index])
    None => None
  }
}

///|
pub fn resample_array_linear(
  coordinates : ArrayView[Double],
  values : ArrayView[Double],
  queries : ArrayView[Double],
) -> Array[Double] {
  queries.map(fn(query) {
    linear_interpolate(coordinates, values, query).unwrap_or(0.0 / 0.0)
  })
}

///|
pub fn resample_array_nearest(
  coordinates : ArrayView[Double],
  values : ArrayView[Double],
  queries : ArrayView[Double],
) -> Array[Double] {
  queries.map(fn(query) {
    nearest_interpolate(coordinates, values, query).unwrap_or(0.0 / 0.0)
  })
}

///|
pub fn finite_linear_interpolate(
  coordinates : ArrayView[Double],
  values : ArrayView[Double],
  query : Double,
) -> Double? {
  let filtered_x : Array[Double] = []
  let filtered_y : Array[Double] = []
  for i in 0.. Array[Double] {
  let output = values.to_owned()
  for i in 0.. output[i] = value
        None => ()
      }
    }
  }
  output
}

///|
pub fn moving_linear_smooth(
  values : ArrayView[Double],
  radius : Int,
) -> Array[Double] {
  let safe = Int::max(0, radius)
  let output : Array[Double] = []
  for i in 0.. Array[Double] {
  let safe = Int::max(0, radius)
  let safe_decay = Double::max(0.000001, decay)
  let output : Array[Double] = []
  for i in 0.. Array[Double] {
  let safe = Int::max(0, radius)
  let safe_sigma = Double::max(0.000001, sigma)
  let output : Array[Double] = []
  let mut total = 0.0
  for i in -safe..<=safe {
    let distance = Double::from_int(i)
    let value = @math.exp(
      -distance * distance / (2.0 * safe_sigma * safe_sigma),
    )
    output.push(value)
    total += value
  }
  if total == 0.0 {
    output
  } else {
    output.map(fn(value) { value / total })
  }
}

///|
pub fn convolve_1d(
  values : ArrayView[Double],
  kernel : ArrayView[Double],
) -> Array[Double] {
  let output : Array[Double] = []
  let centre = kernel.length() / 2
  for i in 0..= 0 && index < values.length() && is_finite(values[index]) {
        total += values[index] * kernel[k]
        weight_total += kernel[k]
      }
    }
    output.push(if weight_total == 0.0 { total } else { total / weight_total })
  }
  output
}

///|
pub fn gaussian_smooth(
  values : ArrayView[Double],
  radius : Int,
  sigma : Double,
) -> Array[Double] {
  convolve_1d(values, gaussian_kernel(radius, sigma))
}

///|
pub fn cubic_hermite(
  y0 : Double,
  y1 : Double,
  y2 : Double,
  y3 : Double,
  fraction : Double,
) -> Double {
  let t2 = fraction * fraction
  let t3 = t2 * fraction
  0.5 *
  (
    2.0 * y1 +
    (-y0 + y2) * fraction +
    (2.0 * y0 - 5.0 * y1 + 4.0 * y2 - y3) * t2 +
    (-y0 + 3.0 * y1 - 3.0 * y2 + y3) * t3
  )
}

///|
pub fn cubic_interpolate(
  coordinates : ArrayView[Double],
  values : ArrayView[Double],
  query : Double,
) -> Double? {
  let count = Int::min(coordinates.length(), values.length())
  if count < 2 {
    None
  } else {
    let index = bracket_index(coordinates[:count], query).unwrap_or(0)
    let left = Int::max(0, index - 1)
    let right = Int::min(count - 1, index + 2)
    if right - left < 3 {
      linear_interpolate(coordinates, values, query)
    } else {
      let fraction = (query - coordinates[index]) /
        (coordinates[index + 1] - coordinates[index])
      Some(
        cubic_hermite(
          values[left],
          values[index],
          values[index + 1],
          values[right],
          fraction,
        ),
      )
    }
  }
}

///|
pub fn monotone_envelope(
  values : ArrayView[Double],
  increasing : Bool,
) -> Array[Double] {
  let output = values.to_owned()
  if output.length() > 0 {
    for i in 1.. Bool {
  for i in 1.. values[i - 1] {
      return false
    }
  }
  true
}

///|
pub fn interpolate_series_cubic(
  series : FieldSeries,
  coordinate : Double,
) -> Double? {
  let ordered = if series_is_sorted(series) {
    series
  } else {
    series_sort(series)
  }
  let coordinates = ordered.points.map(fn(point) { point.coordinate })
  let values = ordered.points.map(fn(point) { point.value })
  cubic_interpolate(coordinates, values, coordinate)
}

///|
pub fn interpolate_grid_value(
  grid : Grid,
  point : GeoPoint,
  nearest : Bool,
) -> Double? {
  if nearest {
    grid.sample_nearest(point)
  } else {
    grid.sample_bilinear(point)
  }
}

///|
pub fn resample_grid_like(source : Grid, target : Grid, nearest : Bool) -> Grid {
  target.map_cells(fn(cell) {
    interpolate_grid_value(source, cell.point, nearest).unwrap_or(0.0 / 0.0)
  })
}

///|
pub fn interpolate_grid_profile(
  grid : Grid,
  points : ArrayView[GeoPoint],
  nearest : Bool,
) -> Array[PointEstimate] {
  let output : Array[PointEstimate] = []
  for point in points {
    output.push({
      point,
      value: interpolate_grid_value(grid, point, nearest).unwrap_or(0.0 / 0.0),
      neighbours: if nearest {
        1
      } else {
        4
      },
    })
  }
  output
}

///|
pub fn inverse_distance_weights(
  distances : ArrayView[Double],
  power : Double,
) -> Array[Double] {
  let safe_power = Double::max(0.1, power)
  let raw : Array[Double] = []
  let mut total = 0.0
  for distance in distances {
    let weight = if distance == 0.0 {
      1.0e12
    } else {
      1.0 / @math.pow(distance, safe_power)
    }
    raw.push(weight)
    total += weight
  }
  if total == 0.0 {
    raw.map(fn(_) { 0.0 })
  } else {
    raw.map(fn(value) { value / total })
  }
}

///|
pub fn weighted_interpolate(
  values : ArrayView[Double],
  weights : ArrayView[Double],
) -> Double {
  let mut total = 0.0
  let mut weight_total = 0.0
  for i in 0..