///|
pub fn validate_k_values(k_values : Array[Double]) -> Unit raise VleError {
  if k_values.length() == 0 {
    raise VleError::EmptyMixture
  }
  for value in k_values {
    if value <= 0.0 {
      raise VleError::InvalidParameter("K-values must be positive")
    }
  }
}

///|
pub fn classify_phase_regime(
  k_values : Array[Double],
  feed : Array[Double],
) -> PhaseRegime raise VleError {
  validate_k_values(k_values)
  assert_same_length(k_values.length(), feed.length())
  let z = normalize(feed)
  let f0 = rachford_rice(k_values, z, 0.0)
  let f1 = rachford_rice(k_values, z, 1.0)
  if f0 <= 0.0 {
    PhaseRegime::SingleLiquid
  } else if f1 >= 0.0 {
    PhaseRegime::SingleVapor
  } else {
    PhaseRegime::TwoPhase
  }
}

///|
pub fn phase_regime_label(regime : PhaseRegime) -> String {
  match regime {
    SingleLiquid => "SingleLiquid"
    TwoPhase => "TwoPhase"
    SingleVapor => "SingleVapor"
  }
}

///|
pub fn k_value_stability_margin(
  k_values : Array[Double],
  feed : Array[Double],
) -> Double raise VleError {
  validate_k_values(k_values)
  assert_same_length(k_values.length(), feed.length())
  let z = normalize(feed)
  let f0 = abs_double(rachford_rice(k_values, z, 0.0))
  let f1 = abs_double(rachford_rice(k_values, z, 1.0))
  if f0 < f1 {
    f0
  } else {
    f1
  }
}

///|
pub fn k_value_diagnostics(
  k_values : Array[Double],
  feed : Array[Double],
) -> KValueDiagnostics raise VleError {
  validate_k_values(k_values)
  let regime = classify_phase_regime(k_values, feed)
  let z = normalize(feed)
  let f0 = rachford_rice(k_values, z, 0.0)
  let f1 = rachford_rice(k_values, z, 1.0)
  let minimum = for i = 1, value = k_values[0]; i < k_values.length(); {
    continue i + 1, if k_values[i] < value { k_values[i] } else { value }
  } nobreak {
    value
  }
  let maximum = for i = 1, value = k_values[0]; i < k_values.length(); {
    continue i + 1, if k_values[i] > value { k_values[i] } else { value }
  } nobreak {
    value
  }
  KValueDiagnostics::{
    regime,
    f_at_zero: f0,
    f_at_one: f1,
    minimum_k: minimum,
    maximum_k: maximum,
    stability_margin: k_value_stability_margin(k_values, z),
  }
}

///|
pub fn k_value_span(k_values : Array[Double]) -> Double raise VleError {
  validate_k_values(k_values)
  k_value_diagnostics(k_values, [ for _ in k_values => 1.0 ]).maximum_k -
  k_value_diagnostics(k_values, [ for _ in k_values => 1.0 ]).minimum_k
}

///|
pub fn mean_k_value(k_values : Array[Double]) -> Double raise VleError {
  validate_k_values(k_values)
  sum(k_values) / k_values.length().to_double()
}

///|
pub fn phase_split_residual(
  feed : Array[Double],
  liquid : Array[Double],
  vapor : Array[Double],
  vapor_fraction : Double,
) -> Double raise VleError {
  if vapor_fraction < 0.0 || vapor_fraction > 1.0 {
    raise VleError::InvalidParameter("vapor fraction must be bounded")
  }
  assert_same_length(feed.length(), liquid.length())
  assert_same_length(feed.length(), vapor.length())
  let z = normalize(feed)
  let x = normalize(liquid)
  let y = normalize(vapor)
  let reconstructed = [
    for i in 0.. {
      (1.0 - vapor_fraction) * x[i] + vapor_fraction * y[i]
    }
  ]
  vector_distance(z, reconstructed)
}

///|
pub fn phase_split_is_balanced(
  feed : Array[Double],
  liquid : Array[Double],
  vapor : Array[Double],
  vapor_fraction : Double,
  tolerance : Double,
) -> Bool raise VleError {
  if tolerance <= 0.0 {
    raise VleError::InvalidParameter("balance tolerance must be positive")
  }
  phase_split_residual(feed, liquid, vapor, vapor_fraction) <= tolerance
}

///|
pub fn reconstruct_phase_feed(
  liquid : Array[Double],
  vapor : Array[Double],
  vapor_fraction : Double,
) -> Array[Double] raise VleError {
  if vapor_fraction < 0.0 || vapor_fraction > 1.0 {
    raise VleError::InvalidParameter("vapor fraction must be bounded")
  }
  assert_same_length(liquid.length(), vapor.length())
  let x = normalize(liquid)
  let y = normalize(vapor)
  normalize(
    [
      for i in 0.. {
        (1.0 - vapor_fraction) * x[i] + vapor_fraction * y[i]
      }
    ],
  )
}

