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 a. (forall s. ST s a) -> a
runST forall a b. (a -> b) -> a -> b
$ do
let (Int
m,Int
n) = Matrix -> (Int, Int)
dimension Matrix
mat
MMatrix s
r <- forall s. Int -> Int -> Double -> ST s (MMatrix s)
M.replicate Int
n Int
n Double
0
MMatrix s
a <- forall s. Matrix -> ST s (MMatrix s)
M.thaw Matrix
mat
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
for Int
0 Int
n forall a b. (a -> b) -> a -> b
$ \Int
j -> do
Double
cn <- forall a s. NFData a => MMatrix s -> (Matrix -> a) -> ST s a
M.immutably MMatrix s
a forall a b. (a -> b) -> a -> b
$ \Matrix
aa -> Vector -> Double
norm (Matrix -> Int -> Vector
column Matrix
aa Int
j)
forall s. MMatrix s -> Int -> Int -> Double -> ST s ()
M.unsafeWrite MMatrix s
r Int
j Int
j Double
cn
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
for Int
0 Int
m forall a b. (a -> b) -> a -> b
$ \Int
i -> forall s. MMatrix s -> Int -> Int -> (Double -> Double) -> ST s ()
M.unsafeModify MMatrix s
a Int
i Int
j (forall a. Fractional a => a -> a -> a
/ Double
cn)
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
for (Int
jforall a. Num a => a -> a -> a
+Int
1) Int
n forall a b. (a -> b) -> a -> b
$ \Int
jj -> do
Double
p <- forall s. MMatrix s -> Int -> Int -> ST s Double
innerProduct MMatrix s
a Int
j Int
jj
forall s. MMatrix s -> Int -> Int -> Double -> ST s ()
M.unsafeWrite MMatrix s
r Int
j Int
jj Double
p
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
for Int
0 Int
m forall a b. (a -> b) -> a -> b
$ \Int
i -> do
Double
aij <- forall s. MMatrix s -> Int -> Int -> ST s Double
M.unsafeRead MMatrix s
a Int
i Int
j
forall s. MMatrix s -> Int -> Int -> (Double -> Double) -> ST s ()
M.unsafeModify MMatrix s
a Int
i Int
jj forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a -> a
subtract (Double
p forall a. Num a => a -> a -> a
* Double
aij)
(,) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall s. MMatrix s -> ST s Matrix
M.unsafeFreeze MMatrix s
a forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall s. MMatrix s -> ST s Matrix
M.unsafeFreeze MMatrix s
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 = forall a s. NFData a => MMatrix s -> (Matrix -> a) -> ST s a
M.immutably MMatrix s
mmat forall a b. (a -> b) -> a -> b
$ \Matrix
mat ->
forall (v :: * -> *) s.
(Vector v Double, Summation s) =>
(s -> Double) -> v Double -> Double
sumVector KBNSum -> Double
kbn forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
U.zipWith forall a. Num a => a -> a -> a
(*) (Matrix -> Int -> Vector
column Matrix
mat Int
j) (Matrix -> Int -> Vector
column Matrix
mat Int
k)