module Statistics.Matrix.Algorithms
(
qr
) where
import Control.Applicative ((<$>), (<*>))
import Control.Monad.ST (ST, runST)
import Prelude hiding (replicate)
import Numeric.Sum (sumVector,kbn)
import Statistics.Matrix (Matrix, column, dimension, for, norm)
import qualified Statistics.Matrix.Mutable as M
import qualified Data.Vector.Unboxed as U
qr :: Matrix -> (Matrix, Matrix)
qr :: Matrix -> (Matrix, Matrix)
qr Matrix
mat = (forall s. ST s (Matrix, Matrix)) -> (Matrix, Matrix)
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Matrix, Matrix)) -> (Matrix, Matrix))
-> (forall s. ST s (Matrix, Matrix)) -> (Matrix, Matrix)
forall a b. (a -> b) -> a -> b
$ do
let (Int
m,Int
n) = Matrix -> (Int, Int)
dimension Matrix
mat
r <- Int -> Int -> Double -> ST s (MMatrix s)
forall s. Int -> Int -> Double -> ST s (MMatrix s)
M.replicate Int
n Int
n Double
0
a <- M.thaw mat
for 0 n $ \Int
j -> do
cn <- MMatrix s -> (Matrix -> Double) -> ST s Double
forall a s. NFData a => MMatrix s -> (Matrix -> a) -> ST s a
M.immutably MMatrix s
a ((Matrix -> Double) -> ST s Double)
-> (Matrix -> Double) -> ST s Double
forall a b. (a -> b) -> a -> b
$ \Matrix
aa -> Vector -> Double
norm (Matrix -> Int -> Vector
column Matrix
aa Int
j)
M.unsafeWrite r j j cn
for 0 m $ \Int
i -> MMatrix s -> Int -> Int -> (Double -> Double) -> ST s ()
forall s. MMatrix s -> Int -> Int -> (Double -> Double) -> ST s ()
M.unsafeModify MMatrix s
a Int
i Int
j (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
cn)
for (j+1) n $ \Int
jj -> do
p <- MMatrix s -> Int -> Int -> ST s Double
forall s. MMatrix s -> Int -> Int -> ST s Double
innerProduct MMatrix s
a Int
j Int
jj
M.unsafeWrite r j jj p
for 0 m $ \Int
i -> do
aij <- MMatrix s -> Int -> Int -> ST s Double
forall s. MMatrix s -> Int -> Int -> ST s Double
M.unsafeRead MMatrix s
a Int
i Int
j
M.unsafeModify a i jj $ subtract (p * aij)
(,) <$> M.unsafeFreeze a <*> M.unsafeFreeze r
innerProduct :: M.MMatrix s -> Int -> Int -> ST s Double
innerProduct :: forall s. MMatrix s -> Int -> Int -> ST s Double
innerProduct MMatrix s
mmat Int
j Int
k = MMatrix s -> (Matrix -> Double) -> ST s Double
forall a s. NFData a => MMatrix s -> (Matrix -> a) -> ST s a
M.immutably MMatrix s
mmat ((Matrix -> Double) -> ST s Double)
-> (Matrix -> Double) -> ST s Double
forall a b. (a -> b) -> a -> b
$ \Matrix
mat ->
(KBNSum -> Double) -> Vector -> Double
forall (v :: * -> *) s.
(Vector v Double, Summation s) =>
(s -> Double) -> v Double -> Double
sumVector KBNSum -> Double
kbn (Vector -> Double) -> Vector -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double) -> Vector -> Vector -> Vector
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
U.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) (Matrix -> Int -> Vector
column Matrix
mat Int
j) (Matrix -> Int -> Vector
column Matrix
mat Int
k)