module Data.Matrix.Dense.Generic
(
Matrix(..)
, MG.dim
, MG.rows
, MG.cols
, MG.unsafeIndex
, (MG.!)
, MG.takeRow
, MG.takeColumn
, MG.takeDiag
, MG.unsafeFromVector
, MG.fromVector
, MG.matrix
, MG.fromLists
, MG.fromRows
, fromColumns
, MG.empty
, MG.flatten
, MG.toRows
, MG.toColumns
, MG.toList
, MG.toLists
, convert
, tr
, subMatrix
, ident
, diag
, diagRect
, fromBlocks
, isSymmetric
, force
, Data.Matrix.Dense.Generic.foldl
, imap
, Data.Matrix.Dense.Generic.map
, mapM
, mapM_
, forM
, forM_
, Data.Matrix.Dense.Generic.zipWith
, Data.Matrix.Dense.Generic.zipWith3
, zipWith4
, zipWith5
, zipWith6
, izipWith
, izipWith3
, izipWith4
, izipWith5
, izipWith6
, Data.Matrix.Dense.Generic.zip
, Data.Matrix.Dense.Generic.zip3
, zip4
, zip5
, zip6
, zipWithM
, zipWithM_
, Data.Matrix.Dense.Generic.unzip
, Data.Matrix.Dense.Generic.unzip3
, unzip4
, unzip5
, unzip6
, Data.Matrix.Dense.Generic.sequence
, Data.Matrix.Dense.Generic.sequence_
, generate
, MG.thaw
, MG.unsafeThaw
, MG.freeze
, MG.unsafeFreeze
, MG.create
) where
import Prelude hiding (mapM_, mapM)
import Control.Arrow ((***), (&&&))
import Control.Monad (liftM, foldM, foldM_)
import qualified Data.Foldable as F
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as GM
import qualified Data.Matrix.Generic as MG
import Data.Matrix.Dense.Generic.Mutable (MMatrix(..))
type instance MG.Mutable Matrix = MMatrix
data Matrix v a = Matrix !Int
!Int
!Int
!Int
!(v a)
deriving (Show)
instance G.Vector v a => MG.Matrix Matrix v a where
dim (Matrix r c _ _ _) = (r,c)
unsafeIndex (Matrix _ _ tda offset vec) (i,j) = vec `G.unsafeIndex` idx
where
idx = offset + i * tda + j
unsafeFromVector (r,c) = Matrix r c c 0
unsafeTakeRow (Matrix _ c tda offset vec) i = G.slice i' c vec
where
i' = offset + i * tda
flatten (Matrix r c tda offset vec)
| c == tda = G.slice offset (r*c) vec
| otherwise = G.generate (r*c) $ \i ->
vec `G.unsafeIndex` (offset + (i `div` c) * tda + (i `mod` c))
thaw (Matrix r c tda offset v) = MMatrix r c tda offset `liftM` G.thaw v
unsafeThaw (Matrix r c tda offset v) = MMatrix r c tda offset `liftM` G.unsafeThaw v
freeze (MMatrix r c tda offset v) = Matrix r c tda offset `liftM` G.freeze v
unsafeFreeze (MMatrix r c tda offset v) = Matrix r c tda offset `liftM` G.unsafeFreeze v
fromColumns :: G.Vector v a => [v a] -> Matrix v a
fromColumns = tr . MG.fromRows
convert :: (G.Vector v a, G.Vector w a) => Matrix v a -> Matrix w a
convert (Matrix r c tda offset vec) = Matrix r c tda offset . G.convert $ vec
subMatrix :: G.Vector v a
=> (Int, Int)
-> (Int, Int)
-> Matrix v a -> Matrix v a
subMatrix (i,j) (i',j') (Matrix _ n tda offset vec)
| m' <= 0 || n' <= 0 = MG.empty
| otherwise = Matrix m' n' tda offset' vec
where
m' = i' i + 1
n' = j' j + 1
offset' = offset + i * n + j
tr :: G.Vector v a => Matrix v a -> Matrix v a
tr (Matrix r c tda offset vec) = MG.fromVector (c,r) $ G.generate (r*c) f
where
f i = vec G.! (offset + i `mod` r * tda + i `div` r)
ident :: (Num a, G.Vector v a) => Int -> Matrix v a
ident n = diagRect 0 (n,n) $ replicate n 1
diag :: (Num a, G.Vector v a, F.Foldable t)
=> t a
-> Matrix v a
diag d = diagRect 0 (n,n) d
where n = length . F.toList $ d
diagRect :: (G.Vector v a, F.Foldable t)
=> a
-> (Int, Int)
-> t a
-> Matrix v a
diagRect z0 (r,c) d = MG.fromVector (r,c) $ G.create $ GM.replicate n z0 >>= go d c
where
go xs c' v = F.foldlM f 0 xs >> return v
where
f !i x = GM.unsafeWrite v (i*(c'+1)) x >> return (i+1)
n = r * c
fromBlocks :: G.Vector v a
=> a
-> [[Matrix v a]]
-> Matrix v a
fromBlocks d ms = MG.fromVector (m,n) $ G.create $ GM.replicate (m*n) d >>= go n ms
where
go n' xss v = foldM_ f 0 xss >> return v
where
f !cr xs = do (r', _) <- foldM g (0, 0) xs
return $ cr + r'
where
g (!maxR, !cc) x = do
let (r,c) = MG.dim x
vec = MG.flatten x
step i u = do
GM.unsafeWrite v ((cr + i `div` c) * n' + i `mod` c + cc) u
return (i+1)
G.foldM'_ step (0::Int) vec
return (max maxR r, cc + c)
(m, n) = (sum *** maximum) . Prelude.unzip . Prelude.map ((maximum *** sum) .
Prelude.unzip . Prelude.map (MG.rows &&& MG.cols)) $ ms
isSymmetric :: (Eq a, G.Vector v a) => Matrix v a -> Bool
isSymmetric m@(Matrix r c _ _ _) | r /= c = False
| otherwise = all f [0 .. r1]
where
f i = all g [i + 1 .. c1]
where g j = m MG.! (i,j) == m MG.! (j,i)
force :: G.Vector v a => Matrix v a -> Matrix v a
force m@(Matrix r c _ _ _) = MG.fromVector (r,c) . G.force . MG.flatten $ m
imap :: (G.Vector v a, G.Vector v b) => ((Int, Int) -> a -> b) -> Matrix v a -> Matrix v b
imap f m@(Matrix r c _ _ _) = MG.fromVector (r,c) $ G.imap f' . MG.flatten $ m
where
f' i = f (i `div` c, i `mod` c)
map :: (G.Vector v a, G.Vector v b) => (a -> b) -> Matrix v a -> Matrix v b
map f m@(Matrix r c _ _ _) = MG.fromVector (r,c) $ G.map f . MG.flatten $ m
foldl :: G.Vector v b => (a -> b -> a) -> a -> Matrix v b -> a
foldl f acc m = G.foldl f acc . MG.flatten $ m
mapM :: (G.Vector v a, G.Vector v b, Monad m) => (a -> m b) -> Matrix v a -> m (Matrix v b)
mapM f m@(Matrix r c _ _ _) = liftM (MG.fromVector (r,c)) . G.mapM f . MG.flatten $ m
mapM_ :: (G.Vector v a, Monad m) => (a -> m b) -> Matrix v a -> m ()
mapM_ f = G.mapM_ f . MG.flatten
forM :: (G.Vector v a, G.Vector v b, Monad m) => Matrix v a -> (a -> m b) -> m (Matrix v b)
forM = flip mapM
forM_ :: (G.Vector v a, Monad m) => Matrix v a -> (a -> m b) -> m ()
forM_ = flip mapM_
zipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)
=> (a -> b -> c) -> Matrix v a -> Matrix v b -> Matrix v c
zipWith f m1 m2
| MG.dim m1 /= MG.dim m2 = error "zipWith: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zipWith f (MG.flatten m1) $ MG.flatten m2
zipWith3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d)
=> (a -> b -> c -> d) -> Matrix v a -> Matrix v b -> Matrix v c
-> Matrix v d
zipWith3 f m1 m2 m3
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 = error "zipWith3: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zipWith3 f (MG.flatten m1) (MG.flatten m2) $ MG.flatten m3
zipWith4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v e)
=> (a -> b -> c -> d -> e) -> Matrix v a -> Matrix v b -> Matrix v c
-> Matrix v d -> Matrix v e
zipWith4 f m1 m2 m3 m4
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 = error "zipWith4: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zipWith4 f (MG.flatten m1) (MG.flatten m2)
(MG.flatten m3) $ MG.flatten m4
zipWith5 :: ( G.Vector v a, G.Vector v b, G.Vector v c,G.Vector v d
, G.Vector v e, G.Vector v f )
=> (a -> b -> c -> d -> e -> f) -> Matrix v a -> Matrix v b
-> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
zipWith5 f m1 m2 m3 m4 m5
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 ||
MG.dim m4 /= MG.dim m5 = error "zipWith5: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zipWith5 f (MG.flatten m1) (MG.flatten m2)
(MG.flatten m3) (MG.flatten m4) $ MG.flatten m5
zipWith6 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
, G.Vector v e, G.Vector v f, G.Vector v g )
=> (a -> b -> c -> d -> e -> f -> g) -> Matrix v a -> Matrix v b
-> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f -> Matrix v g
zipWith6 f m1 m2 m3 m4 m5 m6
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 ||
MG.dim m4 /= MG.dim m5 ||
MG.dim m5 /= MG.dim m6 = error "zipWith6: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zipWith6 f (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
(MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
izipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)
=> ((Int, Int) -> a -> b -> c) -> Matrix v a -> Matrix v b -> Matrix v c
izipWith f m1 m2
| MG.dim m1 /= MG.dim m2 = error "izipWith: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.izipWith f' (MG.flatten m1) $ MG.flatten m2
where
c = MG.cols m1
f' i = f (i `div` c, i `mod` c)
izipWith3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d)
=> ((Int, Int) -> a -> b -> c -> d) -> Matrix v a -> Matrix v b
-> Matrix v c -> Matrix v d
izipWith3 f m1 m2 m3
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 = error "izipWith3: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.izipWith3 f' (MG.flatten m1) (MG.flatten m2) $ MG.flatten m3
where
c = MG.cols m1
f' i = f (i `div` c, i `mod` c)
izipWith4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v e)
=> ((Int, Int) -> a -> b -> c -> d -> e) -> Matrix v a -> Matrix v b
-> Matrix v c -> Matrix v d -> Matrix v e
izipWith4 f m1 m2 m3 m4
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 = error "izipWith4: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.izipWith4 f' (MG.flatten m1) (MG.flatten m2)
(MG.flatten m3) $ MG.flatten m4
where
c = MG.cols m1
f' i = f (i `div` c, i `mod` c)
izipWith5 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
, G.Vector v e, G.Vector v f )
=> ((Int, Int) -> a -> b -> c -> d -> e -> f) -> Matrix v a
-> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
izipWith5 f m1 m2 m3 m4 m5
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 ||
MG.dim m4 /= MG.dim m5 = error "izipWith5: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.izipWith5 f' (MG.flatten m1) (MG.flatten m2)
(MG.flatten m3) (MG.flatten m4) $ MG.flatten m5
where
c = MG.cols m1
f' i = f (i `div` c, i `mod` c)
izipWith6 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
, G.Vector v e, G.Vector v f, G.Vector v g )
=> ((Int, Int) -> a -> b -> c -> d -> e -> f -> g) -> Matrix v a
-> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
-> Matrix v g
izipWith6 f m1 m2 m3 m4 m5 m6
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 ||
MG.dim m4 /= MG.dim m5 ||
MG.dim m5 /= MG.dim m6 = error "izipWith6: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.izipWith6 f' (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
(MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
where
c = MG.cols m1
f' i = f (i `div` c, i `mod` c)
zip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))
=> Matrix v a -> Matrix v b -> Matrix v (a,b)
zip m1 m2
| MG.dim m1 /= MG.dim m2 = error "zip: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zip (MG.flatten m1) $ MG.flatten m2
zip3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c))
=> Matrix v a -> Matrix v b -> Matrix v c -> Matrix v (a,b,c)
zip3 m1 m2 m3
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 = error "zip3: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zip3 (MG.flatten m1) (MG.flatten m2) $ MG.flatten m3
zip4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v (a,b,c,d))
=> Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v (a,b,c,d)
zip4 m1 m2 m3 m4
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 = error "zip4: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zip4 (MG.flatten m1) (MG.flatten m2)
(MG.flatten m3) $ MG.flatten m4
zip5 :: ( G.Vector v a, G.Vector v b, G.Vector v c
, G.Vector v d, G.Vector v e, G.Vector v (a,b,c,d,e) )
=> Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e
-> Matrix v (a,b,c,d,e)
zip5 m1 m2 m3 m4 m5
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 ||
MG.dim m4 /= MG.dim m5 = error "zip5: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zip5 (MG.flatten m1) (MG.flatten m2)
(MG.flatten m3) (MG.flatten m4) $ MG.flatten m5
zip6 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v e
, G.Vector v f, G.Vector v (a,b,c,d,e,f) )
=> Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e
-> Matrix v f -> Matrix v (a,b,c,d,e,f)
zip6 m1 m2 m3 m4 m5 m6
| MG.dim m1 /= MG.dim m2 ||
MG.dim m2 /= MG.dim m3 ||
MG.dim m3 /= MG.dim m4 ||
MG.dim m4 /= MG.dim m5 ||
MG.dim m5 /= MG.dim m6 = error "zip6: Dimensions don't match."
| otherwise = MG.fromVector (MG.dim m1) $
G.zip6 (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
(MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
zipWithM :: (Monad m, G.Vector v a, G.Vector v b, G.Vector v c)
=> (a -> b -> m c) -> Matrix v a -> Matrix v b -> m (Matrix v c)
zipWithM f m1 m2
| MG.dim m1 /= MG.dim m2 = error "zipWithM: Dimensions don't match."
| otherwise = liftM (MG.fromVector $ MG.dim m1) $
G.zipWithM f (MG.flatten m1) $ MG.flatten m2
zipWithM_ :: (Monad m, G.Vector v a, G.Vector v b)
=> (a -> b -> m c) -> Matrix v a -> Matrix v b -> m ()
zipWithM_ f m1 m2
| MG.dim m1 /= MG.dim m2 = error "zipWithM_: Dimensions don't match."
| otherwise = G.zipWithM_ f (MG.flatten m1) $ MG.flatten m2
unzip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))
=> Matrix v (a,b) -> (Matrix v a, Matrix v b )
unzip m = (MG.fromVector d v1, MG.fromVector d v2)
where
d = MG.dim m
(v1, v2) = G.unzip $ MG.flatten m
unzip3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c))
=> Matrix v (a,b, c) -> (Matrix v a, Matrix v b, Matrix v c)
unzip3 m = (MG.fromVector d v1, MG.fromVector d v2, MG.fromVector d v3)
where
d = MG.dim m
(v1, v2, v3) = G.unzip3 $ MG.flatten m
unzip4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v (a,b,c,d))
=> Matrix v (a,b,c,d) -> (Matrix v a, Matrix v b, Matrix v c, Matrix v d)
unzip4 m = ( MG.fromVector d v1
, MG.fromVector d v2
, MG.fromVector d v3
, MG.fromVector d v4
)
where
d = MG.dim m
(v1, v2, v3, v4) = G.unzip4 $ MG.flatten m
unzip5 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
, G.Vector v e, G.Vector v (a,b,c,d,e) )
=> Matrix v (a,b,c,d,e)
-> (Matrix v a, Matrix v b, Matrix v c, Matrix v d, Matrix v e)
unzip5 m = ( MG.fromVector d v1
, MG.fromVector d v2
, MG.fromVector d v3
, MG.fromVector d v4
, MG.fromVector d v5
)
where
d = MG.dim m
(v1, v2, v3, v4, v5) = G.unzip5 $ MG.flatten m
unzip6 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
, G.Vector v e, G.Vector v f, G.Vector v (a,b,c,d,e,f) )
=> Matrix v (a,b,c,d,e,f)
-> (Matrix v a, Matrix v b, Matrix v c, Matrix v d, Matrix v e, Matrix v f)
unzip6 m = ( MG.fromVector d v1
, MG.fromVector d v2
, MG.fromVector d v3
, MG.fromVector d v4
, MG.fromVector d v5
, MG.fromVector d v6
)
where
d = MG.dim m
(v1, v2, v3, v4, v5, v6) = G.unzip6 $ MG.flatten m
sequence :: (G.Vector v a, G.Vector v (m a), Monad m)
=> Matrix v (m a) -> m (Matrix v a)
sequence (Matrix r c tda offset vec) = liftM (Matrix r c tda offset) . G.sequence $ vec
sequence_ :: (G.Vector v (m a), Monad m)
=> Matrix v (m a) -> m ()
sequence_ (Matrix _ _ _ _ vec) = G.sequence_ vec
generate :: G.Vector v a => (Int, Int) -> ((Int, Int) -> a) -> Matrix v a
generate (r,c) f = MG.fromVector (r,c) . G.generate (r*c) $ \i -> f (i `div` c, i `mod` c)