///|
pub struct Pca {
  features : Array[String]
  basis : Basis
  location : Array[Double]
  axes : Array[Array[Double]]
  variances : Array[Double]
  total : Double
}

///|
fn coordinate_mean(rows : Array[Array[Double]]) -> Array[Double] {
  let mean = Array::make(rows[0].length(), 0.0)
  for i, row in rows {
    for j in 0.. Unit raise CodaError {
  if actual.length() != expected.length() {
    raise Shape("feature count differs from fitted model")
  }
  for i in 0.. Pca raise CodaError {
  if components < 1 || components >= data.part_count() {
    raise Invalid("PCA components outside 1..D-1")
  }
  let rows = data.ilr_rows(basis)
  let spectrum = symmetric_eigen(covariance(rows, 1))
  let mut total = 0.0
  let values = spectrum.values.map(v => {
    if v < -1.0e-10 {
      raise Numerical("negative covariance eigenvalue")
    }
    let nonnegative = if v < 0.0 { 0.0 } else { v }
    total += nonnegative
    nonnegative
  })
  {
    features: data.features.copy(),
    basis,
    location: coordinate_mean(rows),
    axes: Array::makei(components, i => spectrum.vectors[i].copy()),
    variances: Array::makei(components, i => values[i]),
    total,
  }
}

///|
pub fn Pca::eigenvalues(self : Pca) -> Array[Double] {
  self.variances.copy()
}

///|
pub fn Pca::loadings(self : Pca) -> Array[Array[Double]] {
  self.axes.map(r => r.copy())
}

///|
/// Fractions refer to total pre-truncation variance; constant data returns zeros.
pub fn Pca::explained_ratio(self : Pca) -> Array[Double] {
  self.variances.map(v => if self.total == 0.0 { 0.0 } else { v / self.total })
}

///|
pub fn Pca::center(self : Pca) -> Array[Double] raise CodaError {
  ilr_inverse(self.location, self.basis)
}

///|
pub fn Pca::transform(
  self : Pca,
  data : Dataset,
) -> Array[Array[Double]] raise CodaError {
  matching_features(data.features, self.features)
  data
  .ilr_rows(self.basis)
  .map(row => {
    self.axes.map(axis => {
      let mut score = 0.0
      for j in 0.. Array[Array[Double]] raise CodaError {
  if scores.length() > 1024 {
    raise Budget("PCA reconstruction limited to 1024 rows")
  }
  scores.map(row => {
    if row.length() != self.axes.length() {
      raise Shape("PCA score width mismatch")
    }
    real_vector(row)
    let z = self.location.copy()
    for i in 0..