///|
pub fn k_value_tangent_plane_proxy(
  k_values : Array[Double],
  composition : Array[Double],
) -> Double raise VleError {
  validate_k_values(k_values)
  assert_same_length(k_values.length(), composition.length())
  let z = normalize(composition)
  for i = 0, value = 0.0; i < z.length(); {
    continue i + 1, value + z[i] * @math.ln(clamp_positive(k_values[i]))
  } nobreak {
    value
  }
}

///|
pub fn k_value_is_near_ideal(
  k_values : Array[Double],
  tolerance : Double,
) -> Bool raise VleError {
  validate_k_values(k_values)
  if tolerance <= 0.0 {
    raise VleError::InvalidParameter("K-value tolerance must be positive")
  }
  for value in k_values {
    if abs_double(value - 1.0) > tolerance {
      return false
    }
  }
  true
}

///|
pub fn validate_square_matrix(
  matrix : Array[Array[Double]],
  dimension : Int,
) -> Bool raise VleError {
  if dimension <= 0 || matrix.length() != dimension {
    raise VleError::LengthMismatch(expected=dimension, actual=matrix.length())
  }
  for row in matrix {
    if row.length() != dimension {
      raise VleError::LengthMismatch(expected=dimension, actual=row.length())
    }
  }
  true
}

///|
pub fn matrix_statistics(
  matrix : Array[Array[Double]],
) -> MatrixStatistics raise VleError {
  if matrix.length() == 0 {
    raise VleError::EmptyMixture
  }
  ignore(validate_square_matrix(matrix, matrix.length()))
  let minimum = matrix[0][0]
  let maximum = matrix[0][0]
  let (min_value, max_value, asymmetry, diagonal_error) = for i = 0, minimum = minimum, maximum = maximum, asymmetry = 0.0, diagonal_error = 0.0; i <
                                                             matrix.length(); {
    let (next_min, next_max) = for j = 0, low = minimum, high = maximum; j <
                                  matrix.length(); {
      continue j + 1,
        if matrix[i][j] < low {
          matrix[i][j]
        } else {
          low
        },
        if matrix[i][j] > high {
          matrix[i][j]
        } else {
          high
        }
    } nobreak {
      (low, high)
    }
    let next_asymmetry = for j = 0, value = asymmetry; j < matrix.length(); {
      let difference = abs_double(matrix[i][j] - matrix[j][i])
      continue j + 1, if difference > value { difference } else { value }
    } nobreak {
      value
    }
    let next_diagonal = diagonal_error + abs_double(matrix[i][i] - 1.0)
    continue i + 1, next_min, next_max, next_asymmetry, next_diagonal
  } nobreak {
    (minimum, maximum, asymmetry, diagonal_error)
  }
  MatrixStatistics::{
    dimension: matrix.length(),
    minimum: min_value,
    maximum: max_value,
    maximum_asymmetry: asymmetry,
    diagonal_error,
  }
}

///|
pub fn matrix_trace(matrix : Array[Array[Double]]) -> Double raise VleError {
  if matrix.length() == 0 {
    raise VleError::EmptyMixture
  }
  ignore(validate_square_matrix(matrix, matrix.length()))
  for i = 0, value = 0.0; i < matrix.length(); {
    continue i + 1, value + matrix[i][i]
  } nobreak {
    value
  }
}

///|
pub fn matrix_frobenius_norm(
  matrix : Array[Array[Double]],
) -> Double raise VleError {
  if matrix.length() == 0 {
    raise VleError::EmptyMixture
  }
  ignore(validate_square_matrix(matrix, matrix.length()))
  @math.pow(
    for i = 0, value = 0.0; i < matrix.length(); {
      let row_value = for j = 0, row_sum = 0.0; j < matrix.length(); {
        continue j + 1, row_sum + matrix[i][j] * matrix[i][j]
      } nobreak {
        row_sum
      }
      continue i + 1, value + row_value
    } nobreak {
      value
    },
    0.5,
  )
}

///|
pub fn matrix_row_sums(
  matrix : Array[Array[Double]],
) -> Array[Double] raise VleError {
  if matrix.length() == 0 {
    raise VleError::EmptyMixture
  }
  ignore(validate_square_matrix(matrix, matrix.length()))
  [
    for row in matrix => sum(row)
  ]
}

///|
pub fn matrix_column_sums(
  matrix : Array[Array[Double]],
) -> Array[Double] raise VleError {
  if matrix.length() == 0 {
    raise VleError::EmptyMixture
  }
  ignore(validate_square_matrix(matrix, matrix.length()))
  [
    for j in 0.. {
      for i = 0, value = 0.0; i < matrix.length(); {
        continue i + 1, value + matrix[i][j]
      } nobreak {
        value
      }
    }
  ]
}

///|
pub fn matrix_is_symmetric(
  matrix : Array[Array[Double]],
  tolerance : Double,
) -> Bool raise VleError {
  if tolerance <= 0.0 {
    raise VleError::InvalidParameter("matrix tolerance must be positive")
  }
  matrix_statistics(matrix).maximum_asymmetry <= tolerance
}