///|
fn zero_matrix(n : Int, p : Int) -> Array[Array[Double]] {
  Array::makei(n, _ => Array::make(p, 0.0))
}

///|
/// Multivariate Welford accumulator; rows must already be validated coordinates.
fn covariance(
  rows : Array[Array[Double]],
  ddof : Int,
) -> Array[Array[Double]] raise CodaError {
  if ddof < 0 || ddof > 1 || rows.length() <= ddof {
    raise Invalid("covariance requires n > ddof with ddof 0 or 1")
  }
  let p = rows[0].length()
  let mean = Array::make(p, 0.0)
  let m2 = zero_matrix(p, p)
  for r, row in rows {
    let delta = Array::makei(p, j => row[j] - mean[j])
    for j in 0..
 {
    Array::makei(p, j => (m2[i][j] + m2[j][i]) / (2.0 * denom))
  })
}

///|
pub fn clr_covariance(
  data : Dataset,
  ddof? : Int = 1,
) -> Array[Array[Double]] raise CodaError {
  covariance(data.clr_rows(), ddof)
}

///|
pub fn ilr_covariance(
  data : Dataset,
  b : Basis,
  ddof? : Int = 1,
) -> Array[Array[Double]] raise CodaError {
  covariance(data.ilr_rows(b), ddof)
}

///|
/// Entry (i,j) is sample variance of log(x_i/x_j).
pub fn variation_matrix(
  data : Dataset,
  ddof? : Int = 1,
) -> Array[Array[Double]] raise CodaError {
  let c = clr_covariance(data, ddof~)
  let d = c.length()
  Array::makei(d, i => {
    Array::makei(d, j => {
      let v = c[i][i] + c[j][j] - 2.0 * c[i][j]
      if v < 0.0 {
        0.0
      } else {
        v
      }
    })
  })
}

///|
pub fn total_variance(
  data : Dataset,
  ddof? : Int = 1,
) -> Double raise CodaError {
  let c = clr_covariance(data, ddof~)
  let mut total = 0.0
  for i in 0.. Array[Array[Double]] raise CodaError {
  covariance(data.cells.map(r => alr(r, reference)), ddof)
}