{-# OPTIONS_GHC -fglasgow-exts #-} ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Matrix -- Copyright : (c) Alberto Ruiz 2007 -- License : GPL-style -- -- Maintainer : Alberto Ruiz -- Stability : provisional -- Portability : portable -- -- A Matrix representation suitable for numerical computations using LAPACK and GSL. -- ----------------------------------------------------------------------------- module Data.Packed.Matrix ( Field, Matrix,rows,cols, (><), reshape, flatten, fromLists, toLists, (@@>), trans, conjTrans, asRow, asColumn, fromRows, toRows, fromColumns, toColumns, fromBlocks, repmat, flipud, fliprl, subMatrix, takeRows, dropRows, takeColumns, dropColumns, extractRows, ident, diag, diagRect, takeDiag, liftMatrix, liftMatrix2, format, dispR, readMatrix, fromFile, fromArray2D ) where import Data.Packed.Internal import Foreign(Storable) import Complex import Data.Packed.Vector import Numeric(showGFloat) import Data.List(transpose,intersperse) import Data.Array -- | creates a matrix from a vertical list of matrices joinVert :: Field t => [Matrix t] -> Matrix t joinVert ms = case common cols ms of Nothing -> error "joinVert on matrices with different number of columns" Just c -> reshape c $ join (map cdat ms) -- | creates a matrix from a horizontal list of matrices joinHoriz :: Field t => [Matrix t] -> Matrix t joinHoriz ms = trans. joinVert . map trans $ ms {- | Creates a matrix from blocks given as a list of lists of matrices: @\> let a = 'diag' $ 'fromList' [5,7,2] \> let b = 'reshape' 4 $ 'constant' (-1) 12 \> fromBlocks [[a,b],[b,a]] (6><7) [ 5.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0 , 0.0, 7.0, 0.0, -1.0, -1.0, -1.0, -1.0 , 0.0, 0.0, 2.0, -1.0, -1.0, -1.0, -1.0 , -1.0, -1.0, -1.0, -1.0, 5.0, 0.0, 0.0 , -1.0, -1.0, -1.0, -1.0, 0.0, 7.0, 0.0 , -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 2.0 ]@ -} fromBlocks :: Field t => [[Matrix t]] -> Matrix t fromBlocks = joinVert . map joinHoriz -- | Reverse rows flipud :: Field t => Matrix t -> Matrix t flipud m = fromRows . reverse . toRows $ m -- | Reverse columns fliprl :: Field t => Matrix t -> Matrix t fliprl m = fromColumns . reverse . toColumns $ m ------------------------------------------------------------ {- | creates a rectangular diagonal matrix @> diagRect (constant 5 3) 3 4 (3><4) [ 5.0, 0.0, 0.0, 0.0 , 0.0, 5.0, 0.0, 0.0 , 0.0, 0.0, 5.0, 0.0 ]@ -} diagRect :: (Field t, Num t) => Vector t -> Int -> Int -> Matrix t diagRect s r c | dim s < min r c = error "diagRect" | r == c = diag s | r < c = trans $ diagRect s c r | r > c = joinVert [diag s , zeros (r-c,c)] where zeros (r,c) = reshape c $ constantD 0 (r*c) -- | extracts the diagonal from a rectangular matrix takeDiag :: (Field t) => Matrix t -> Vector t takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] -- | creates the identity matrix of given dimension ident :: Field a => Int -> Matrix a ident n = diag (constant 1 n) ------------------------------------------------------------ {- | An easy way to create a matrix: @\> (2><3)[1..6] (2><3) [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 ]@ This is the format produced by the instances of Show (Matrix a), which can also be used for input. -} (><) :: (Field 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 l ---------------------------------------------------------------- -- | Creates a matrix with the first n rows of another matrix takeRows :: Field t => Int -> Matrix t -> Matrix t takeRows n mat = subMatrix (0,0) (n, cols mat) mat -- | Creates a copy of a matrix without the first n rows dropRows :: Field t => Int -> Matrix t -> Matrix t dropRows n mat = subMatrix (n,0) (rows mat - n, cols mat) mat -- |Creates a matrix with the first n columns of another matrix takeColumns :: Field t => Int -> Matrix t -> Matrix t takeColumns n mat = subMatrix (0,0) (rows mat, n) mat -- | Creates a copy of a matrix without the first n columns dropColumns :: Field t => Int -> Matrix t -> Matrix t dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat ---------------------------------------------------------------- {- | Creates a vector by concatenation of rows @\> flatten ('ident' 3) 9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@ -} flatten :: Field t => Matrix t -> Vector t flatten = cdat {- | Creates a 'Matrix' from a list of lists (considered as rows). @\> fromLists [[1,2],[3,4],[5,6]] (3><2) [ 1.0, 2.0 , 3.0, 4.0 , 5.0, 6.0 ]@ -} fromLists :: Field t => [[t]] -> Matrix t fromLists = fromRows . map fromList -- | conjugate transpose conjTrans :: Matrix (Complex Double) -> Matrix (Complex Double) conjTrans = trans . liftMatrix conj -- | creates a 1-row matrix from a vector asRow :: Field a => Vector a -> Matrix a asRow v = reshape (dim v) v -- | creates a 1-column matrix from a vector asColumn :: Field a => Vector a -> Matrix a asColumn v = reshape 1 v ----------------------------------------------------- fromArray2D :: (Field e) => Array (Int, Int) e -> Matrix e fromArray2D m = (r> Int -> a -> String shf dec n | abs n < 1e-10 = "0." | abs (n - (fromIntegral.round $ n)) < 1e-10 = show (round n) ++"." | otherwise = showGFloat (Just dec) n "" -- shows a Complex Double as a pair, with n digits after the decimal point shfc n z@ (a:+b) | magnitude z <1e-10 = "0." | abs b < 1e-10 = shf n a | abs a < 1e-10 = shf n b ++"i" | b > 0 = shf n a ++"+"++shf n b ++"i" | otherwise = shf n a ++shf n b ++"i" dsp' :: String -> [[String]] -> String dsp' sep as = unlines . map unwords' $ transpose mtp where mt = transpose as longs = map (maximum . map length) mt mtp = zipWith (\a b -> map (pad a) b) longs mt pad n str = replicate (n - length str) ' ' ++ str unwords' = concat . intersperse sep format :: (Field t) => String -> (t -> String) -> Matrix t -> String format sep f m = dsp' sep . map (map f) . toLists $ m disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m dispR :: Int -> Matrix Double -> IO () dispR d m = disp m (shf d) dispC :: Int -> Matrix (Complex Double) -> IO () dispC d m = disp m (shfc d) -- | creates a matrix from a table of numbers. readMatrix :: String -> Matrix Double readMatrix = fromLists . map (map read). map words . filter (not.null) . lines -- | rearranges the rows of a matrix according to the order given in a list of integers. extractRows :: Field t => [Int] -> Matrix t -> Matrix t extractRows l m = fromRows $ extract (toRows $ m) l where extract l is = [l!!i |i<-is] {- | creates matrix by repetition of a matrix a given number of rows and columns @> repmat (ident 2) 2 3 :: Matrix Double (4><6) [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]@ -} repmat :: (Field t) => Matrix t -> Int -> Int -> Matrix t repmat m r c = fromBlocks $ partit c $ replicate (r*c) m