-- | Compute a Gram-Schmidt orthogonal basis
module Math.LinearAlgebra.GramSchmidt (
    gramSchmidtBasis,
    gramSchmidtOrthogonalization
) where

import           Math.Algebra.LinearAlgebra

-- | Given a basis, return the Gram-Schmidt orhthogonal basis
gramSchmidtBasis :: Fractional a => [[a]] -> [[a]]
gramSchmidtBasis a = fst $ gramSchmidtOrthogonalization a

-- | Given a basis, return the Gram-Schmidt orthogonalization, which is a tuple with the Gram-Schmidt orthogonal basis first, and the
--   $\mu_{i,j} = \langle b_i, b^*_j \rangle / \langle b^*_j, b^*_j \rangle$ triangular matrix second, for $1 \leq j < i < n$.
gramSchmidtOrthogonalization :: Fractional a => [[a]] -> ([[a]], [[a]])
gramSchmidtOrthogonalization (b0:bs) = gs bs [b0] []

-- TODO get rid of the (++) used like this, to make it faster
-- | Perform actual Gram-Schmidt reduction
gs []       b' mu = (b', mu)
gs (b_i:bs) b' mu = gs bs b'' mu'
    where
        mu_i   = mu_row b' b_i
        mu'    = mu ++ [mu_i]
        tosum  = zipWith (*>) mu_i b'
        offset = foldl1 (<+>) tosum
        b'_i   = b_i <-> offset
        b''    = b' ++ [b'_i]

-- | Compute a (partial) row of the $\mu_{i,j}$ matrix. This is based on the previously orthogonalized vectors $b^*_j$, and the current vector $b_i$.
--   This assumes that 'b_i is of length 'i
mu_row b' b_i = flip map b' $ \b'_j -> (b_i <.> b'_j) / (b'_j <.> b'_j)