///|
/// Rows are orthonormal, zero-sum contrasts; columns follow feature order.
pub struct Basis {
  contrasts : Array[Array[Double]]
  parts : Int
}

///|
pub fn basis(rows : Array[Array[Double]]) -> Basis raise CodaError {
  let d = rows.length() + 1
  vector_size(d)
  for i, row in rows {
    if row.length() != d {
      raise Shape("basis must be (D-1) by D")
    }
    real_vector(row)
    let mut total = 0.0
    for v in row {
      total = checked(total + v)
    }
    if total.abs() > 1.0e-10 {
      raise Invalid("basis contrast is not zero-sum")
    }
    for j in 0..<(i + 1) {
      let mut dot = 0.0
      for k in 0.. 1.0e-10 {
        raise Invalid("basis is not orthonormal")
      }
    }
  }
  { contrasts: rows.map(r => r.copy()), parts: d, }
}

///|
pub fn helmert_basis(parts : Int) -> Basis raise CodaError {
  vector_size(parts)
  let rows = []
  for k in 1.. {
        if i < k {
          1.0 / denom
        } else if i == k {
          -k.to_double() / denom
        } else {
          0.0
        }
      }),
    )
  }
  basis(rows)
}

///|
pub fn Basis::rows(self : Basis) -> Array[Array[Double]] {
  self.contrasts.map(r => r.copy())
}

///|
pub fn Basis::dimension(self : Basis) -> Int {
  self.parts
}

///|
pub fn ilr(x : Array[Double], b : Basis) -> Array[Double] raise CodaError {
  if x.length() != b.parts {
    raise Shape("ILR basis dimension mismatch")
  }
  let z = clr(x)
  b.contrasts.map(row => {
    let mut value = 0.0
    for i in 0.. Array[Double] raise CodaError {
  if z.length() != b.parts - 1 {
    raise Shape("ILR coordinate count mismatch")
  }
  real_vector(z)
  let logits = Array::make(b.parts, 0.0)
  for j, v in z {
    for i in 0..