///|
/// QR decomposition by modified Gram-Schmidt. The returned pair is `(Q, R)`
/// with orthonormal columns in Q and an upper-triangular R.
pub fn Matrix::qr_decompose(self : Matrix) -> (Matrix, Matrix)? {
  if self.rows() < self.cols() || self.is_empty() {
    return None
  }
  let q = Matrix::zeros(self.rows(), self.cols())
  let r = Matrix::zeros(self.cols(), self.cols())
  for column in 0.. ignore
      for row in 0.. ignore
    for row in 0.. ignore
    }
  }
  Some((q, r))
}

///|
/// Least-squares solution using QR decomposition.
pub fn Matrix::least_squares(
  self : Matrix,
  rhs : Array[Double],
) -> MatrixSolveResult {
  if self.rows() < self.cols() || rhs.length() != self.rows() {
    return InvalidShape
  }
  match self.qr_decompose() {
    None => Singular
    Some((q, r)) => {
      let projected = q.transpose().multiply_vector(rhs)
      r.solve(projected)
    }
  }
}

///|
/// Tikhonov-regularized least squares for rank-deficient design matrices.
pub fn Matrix::regularized_least_squares(
  self : Matrix,
  rhs : Array[Double],
  regularization : Double,
) -> MatrixSolveResult {
  if rhs.length() != self.rows() {
    return InvalidShape
  }
  let transpose = self.transpose()
  let normal = transpose
    .multiply(self)
    .add_diagonal(if regularization < 0.0 { 0.0 } else { regularization })
  normal.solve(transpose.multiply_vector(rhs))
}

///|
pub fn Matrix::gram(self : Matrix) -> Matrix {
  self.transpose().multiply(self).symmetric_part()
}

///|
pub fn Matrix::kronecker(self : Matrix, other : Matrix) -> Matrix {
  let result = Matrix::zeros(
    self.rows() * other.rows(),
    self.cols() * other.cols(),
  )
  for i in 0.. ignore
        }
      }
    }
  }
  result
}

///|
pub fn Matrix::horizontal_concat(self : Matrix, other : Matrix) -> Matrix {
  if self.rows() != other.rows() {
    return Matrix::zeros(0, 0)
  }
  let result = Matrix::zeros(self.rows(), self.cols() + other.cols())
  for i in 0.. ignore
    }
    for j in 0.. ignore
    }
  }
  result
}

///|
pub fn Matrix::vertical_concat(self : Matrix, other : Matrix) -> Matrix {
  if self.cols() != other.cols() {
    return Matrix::zeros(0, 0)
  }
  let result = Matrix::zeros(self.rows() + other.rows(), self.cols())
  for i in 0.. ignore
    }
  }
  for i in 0.. ignore
    }
  }
  result
}

///|
pub fn Matrix::row_sums(self : Matrix) -> Array[Double] {
  Array::makei(self.rows(), i => {
    let mut total = 0.0
    for j in 0.. Array[Double] {
  Array::makei(self.cols(), j => {
    let mut total = 0.0
    for i in 0.. Bool {
  if !self.is_square() {
    return false
  }
  for i in 0.. Matrix {
  if !self.is_square() {
    return Matrix::zeros(0, 0)
  }
  let low = if lower < upper { lower } else { upper }
  let high = if lower < upper { upper } else { lower }
  let result = self.copy()
  for i in 0.. high {
        high
      } else {
        value
      },
    )
    |> ignore
  }
  result
}

///|
pub fn Matrix::power(self : Matrix, exponent : Int) -> Matrix {
  if !self.is_square() || exponent < 0 {
    return Matrix::zeros(0, 0)
  }
  let mut remaining = exponent
  let mut base = self.copy()
  let mut result = Matrix::identity(self.rows())
  while remaining > 0 {
    if remaining % 2 == 1 {
      result = result.multiply(base)
    }
    base = base.multiply(base)
    remaining = remaining / 2
  }
  result
}

///|
/// One power-iteration estimate of the dominant eigenpair.
pub struct DominantEigenpair {
  value : Double
  vector : Array[Double]
  iterations : Int
  converged : Bool
} derive(Debug)

///|
pub fn DominantEigenpair::value(self : DominantEigenpair) -> Double {
  self.value
}

///|
pub fn DominantEigenpair::vector(self : DominantEigenpair) -> Array[Double] {
  self.vector.copy()
}

///|
pub fn DominantEigenpair::iterations(self : DominantEigenpair) -> Int {
  self.iterations
}

///|
pub fn DominantEigenpair::converged(self : DominantEigenpair) -> Bool {
  self.converged
}

///|
pub fn Matrix::dominant_eigenpair(
  self : Matrix,
  iterations : Int,
) -> DominantEigenpair {
  if !self.is_square() || self.rows() == 0 {
    return { value: 0.0, vector: [], iterations: 0, converged: false }
  }
  let vector = vector_normalize(
    Array::makei(self.rows(), i => 1.0 + i.to_double()),
  )
  let mut estimate = 0.0
  let mut converged = false
  let limit = if iterations < 1 { 1 } else { iterations }
  let mut completed = 0
  for _ in 0.. Double {
  if !self.is_square() {
    return 0.0
  }
  match self.inverse() {
    None => 0.0
    Some(inverse) => self.one_norm() * inverse.one_norm()
  }
}

///|
pub fn Matrix::finite_or_zero(self : Matrix) -> Matrix {
  self.map(value => if value.is_nan() || value.is_inf() { 0.0 } else { value })
}