///|
/// Local polynomial smoothing and derivatives on physical coordinates.
/// Uniform grids reproduce Savitzky-Golay mode=interp; irregular grids use
/// an actual local fit. At edges a complete shifted window is fitted.
pub fn local_polynomial(
  input : Spectrum,
  window~ : Int32 = 7,
  degree~ : Int32 = 2,
  derivative~ : Int32 = 0,
) -> Spectrum raise SpectralError {
  if window < 3 || window > 101 || window % 2 == 0 || window > input.length() {
    raise Invalid("local window must be odd, 3..101, and <= sample count")
  }
  if degree < 0 || degree > 8 || degree >= window {
    raise Invalid("local degree must be 0..8 and less than window")
  }
  if derivative < 0 || derivative > degree {
    raise Invalid("derivative exceeds fitted degree")
  }
  let n = input.length()
  let half = window / 2
  let output = Array::make(n, 0.0)
  for i in 0.. n {
      start = n - window
    }
    let x = Array::makei(window, fn(k) { input.axis[start + k] })
    let y = Array::makei(window, fn(k) { input.values[start + k] })
    let fit = fit_polynomial(x, y, degree)
    output[i] = fit.evaluate(input.axis[i], derivative~)
  }
  spectrum(input.axis, output)
}

///|
/// Median filter with a shifted complete window at edges.
pub fn median_filter(
  input : Spectrum,
  window~ : Int32 = 5,
) -> Spectrum raise SpectralError {
  if window < 3 || window > 101 || window % 2 == 0 || window > input.length() {
    raise Invalid("median window must be odd, 3..101, and <= sample count")
  }
  let output = Array::make(input.length(), 0.0)
  for i in 0.. input.length() {
      start = input.length() - window
    }
    let values = Array::makei(window, fn(j) { input.values[start + j] })
    output[i] = quantile(values, 0.5)
  }
  spectrum(input.axis, output)
}