///|
/// Gaussian density is with respect to ILR coordinates, NOT simplex volume.
pub struct Gaussian {
  features : Array[String]
  basis : Basis
  location : Array[Double]
  covariance : Array[Array[Double]]
  factor : Cholesky
  ridge : Double
}

///|
/// Maximum-likelihood covariance (ddof=0) plus absolute ILR ridge > 0.
pub fn gaussian(
  data : Dataset,
  ridge : Double,
  basis? : Basis = helmert_basis(data.part_count()),
) -> Gaussian raise CodaError {
  ignore(checked(ridge))
  if ridge <= 0.0 || ridge > 1.0e12 {
    raise Invalid("ridge must be in (0,1e12]")
  }
  let rows = data.ilr_rows(basis)
  let c = covariance(rows, 0)
  for j in 0.. Double {
  self.ridge
}

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

///|
pub fn Gaussian::coordinate_covariance(self : Gaussian) -> Array[Array[Double]] {
  self.covariance.map(r => r.copy())
}

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

///|
pub fn Gaussian::squared_mahalanobis(
  self : Gaussian,
  data : Dataset,
) -> Array[Double] raise CodaError {
  matching_features(data.features, self.features)
  data
  .ilr_rows(self.basis)
  .map(row => {
    let delta = Array::makei(row.length(), j => row[j] - self.location[j])
    self.factor.quadratic(delta)
  })
}

///|
/// Negative log density can be negative; it is not a calibrated anomaly probability.
pub fn Gaussian::negative_log_density(
  self : Gaussian,
  data : Dataset,
) -> Array[Double] raise CodaError {
  let normalizer = self.location.length().to_double() *
    @math.ln(2.0 * 3.141592653589793) +
    self.factor.logdet
  self.squared_mahalanobis(data).map(q => 0.5 * (normalizer + q))
}

///|
/// Rank by descending squared Mahalanobis distance; ties keep input order.
pub fn Gaussian::rank(
  self : Gaussian,
  data : Dataset,
) -> Array[Int] raise CodaError {
  let scores = self.squared_mahalanobis(data)
  let order = Array::makei(scores.length(), i => i)
  for i in 1.. 0 && scores[order[j - 1]] < scores[item] {
      order[j] = order[j - 1]
      j -= 1
    }
    order[j] = item
  }
  order
}