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

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

///|
pub fn matrix_add(
  first : Array[Array[Double]],
  second : Array[Array[Double]],
) -> Array[Array[Double]] raise VleError {
  if first.length() == 0 {
    raise VleError::EmptyMixture
  }
  ignore(validate_square_matrix(first, first.length()))
  ignore(validate_square_matrix(second, first.length()))
  [
    for i in 0.. {
      [
        for j in 0.. first[i][j] + second[i][j]
      ]
    }
  ]
}

///|
pub fn matrix_subtract(
  first : Array[Array[Double]],
  second : Array[Array[Double]],
) -> Array[Array[Double]] raise VleError {
  if first.length() == 0 {
    raise VleError::EmptyMixture
  }
  ignore(validate_square_matrix(first, first.length()))
  ignore(validate_square_matrix(second, first.length()))
  [
    for i in 0.. {
      [
        for j in 0.. first[i][j] - second[i][j]
      ]
    }
  ]
}

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

///|
pub fn matrix_determinant_2x2(
  matrix : Array[Array[Double]],
) -> Double raise VleError {
  ignore(validate_square_matrix(matrix, 2))
  matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]
}

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

///|
pub fn matrix_identity(dimension : Int) -> Array[Array[Double]] raise VleError {
  if dimension <= 0 {
    raise VleError::InvalidParameter("identity dimension must be positive")
  }
  [
    for i in 0.. {
      [
        for j in 0.. if i == j { 1.0 } else { 0.0 }
      ]
    }
  ]
}

///|
pub fn matrix_quadratic_form(
  matrix : Array[Array[Double]],
  vector : Array[Double],
) -> Double raise VleError {
  let product = matrix_vector_product(matrix, vector)
  for i = 0, value = 0.0; i < vector.length(); {
    continue i + 1, value + vector[i] * product[i]
  } nobreak {
    value
  }
}

///|
pub fn audit_equilibrium_point(point : EquilibriumPoint) -> EquilibriumAudit {
  let liquid_error = abs_double(sum(point.liquid) - 1.0)
  let vapor_error = abs_double(sum(point.vapor) - 1.0)
  let composition_error = if liquid_error > vapor_error {
    liquid_error
  } else {
    vapor_error
  }
  let pressure_error = if point.pressure_bar > 0.0 { 0.0 } else { 1.0 }
  let valid = composition_error <= 1.0E-8 && pressure_error == 0.0
  let messages = if valid {
    []
  } else {
    ["equilibrium point failed normalization or pressure checks"]
  }
  EquilibriumAudit::{
    is_valid: valid,
    maximum_composition_error: composition_error,
    pressure_error,
    material_balance_error: 0.0,
    messages,
  }
}

///|
pub fn audit_flash_result(
  result : FlashResult,
  feed : Array[Double],
) -> EquilibriumAudit raise VleError {
  let composition_error = phase_split_residual(
    feed,
    result.liquid,
    result.vapor,
    result.vapor_fraction,
  )
  let valid = composition_error <= 1.0E-8 &&
    result.pressure_bar > 0.0 &&
    result.vapor_fraction >= 0.0 &&
    result.vapor_fraction <= 1.0
  EquilibriumAudit::{
    is_valid: valid,
    maximum_composition_error: 0.0,
    pressure_error: if result.pressure_bar > 0.0 {
      0.0
    } else {
      1.0
    },
    material_balance_error: composition_error,
    messages: if valid {
      []
    } else {
      ["flash material balance failed"]
    },
  }
}

///|
pub fn audit_activity_flash(
  result : ActivityFlashResult,
  feed : Array[Double],
) -> EquilibriumAudit raise VleError {
  let balance = activity_flash_material_balance_error(result, feed)
  let valid = balance <= 1.0E-8 && result.residual <= 1.0E-6
  EquilibriumAudit::{
    is_valid: valid,
    maximum_composition_error: 0.0,
    pressure_error: if result.pressure_bar > 0.0 {
      0.0
    } else {
      1.0
    },
    material_balance_error: balance,
    messages: if valid {
      []
    } else {
      ["activity flash audit failed"]
    },
  }
}

///|
pub fn equilibrium_point_is_valid(point : EquilibriumPoint) -> Bool {
  audit_equilibrium_point(point).is_valid
}

///|
pub fn audit_error_budget(audit : EquilibriumAudit, tolerance : Double) -> Bool {
  audit.is_valid &&
  audit.maximum_composition_error <= tolerance &&
  audit.material_balance_error <= tolerance
}