///|
/// Convert a complete sequential binary partition to orthonormal balances.
/// A row uses +1 and -1 for its two groups and 0 for parts outside that split.
pub fn sbp_basis(partition : Array[Array[Int]]) -> Basis raise CodaError {
  let d = partition.length() + 1
  vector_size(d)
  let rows = []
  for split in partition {
    if split.length() != d {
      raise Shape("SBP must be (D-1) by D")
    }
    let mut plus = 0
    let mut minus = 0
    for sign in split {
      match sign {
        1 => plus += 1
        -1 => minus += 1
        0 => ()
        _ => raise Invalid("SBP signs must be -1, 0 or 1")
      }
    }
    if plus == 0 || minus == 0 {
      raise Invalid("SBP split requires both groups")
    }
    let p = plus.to_double()
    let m = minus.to_double()
    let positive_scale = (m / (p * (p + m))).sqrt()
    let negative_scale = -(p / (m * (p + m))).sqrt()
    rows.push(
      split.map(s => {
        if s == 1 {
          positive_scale
        } else if s == -1 {
          negative_scale
        } else {
          0.0
        }
      }),
    )
  }
  basis(rows)
}

///|
/// One interpretable balance, usable without completing a whole SBP.
pub fn balance(
  x : Array[Double],
  numerator : Array[Int],
  denominator : Array[Int],
) -> Double raise CodaError {
  positive(x)
  if numerator.length() == 0 || denominator.length() == 0 {
    raise Shape("balance groups must be nonempty")
  }
  let seen = Array::make(x.length(), false)
  let mut a = 0.0
  let mut b = 0.0
  for i in numerator {
    if i < 0 || i >= x.length() {
      raise Shape("balance index outside composition")
    }
    if seen[i] {
      raise Invalid("repeated balance index")
    }
    seen[i] = true
    a += @math.ln(x[i]) / numerator.length().to_double()
  }
  for i in denominator {
    if i < 0 || i >= x.length() {
      raise Shape("balance index outside composition")
    }
    if seen[i] {
      raise Invalid("overlapping balance groups")
    }
    seen[i] = true
    b += @math.ln(x[i]) / denominator.length().to_double()
  }
  let p = numerator.length().to_double()
  let m = denominator.length().to_double()
  (p * m / (p + m)).sqrt() * (a - b)
}