///|
pub(all) struct Fit {
  coefficients : Array[Double]
  fitted : Array[Double]
  residuals : Array[Double]
  sse : Double
  dof : Int32
  variance : Double?
} derive(Debug, Eq)

///|
pub extend Fit with @debug.Debug::{to_repr}

///|
pub extend Fit with Eq::{equal, not_equal}

///|
pub fn least_squares(
  a : Array[Array[Double]],
  responses : Array[Double],
  weights~ : Array[Double] = [],
  tol~ : Double = 1.0e-10,
) -> Fit raise DesignError {
  tolerance(tol)
  let (n, p) = matrix_shape(a)
  if p > 64 || n < p {
    raise Budget("least squares requires rows >= columns and <=64 terms")
  }
  if responses.length() != n || (weights.length() != 0 && weights.length() != n) {
    raise Invalid("response/weight count mismatch")
  }
  let work = zeros(n, p)
  let y = Array::make(n, 0.0)
  for i in 0.. scales[j] {
        scales[j] = work[i][j].abs()
      }
    }
    if scales[j] == 0.0 {
      raise Singular("zero model column")
    }
    for i in 0.. j ]
  let r = zeros(p, p)
  let qty = Array::make(p, 0.0)
  for k in 0..
 best {
        best = norm
        pivot = j
      }
    }
    if best <= tol * tol {
      raise Singular("model is not identifiable at requested tolerance")
    }
    if pivot != k {
      for i in 0..= 0; k = k - 1 {
    let mut value = qty[k]
    for j in (k + 1)..
 0 {
      Some(sse / dof.to_double())
    } else {
      None
    },
  }
}

///|
pub fn fit_plan(
  p : Plan,
  model : Array[Term],
  responses : Array[Double],
) -> Fit raise DesignError {
  least_squares(design_matrix(p, model), responses)
}