///|
/// Small, dependency-free regression primitives for measured equilibrium data.
/// The implementation deliberately keeps the normal equations visible so that
/// a process engineer can audit the fit and its residuals without a black box.
pub(all) struct Observation {
  x : Double
  y : Double
  weight : Double
} derive(Debug, Eq, ToJson)

///|
pub fn Observation::new(
  x~ : Double,
  y~ : Double,
  weight? : Double = 1.0,
) -> Observation raise @domain.ExtractionError {
  @domain.require_finite("observation x", x)
  @domain.require_finite("observation y", y)
  @domain.require_positive("observation weight", weight)
  { x, y, weight }
}

///|
pub(all) struct RegressionModel {
  coefficients : Array[Double]
  degree : Int
  observations : Int
  weighted_mean_y : Double
  residual_sum_squares : Double
  total_sum_squares : Double
  r_squared : Double
  root_mean_squared_error : Double
  max_absolute_error : Double
} derive(Debug, Eq, ToJson)

///|
pub fn RegressionModel::predict(self : RegressionModel, x : Double) -> Double {
  let mut result = 0.0
  for coefficient in self.coefficients.rev() {
    result = result * x + coefficient
  }
  result
}

///|
pub fn RegressionModel::residual(
  self : RegressionModel,
  observation : Observation,
) -> Double {
  observation.y - self.predict(observation.x)
}

///|
pub fn RegressionModel::coefficient(
  self : RegressionModel,
  index : Int,
) -> Double raise @domain.ExtractionError {
  match self.coefficients.get(index) {
    Some(value) => value
    None =>
      raise @domain.ExtractionError::InvalidArgument(
        message="coefficient index out of range",
      )
  }
}

///|
pub fn RegressionModel::summary(self : RegressionModel) -> String {
  "degree=\{self.degree},n=\{self.observations},r2=\{self.r_squared},rmse=\{self.root_mean_squared_error},max_abs_error=\{self.max_absolute_error}"
}

///|
fn solve_linear_system(
  matrix : Array[Array[Double]],
  rhs : Array[Double],
) -> Array[Double] raise @domain.ExtractionError {
  let n = rhs.length()
  let a = matrix
  let b = rhs
  for pivot in 0.. magnitude {
        best = row
        magnitude = a[row][pivot].abs()
      }
    }
    if magnitude <= 1.0e-14 {
      raise @domain.ExtractionError::InvalidArgument(
        message="regression design matrix is singular",
      )
    }
    if best != pivot {
      let row = a[pivot]
      a[pivot] = a[best]
      a[best] = row
      let value = b[pivot]
      b[pivot] = b[best]
      b[best] = value
    }
    for row in (pivot + 1).. 0.0 ]
  for offset in 0.. RegressionModel raise @domain.ExtractionError {
  if observations.length() < degree + 1 {
    raise @domain.ExtractionError::InsufficientData(
      required=degree + 1,
      actual=observations.length(),
    )
  }
  if degree < 0 || degree > 8 {
    raise @domain.ExtractionError::InvalidArgument(
      message="polynomial degree must be between 0 and 8",
    )
  }
  let size = degree + 1
  let normal : Array[Array[Double]] = [
    for _ in 0.. [ for _ in 0.. 0.0 ]
  ]
  let rhs : Array[Double] = [ for _ in 0.. 0.0 ]
  let mut weight_total = 0.0
  let mut weighted_y = 0.0
  for item in observations {
    let powers : Array[Double] = [ for _ in 0..<(2 * degree + 1) => 1.0 ]
    for index in 1.. RegressionModel raise @domain.ExtractionError {
  fit_polynomial(observations, 1)
}

///|
pub(all) struct CrossValidation {
  folds : Int
  mean_squared_error : Double
  worst_squared_error : Double
  predictions : Array[Double]
} derive(Debug, Eq, ToJson)

///|
/// Leave-one-out validation is useful for small laboratory equilibrium tables.
pub fn leave_one_out(
  observations : ArrayView[Observation],
  degree : Int,
) -> CrossValidation raise @domain.ExtractionError {
  if observations.length() < degree + 2 {
    raise @domain.ExtractionError::InsufficientData(
      required=degree + 2,
      actual=observations.length(),
    )
  }
  let predictions : Array[Double] = [
    for _ in 0.. 0.0
  ]
  let mut sum = 0.0
  let mut worst = 0.0
  for held_out in 0.. EfficiencyCalibration raise @domain.ExtractionError {
  if ideal_stages.length() != actual_stages.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="ideal and actual stage arrays must have equal length",
    )
  }
  if ideal_stages.length() < 2 {
    raise @domain.ExtractionError::InsufficientData(
      required=2,
      actual=ideal_stages.length(),
    )
  }
  let observations = []
  for index, ideal in ideal_stages {
    @domain.require_positive("ideal stage count", ideal)
    @domain.require_positive("actual stage count", actual_stages[index])
    observations.push(Observation::new(x=ideal, y=ideal / actual_stages[index]))
  }
  let model = fit_linear(observations)
  let predictions = []
  let mut clamped = 0
  for ideal in ideal_stages {
    let raw = model.predict(ideal)
    if raw <= 0.0 {
      clamped += 1
      predictions.push(1.0e-12)
    } else if raw > 1.0 {
      clamped += 1
      predictions.push(1.0)
    } else {
      predictions.push(raw)
    }
  }
  { model, predicted_efficiencies: predictions, clamped_count: clamped }
}