///|
/// Validated positive-definite matrix factorization, with triangular solves.
pub struct Cholesky {
  lower : Array[Array[Double]]
  logdet : Double
}

///|
pub fn cholesky(matrix : Array[Array[Double]]) -> Cholesky raise CodaError {
  let n = matrix.length()
  if n < 1 || n > 63 {
    raise Budget("Cholesky dimension must be 1..63")
  }
  for row in matrix {
    if row.length() != n {
      raise Shape("Cholesky input must be square")
    }
    real_vector(row)
  }
  let lower = zero_matrix(n, n)
  let mut logdet = 0.0
  for i in 0.. 1.0e-12 {
        raise Invalid("Cholesky input is not symmetric")
      }
      let mut value = matrix[i][j]
      for k in 0.. Array[Array[Double]] {
  self.lower.map(r => r.copy())
}

///|
pub fn Cholesky::log_determinant(self : Cholesky) -> Double {
  self.logdet
}

///|
fn Cholesky::whiten(
  self : Cholesky,
  rhs : Array[Double],
) -> Array[Double] raise CodaError {
  if rhs.length() != self.lower.length() {
    raise Shape("triangular solve dimension mismatch")
  }
  real_vector(rhs)
  let out = rhs.copy()
  for i in 0.. Array[Double] raise CodaError {
  let out = self.whiten(rhs)
  for offset in 0.. Double raise CodaError {
  let y = self.whiten(rhs)
  let mut result = 0.0
  for value in y {
    result = checked(result + value * value)
  }
  result
}