module Data.Packed.Matrix (
Element,
Matrix,rows,cols,
(><),
trans,
reshape, flatten,
fromLists, toLists, buildMatrix,
(@@>),
asRow, asColumn,
fromRows, toRows, fromColumns, toColumns,
fromBlocks, repmat,
flipud, fliprl,
subMatrix, takeRows, dropRows, takeColumns, dropColumns,
extractRows,
ident, diag, diagRect, takeDiag,
liftMatrix, liftMatrix2, liftMatrix2Auto,
dispf, disps, dispcf, latexFormat, format,
loadMatrix, saveMatrix, fromFile, fileDimensions,
readMatrix, fromArray2D
) where
import Data.Packed.Internal
import qualified Data.Packed.ST as ST
import Data.Packed.Vector
import Data.Array
import System.Process(readProcess)
import Text.Printf(printf)
import Data.List(transpose,intersperse)
import Data.Complex
joinVert :: Element t => [Matrix t] -> Matrix t
joinVert ms = case common cols ms of
Nothing -> error "(impossible) joinVert on matrices with different number of columns"
Just c -> reshape c $ join (map flatten ms)
joinHoriz :: Element t => [Matrix t] -> Matrix t
joinHoriz ms = trans. joinVert . map trans $ ms
fromBlocks :: Element t => [[Matrix t]] -> Matrix t
fromBlocks = fromBlocksRaw . adaptBlocks
fromBlocksRaw mms = joinVert . map joinHoriz $ mms
adaptBlocks ms = ms' where
bc = case common length ms of
Just c -> c
Nothing -> error "fromBlocks requires rectangular [[Matrix]]"
rs = map (compatdim . map rows) ms
cs = map (compatdim . map cols) (transpose ms)
szs = sequence [rs,cs]
ms' = splitEvery bc $ zipWith g szs (concat ms)
g [Just nr,Just nc] m
| nr == r && nc == c = m
| r == 1 && c == 1 = reshape nc (constant x (nr*nc))
| r == 1 = fromRows (replicate nr (flatten m))
| otherwise = fromColumns (replicate nc (flatten m))
where
r = rows m
c = cols m
x = m@@>(0,0)
g _ _ = error "inconsistent dimensions in fromBlocks"
flipud :: Element t => Matrix t -> Matrix t
flipud m = fromRows . reverse . toRows $ m
fliprl :: Element t => Matrix t -> Matrix t
fliprl m = fromColumns . reverse . toColumns $ m
diag :: Element a => Vector a -> Matrix a
diag v = ST.runSTMatrix $ do
let d = dim v
m <- ST.newMatrix 0 d d
mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d1]
return m
diagRect :: (Element t, Num t) => Vector t -> Int -> Int -> Matrix t
diagRect v r c
| dim v < min r c = error "diagRect called with dim v < min r c"
| otherwise = ST.runSTMatrix $ do
m <- ST.newMatrix 0 r c
let d = min r c
mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d1]
return m
takeDiag :: (Element t) => Matrix t -> Vector t
takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) 1]]
ident :: Element a => Int -> Matrix a
ident n = diag (constant 1 n)
(><) :: (Element a) => Int -> Int -> [a] -> Matrix a
r >< c = f where
f l | dim v == r*c = matrixFromVector RowMajor c v
| otherwise = error $ "inconsistent list size = "
++show (dim v) ++" in ("++show r++"><"++show c++")"
where v = fromList $ take (r*c) l
takeRows :: Element t => Int -> Matrix t -> Matrix t
takeRows n mt = subMatrix (0,0) (n, cols mt) mt
dropRows :: Element t => Int -> Matrix t -> Matrix t
dropRows n mt = subMatrix (n,0) (rows mt n, cols mt) mt
takeColumns :: Element t => Int -> Matrix t -> Matrix t
takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
dropColumns :: Element t => Int -> Matrix t -> Matrix t
dropColumns n mt = subMatrix (0,n) (rows mt, cols mt n) mt
fromLists :: Element t => [[t]] -> Matrix t
fromLists = fromRows . map fromList
asRow :: Element a => Vector a -> Matrix a
asRow v = reshape (dim v) v
asColumn :: Element a => Vector a -> Matrix a
asColumn v = reshape 1 v
buildMatrix :: Element a => Int -> Int -> ((Int, Int) -> a) -> Matrix a
buildMatrix rc cc f =
fromLists $ map (\x -> map f x)
$ map (\ ri -> map (\ ci -> (ri, ci)) [0 .. (cc 1)]) [0 .. (rc 1)]
fromArray2D :: (Element e) => Array (Int, Int) e -> Matrix e
fromArray2D m = (r><c) (elems m)
where ((r0,c0),(r1,c1)) = bounds m
r = r1r0+1
c = c1c0+1
format :: (Element t) => String -> (t -> String) -> Matrix t -> String
format sep f m = table sep . map (map f) . toLists $ m
disps :: Int -> Matrix Double -> String
disps d x = sdims x ++ " " ++ formatScaled d x
dispf :: Int -> Matrix Double -> String
dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
sdims x = show (rows x) ++ "x" ++ show (cols x)
formatFixed d x = format " " (printf ("%."++show d++"f")) $ x
isInt = all lookslikeInt . toList . flatten
formatScaled dec t = "E"++show o++"\n" ++ ss
where ss = format " " (printf fmt. g) t
g x | o >= 0 = x/10^(o::Int)
| otherwise = x*10^(o)
o = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
fmt = '%':show (dec+3) ++ '.':show dec ++"f"
vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
vecdisp f v
= ((show (dim v) ++ " |> ") ++) . (++"\n")
. unwords . lines . tail . dropWhile (not . (`elem` " \n"))
. f . trans . reshape 1
$ v
latexFormat :: String
-> String
-> String
latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}"
where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines
showComplex :: Int -> Complex Double -> String
showComplex d (a:+b)
| isZero a && isZero b = "0"
| isZero b = sa
| isZero a && isOne b = s2++"i"
| isZero a = sb++"i"
| isOne b = sa++s3++"i"
| otherwise = sa++s1++sb++"i"
where
sa = shcr d a
sb = shcr d b
s1 = if b<0 then "" else "+"
s2 = if b<0 then "-" else ""
s3 = if b<0 then "-" else "+"
shcr d a | lookslikeInt a = printf "%.0f" a
| otherwise = printf ("%."++show d++"f") a
lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
where shx = show x
isZero x = show x `elem` ["0.0","-0.0"]
isOne x = show x `elem` ["1.0","-1.0"]
dispcf :: Int -> Matrix (Complex Double) -> String
dispcf d m = sdims m ++ "\n" ++ format " " (showComplex d) m
readMatrix :: String -> Matrix Double
readMatrix = fromLists . map (map read). map words . filter (not.null) . lines
fileDimensions :: FilePath -> IO (Int,Int)
fileDimensions fname = do
wcres <- readProcess "wc" ["-w",fname] ""
contents <- readFile fname
let tot = read . head . words $ wcres
c = length . head . dropWhile null . map words . lines $ contents
if tot > 0
then return (tot `div` c, c)
else return (0,0)
loadMatrix :: FilePath -> IO (Matrix Double)
loadMatrix file = fromFile file =<< fileDimensions file
fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
extractRows :: Element t => [Int] -> Matrix t -> Matrix t
extractRows l m = fromRows $ extract (toRows $ m) l
where extract l' is = [l'!!i |i<-is]
repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
repmat m r c = fromBlocks $ splitEvery c $ replicate (r*c) m
liftMatrix2Auto :: (Element t, Element a, Element b)
=> (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
liftMatrix2Auto f m1 m2 | compat' m1 m2 = lM f m1 m2
| rows m1 == rows m2 && cols m2 == 1 = lM f m1 (repCols (cols m1) m2)
| rows m1 == rows m2 && cols m1 == 1 = lM f (repCols (cols m2) m1) m2
| cols m1 == cols m2 && rows m2 == 1 = lM f m1 (repRows (rows m1) m2)
| cols m1 == cols m2 && cols m1 == 1 = lM f (repRows (rows m2) m1) m2
| rows m1 == 1 && cols m2 == 1 = lM f (repRows (rows m2) m1) (repCols (cols m1) m2)
| cols m1 == 1 && rows m2 == 1 = lM f (repCols (cols m2) m1) (repRows (rows m1) m2)
| otherwise = error $ "nonconformable matrices in liftMatrix2Auto: " ++ show (size m1) ++ ", " ++ show (size m2)
size m = (rows m, cols m)
lM f m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))
repRows n x = fromRows (replicate n (flatten x))
repCols n x = fromColumns (replicate n (flatten x))
compat' :: Matrix a -> Matrix b -> Bool
compat' m1 m2 = rows m1 == 1 && cols m1 == 1
|| rows m2 == 1 && cols m2 == 1
|| rows m1 == rows m2 && cols m1 == cols m2