///|
/// Validated one-dimensional piecewise-linear table.
pub(all) struct LinearTable {
  xs : Array[Double]
  ys : Array[Double]
  extrapolation : Extrapolation
} derive(Debug, Eq)

///|
pub(all) enum Extrapolation {
  Reject
  Clamp
  Linear
} derive(Debug, Eq)

///|
pub fn LinearTable::new(
  xs : ArrayView[Double],
  ys : ArrayView[Double],
  extrapolation? : Extrapolation = Reject,
) -> LinearTable raise @domain.ExtractionError {
  if xs.length() != ys.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="x and y columns must have the same length",
    )
  }
  if xs.length() < 2 {
    raise @domain.ExtractionError::InsufficientData(
      required=2,
      actual=xs.length(),
    )
  }
  for index, x in xs {
    @domain.require_finite("x[\{index}]", x)
    @domain.require_finite("y[\{index}]", ys[index])
    if index > 0 {
      if x < xs[index - 1] {
        raise @domain.ExtractionError::UnsortedData(index~)
      }
      if x == xs[index - 1] {
        raise @domain.ExtractionError::DuplicateCoordinate(value=x)
      }
    }
  }
  { xs: xs.to_owned(), ys: ys.to_owned(), extrapolation }
}

///|
pub fn LinearTable::lower(self : LinearTable) -> Double {
  self.xs[0]
}

///|
pub fn LinearTable::upper(self : LinearTable) -> Double {
  self.xs[self.xs.length() - 1]
}

///|
fn LinearTable::segment(self : LinearTable, x : Double) -> Int {
  let last = self.xs.length() - 1
  if x <= self.xs[0] {
    0
  } else if x >= self.xs[last] {
    last - 1
  } else {
    for low = 0, high = last; low + 1 < high; {
      let middle = low + (high - low) / 2
      if self.xs[middle] <= x {
        continue middle, high
      } else {
        continue low, middle
      }
    } nobreak {
      low
    }
  }
}

///|
fn interpolate_segment(
  x : Double,
  x0 : Double,
  x1 : Double,
  y0 : Double,
  y1 : Double,
) -> Double {
  let fraction = (x - x0) / (x1 - x0)
  y0 + fraction * (y1 - y0)
}

///|
pub fn LinearTable::evaluate(
  self : LinearTable,
  x : Double,
) -> Double raise @domain.ExtractionError {
  @domain.require_finite("query", x)
  let lower = self.lower()
  let upper = self.upper()
  match self.extrapolation {
    Reject =>
      if x < lower || x > upper {
        raise @domain.ExtractionError::OutOfRange(
          name="query",
          value=x,
          lower~,
          upper~,
        )
      }
    Clamp =>
      if x <= lower {
        return self.ys[0]
      } else if x >= upper {
        return self.ys[self.ys.length() - 1]
      }
    Linear => ()
  }
  let index = self.segment(x)
  interpolate_segment(
    x,
    self.xs[index],
    self.xs[index + 1],
    self.ys[index],
    self.ys[index + 1],
  )
}

///|
pub fn LinearTable::slope_at(
  self : LinearTable,
  x : Double,
) -> Double raise @domain.ExtractionError {
  ignore(self.evaluate(x))
  let index = self.segment(x)
  (self.ys[index + 1] - self.ys[index]) / (self.xs[index + 1] - self.xs[index])
}

///|
/// Integrate the piecewise-linear curve using exact trapezoids per segment.
pub fn LinearTable::integrate(
  self : LinearTable,
  from~ : Double,
  to~ : Double,
) -> Double raise @domain.ExtractionError {
  if from == to {
    return 0.0
  }
  if from > to {
    return -self.integrate(from=to, to=from)
  }
  if from < self.lower() || to > self.upper() {
    raise @domain.ExtractionError::OutOfRange(
      name="integration interval",
      value=if from < self.lower() { from } else { to },
      lower=self.lower(),
      upper=self.upper(),
    )
  }
  let mut cursor = from
  let mut area = 0.0
  for ; cursor < to; {
    let index = self.segment(cursor)
    let right = self.xs[index + 1].min(to)
    let y0 = self.evaluate(cursor)
    let y1 = self.evaluate(right)
    area += 0.5 * (y0 + y1) * (right - cursor)
    cursor = right
  }
  area
}