matrix 0.2.4.0 → 0.3.0.0
raw patch · 4 files changed
+220/−58 lines, 4 filesdep +QuickCheckdep +tastydep +tasty-quickcheck
Dependencies added: QuickCheck, tasty, tasty-quickcheck
Files
- Data/Matrix.hs +104/−54
- license +1/−1
- matrix.cabal +13/−3
- test/Main.hs +102/−0
Data/Matrix.hs view
@@ -1,6 +1,8 @@+ -- | Matrix datatype and operations. -- -- Every provided example has been tested.+-- Run @cabal test@ for further tests. module Data.Matrix ( -- * Matrix type Matrix , prettyMatrix@@ -19,18 +21,21 @@ , getElem , (!) , safeGet , getRow , getCol , getDiag+ , getMatrixAsVector -- * Manipulating matrices , setElem- , transpose , extendTo+ , transpose , setSize , extendTo , mapRow , mapCol -- * Submatrices -- ** Splitting blocks , submatrix , minorMatrix , splitBlocks- -- ** Joining blocks+ -- ** Joining blocks , (<|>) , (<->) , joinBlocks+ -- * Matrix operations+ , elementwise -- * Matrix multiplication -- ** About matrix multiplication -- $mult@@ -46,8 +51,8 @@ , switchRows , switchCols -- * Decompositions- , luDecomp- , luDecomp'+ , luDecomp , luDecompUnsafe+ , luDecomp', luDecompUnsafe' , cholDecomp -- * Properties , trace , diagProd@@ -64,7 +69,7 @@ import Data.Traversable -- Data import Control.Monad.Primitive (PrimMonad, PrimState)-import Data.List (maximumBy)+import Data.List (maximumBy,foldl1') import Data.Ord (comparing) import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV@@ -87,7 +92,7 @@ data Matrix a = M { nrows :: {-# UNPACK #-} !Int -- ^ Number of rows. , ncols :: {-# UNPACK #-} !Int -- ^ Number of columns.- , mvect :: (V.Vector a) -- ^ Content of the matrix as a plain vector.+ , mvect :: V.Vector a -- ^ Content of the matrix as a plain vector. } deriving Eq -- | Just a cool way to output the size of a matrix.@@ -119,8 +124,12 @@ ---- FUNCTOR INSTANCE instance Functor Matrix where+ {-# INLINE fmap #-} fmap f (M n m v) = M n m $ V.map f v +-------------------------------------------------------+-------------------------------------------------------+ -- | /O(rows*cols)/. Map a function over a row. -- Example: --@@ -218,7 +227,7 @@ -> [a] -- ^ List of elements -> Matrix a {-# INLINE fromList #-}-fromList n m = M n m . V.fromList+fromList n m = M n m . V.fromListN (n*m) -- | Create a matrix from an non-empty list of non-empty lists. -- /Each list must have the same number of elements/.@@ -280,7 +289,7 @@ -> Matrix a -- ^ Matrix -> a {-# INLINE getElem #-}-getElem i j (M _ m v) = v V.! encode m (i,j)+getElem i j (M _ m v) = V.unsafeIndex v $ encode m (i,j) -- | Short alias for 'getElem'. {-# INLINE (!) #-}@@ -290,12 +299,12 @@ -- | Safe variant of 'getElem'. safeGet :: Int -> Int -> Matrix a -> Maybe a safeGet i j a@(M n m _)- | i > n || j > m = Nothing+ | i > n || j > m || i < 1 || j < 1 = Nothing | otherwise = Just $ getElem i j a --- | /O(cols)/. Get a row of a matrix as a vector.+-- | /O(1)/. Get a row of a matrix as a vector. getRow :: Int -> Matrix a -> V.Vector a-getRow i (M _ m v) = V.generate m $ \j -> v V.! encode m (i,j+1)+getRow i (M _ m v) = V.slice (m*(i-1)) m v -- | /O(rows)/. Get a column of a matrix as a vector. getCol :: Int -> Matrix a -> V.Vector a@@ -307,6 +316,12 @@ where k = min (nrows m) (ncols m) +-- | /O(1)/. Transform a 'Matrix' to a 'V.Vector' of size /rows*cols/.+-- This is equivalent to get all the rows of the matrix using 'getRow'+-- and then append them, but far more efficient.+getMatrixAsVector :: Matrix a -> V.Vector a+getMatrixAsVector = mvect+ ------------------------------------------------------- ------------------------------------------------------- ---- MANIPULATING MATRICES@@ -330,26 +345,36 @@ transpose :: Matrix a -> Matrix a transpose m = matrix (ncols m) (nrows m) $ \(i,j) -> m ! (j,i) --- | Extend a matrix to a given size adding zeroes.+-- | Extend a matrix to a given size adding a default element. -- If the matrix already has the required size, nothing happens. -- The matrix is /never/ reduced in size. -- Example: ----- > ( 1 2 3 0 0 )--- > ( 1 2 3 ) ( 4 5 6 0 0 )--- > ( 4 5 6 ) ( 7 8 9 0 0 )--- > extendTo 4 5 ( 7 8 9 ) = ( 0 0 0 0 0 )-extendTo :: Num a- => Int -- ^ Minimal number of rows.+-- > ( 1 2 3 0 0 )+-- > ( 1 2 3 ) ( 4 5 6 0 0 )+-- > ( 4 5 6 ) ( 7 8 9 0 0 )+-- > extendTo 0 4 5 ( 7 8 9 ) = ( 0 0 0 0 0 )+extendTo :: a -- ^ Element to add when extending.+ -> Int -- ^ Minimal number of rows. -> Int -- ^ Minimal number of columns. -> Matrix a -> Matrix a-extendTo n m a = a''- where- n' = n - nrows a- a' = if n' <= 0 then a else a <-> zero n' (ncols a)- m' = m - ncols a- a'' = if m' <= 0 then a' else a' <|> zero (nrows a') m'+extendTo e n m a = setSize e (max n $ nrows a) (max m $ ncols a) a +-- | Set the size of a matrix to given parameters. Use a default element+-- for undefined entries if the matrix has been extended.+setSize :: a -- ^ Default element.+ -> Int -- ^ Number of rows.+ -> Int -- ^ Number of columns.+ -> Matrix a+ -> Matrix a+setSize e n m a@(M r c _) = M n m $ V.generate (n*m) $+ \k -> let (i,j) = d k+ in if i <= r && j <= c+ then a ! (i,j)+ else e+ where+ d = decode m+ ------------------------------------------------------- ------------------------------------------------------- ---- WORKING WITH BLOCKS@@ -411,15 +436,12 @@ -> Matrix a -- ^ Matrix to split. -> (Matrix a,Matrix a ,Matrix a,Matrix a) -- ^ (TL,TR,BL,BR)-{-# INLINE splitBlocks #-} splitBlocks i j a@(M n m _) = ( submatrix 1 i 1 j a , submatrix 1 i (j+1) m a , submatrix (i+1) n 1 j a , submatrix (i+1) n (j+1) m a ) + -- | Join blocks of the form detailed in 'splitBlocks'.-joinBlocks :: (Matrix a,Matrix a- ,Matrix a,Matrix a)- -> Matrix a-{-# INLINE joinBlocks #-}+joinBlocks :: (Matrix a,Matrix a,Matrix a,Matrix a) -> Matrix a joinBlocks (tl,tr,bl,br) = (tl <|> tr) <-> (bl <|> br)@@ -458,6 +480,15 @@ ------------------------------------------------------- -------------------------------------------------------+---- MATRIX OPERATIONS++-- | Perform an operation elementwise. The input matrices are assumed+-- to have the same dimensions, but this is not checked.+elementwise :: (a -> b -> c) -> (Matrix a -> Matrix b -> Matrix c)+elementwise f (M n m v) (M _ _ v') = M n m $ V.zipWith f v v'++-------------------------------------------------------+------------------------------------------------------- ---- MATRIX MULTIPLICATION {- $mult@@ -542,12 +573,12 @@ | otherwise = let mx = maximum [n,m,n',m'] n2 = first (>= mx) $ fmap (2^) [(0 :: Int)..]- b1 = extendTo n2 n2 a1- b2 = extendTo n2 n2 a2+ b1 = setSize 0 n2 n2 a1+ b2 = setSize 0 n2 n2 a2 in submatrix 1 n 1 m' $ strassen b1 b2 strmixFactor :: Int-strmixFactor = 100+strmixFactor = 2 ^ (6 :: Int) -- | Strassen's mixed algorithm. strassenMixed :: Num a => Matrix a -> Matrix a -> Matrix a@@ -556,8 +587,8 @@ strassenMixed a@(M r _ _) b | r < strmixFactor = multStd_ a b | odd r = let r' = r + 1- a' = extendTo r' r' a- b' = extendTo r' r' b+ a' = setSize 0 r' r' a+ b' = setSize 0 r' r' b in submatrix 1 r 1 r $ strassenMixed a' b' | otherwise = joinBlocks (c11,c12,c21,c22) where@@ -590,8 +621,8 @@ | otherwise = let mx = maximum [n,m,n',m'] n2 = if even mx then mx else mx+1- b1 = extendTo n2 n2 a1- b2 = extendTo n2 n2 a2+ b1 = setSize 0 n2 n2 a1+ b2 = setSize 0 n2 n2 a2 in submatrix 1 n 1 m' $ strassenMixed b1 b2 -------------------------------------------------------@@ -719,10 +750,10 @@ -- > ( 1 2 0 ) ( 2 0 2 ) ( 1 0 0 ) ( 0 0 1 ) -- > ( 0 2 1 ) ( 0 2 -1 ) ( 1/2 1 0 ) ( 1 0 0 ) -- > luDecomp ( 2 0 2 ) = ( ( 0 0 2 ) , ( 0 1 1 ) , ( 0 1 0 ) , 1 )-luDecomp :: (Ord a, Fractional a) => Matrix a -> (Matrix a,Matrix a,Matrix a,a)+luDecomp :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a,Matrix a,Matrix a,a) luDecomp a = recLUDecomp a i i 1 1 n where- i = (identity $ nrows a)+ i = identity $ nrows a n = min (nrows a) (ncols a) recLUDecomp :: (Ord a, Fractional a)@@ -732,10 +763,11 @@ -> a -- ^ d -> Int -- ^ Current row -> Int -- ^ Total rows- -> (Matrix a,Matrix a,Matrix a,a)+ -> Maybe (Matrix a,Matrix a,Matrix a,a) recLUDecomp u l p d k n =- if k > n then (u,l,p,d)- else recLUDecomp u'' l'' p' d' (k+1) n+ if k > n then Just (u,l,p,d)+ else if ukk == 0 then Nothing+ else recLUDecomp u'' l'' p' d' (k+1) n where -- Pivot strategy: maximum value in absolute value below the current row. i = maximumBy (\x y -> compare (abs $ u ! (x,k)) (abs $ u ! (y,k))) [ k .. n ]@@ -758,6 +790,12 @@ else let x = (u_ ! (j,k)) / ukk in go (combineRows j (-x) k u_) (setElem x (j,k) l_) (j+1) +-- | Unsafe version of 'luDecomp'. It fails when the input matrix is singular.+luDecompUnsafe :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, a)+luDecompUnsafe m = case luDecomp m of+ Just x -> x+ _ -> error "luDecompUnsafe of singular matrix."+ -- | Matrix LU decomposition with /complete pivoting/. -- The result for a matrix /M/ is given in the format /(U,L,P,Q,d,e)/ where: --@@ -765,9 +803,9 @@ -- -- * /L/ is an /unit/ lower triangular matrix. ----- * /P,Q/ is a permutation matrix.+-- * /P,Q/ are permutation matrices. ----- * /d,e/ is the determinant of /P,Q/.+-- * /d,e/ are the determinants of /P/ and /Q/ respectively. -- -- * /PMQ = LU/. --@@ -784,12 +822,18 @@ -- > ( 1 0 ) ( 2 1 ) ( 1 0 0 ) ( 0 0 1 ) -- > ( 0 2 ) ( 0 2 ) ( 0 1 0 ) ( 0 1 0 ) ( 1 0 ) -- > luDecomp' ( 2 1 ) = ( ( 0 0 ) , ( 1/2 -1/4 1 ) , ( 1 0 0 ) , ( 0 1 ) , -1 , 1 )-luDecomp' :: (Ord a, Fractional a) => Matrix a -> (Matrix a,Matrix a,Matrix a,Matrix a,a,a)+luDecomp' :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a,Matrix a,Matrix a,Matrix a,a,a) luDecomp' a = recLUDecomp' a i i (identity $ ncols a) 1 1 1 n where i = identity $ nrows a n = min (nrows a) (ncols a) +-- | Unsafe version of 'luDecomp''. It fails when the input matrix is singular.+luDecompUnsafe' :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, Matrix a, a, a)+luDecompUnsafe' m = case luDecomp' m of+ Just x -> x+ _ -> error "luDecompUnsafe' of singular matrix."+ recLUDecomp' :: (Ord a, Fractional a) => Matrix a -- ^ U -> Matrix a -- ^ L@@ -799,11 +843,13 @@ -> a -- ^ e -> Int -- ^ Current row -> Int -- ^ Total rows- -> (Matrix a,Matrix a,Matrix a,Matrix a,a,a)+ -> Maybe (Matrix a,Matrix a,Matrix a,Matrix a,a,a) recLUDecomp' u l p q d e k n = if k > n || u'' ! (k, k) == 0- then (u,l,p,q,d,e)- else recLUDecomp' u'' l'' p' q' d' e' (k+1) n+ then Just (u,l,p,q,d,e)+ else if ukk == 0+ then Nothing+ else recLUDecomp' u'' l'' p' q' d' e' (k+1) n where -- Pivot strategy: maximum value in absolute value below the current row & col. (i, j) = maximumBy (comparing (\(i0, j0) -> abs $ u ! (i0,j0)))@@ -892,23 +938,27 @@ -- DETERMINANT {-# RULES-"matrix/detOfProduct"+"matrix/detLaplaceProduct" forall a b. detLaplace (a*b) = detLaplace a * detLaplace b -"matrix/detLUOfProduct"+"matrix/detLUProduct" forall a b. detLU (a*b) = detLU a * detLU b #-} -- | Matrix determinant using Laplace expansion. -- If the elements of the 'Matrix' are instance of 'Ord' and 'Fractional' -- consider to use 'detLU' in order to obtain better performance.+-- Function 'detLaplace' is /extremely/ slow. detLaplace :: Num a => Matrix a -> a detLaplace (M 1 1 v) = V.head v-detLaplace m =- sum [ (-1)^(i-1) * m ! (i,1) * detLaplace (minorMatrix i 1 m) | i <- [1 .. nrows m] ]+detLaplace m = sum1 [ (-1)^(i-1) * m ! (i,1) * detLaplace (minorMatrix i 1 m) | i <- [1 .. nrows m] ]+ where+ sum1 = foldl1' (+) -- | Matrix determinant using LU decomposition.+-- It works even when the input matrix is singular. detLU :: (Ord a, Fractional a) => Matrix a -> a-detLU m = d * diagProd u- where- (u,_,_,d) = luDecomp m+detLU m = case luDecomp m of+ Just (u,_,_,d) -> d * diagProd u+ Nothing -> 0+
license view
@@ -1,4 +1,4 @@-Copyright (c)2013, Daniel Díaz +Copyright (c)2014, Daniel Díaz All rights reserved.
matrix.cabal view
@@ -1,5 +1,5 @@ Name: matrix -Version: 0.2.4.0 +Version: 0.3.0.0 Author: Daniel Díaz Category: Math Build-type: Simple @@ -29,7 +29,7 @@ location: git://github.com/Daniel-Diaz/matrix.git Library - Build-depends: base ==4.* + Build-depends: base == 4.* , vector >= 0.10 && < 0.11 , deepseq >= 1.3.0.0 && < 1.4 , primitive >= 0.5 && < 0.6 @@ -40,7 +40,17 @@ type: exitcode-stdio-1.0 hs-source-dirs: bench main-is: mult.hs - build-depends: base ==4.* + build-depends: base == 4.* , matrix , criterion ghc-options: -O2 + +Test-Suite matrix-test + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: Main.hs + build-depends: base == 4.* + , matrix + , tasty + , QuickCheck + , tasty-quickcheck
+ test/Main.hs view
@@ -0,0 +1,102 @@++import Data.Matrix+import Data.Ratio+import Control.Applicative+import Data.Monoid (mconcat)++import Test.Tasty+import qualified Test.Tasty.QuickCheck as QC+import Test.QuickCheck++{- matrix package test set++This program uses QuickCheck to check that the matrix+functions of the matrix package are working properly.++We use the type Rational to have avoid numerical errors+that may cause the test to fail while the algorithm is+correct.++-}++-- | Numbers used in tests.+type R = Rational++newtype I = I { fromI :: Int }++instance Show I where+ show (I n) = show n++instance Arbitrary I where+ arbitrary = I <$> choose (1,9)++instance Arbitrary a => Arbitrary (Matrix a) where+ arbitrary = do+ I n <- arbitrary+ I m <- arbitrary+ genMatrix' n m++genMatrix' :: Arbitrary a => Int -> Int -> Gen (Matrix a)+genMatrix' n m = fromList n m <$> vector (n*m)++genMatrix :: Int -> Int -> Gen (Matrix R)+genMatrix = genMatrix'+++-- | Square matrices+newtype Sq = Sq { fromSq :: Matrix R }++instance Show Sq where+ show (Sq m) = show m++instance Arbitrary Sq where+ arbitrary = do+ I n <- arbitrary+ Sq <$> genMatrix n n++main :: IO ()+main = defaultMain $ testGroup "matrix tests" [+ QC.testProperty "identity * m = m * identity = m"+ $ \(Sq m) -> let n = nrows m in identity n * m == m && m * identity n == m+ , QC.testProperty "getMatrixAsVector m = mconcat [ getRow i m | i <- [1 .. nrows m]]"+ $ \m -> getMatrixAsVector (m :: Matrix R) == mconcat [ getRow i m | i <- [1 .. nrows m] ]+ , QC.testProperty "permMatrix n i j * permMatrix n i j = identity n"+ $ \(I n) -> forAll (choose (1,n))+ $ \i -> forAll (choose (1,n))+ $ \j -> permMatrix n i j * permMatrix n i j == identity n+ , QC.testProperty "setElem (getElem i j m) (i,j) m = m"+ $ \m -> forAll (choose (1,nrows m))+ $ \i -> forAll (choose (1,ncols m))+ $ \j -> setElem (getElem i j m) (i,j) m == (m :: Matrix R)+ , QC.testProperty "transpose (transpose m) = m"+ $ \m -> transpose (transpose m) == (m :: Matrix R)+ , QC.testProperty "getRow i m = getCol i (transpose m)"+ $ \m -> forAll (choose (1,nrows m))+ $ \i -> getRow i (m :: Matrix R) == getCol i (transpose m)+ , QC.testProperty "joinBlocks (splitBlocks i j m) = m"+ $ \m -> forAll (choose (1,nrows m))+ $ \i -> forAll (choose (1,ncols m))+ $ \j -> joinBlocks (splitBlocks i j m) == (m :: Matrix R)+ , QC.testProperty "(+) = elementwise (+)"+ $ \m1 -> forAll (genMatrix (nrows m1) (ncols m1))+ $ \m2 -> m1 + m2 == elementwise (+) m1 m2+ , QC.testProperty "if (u,l,p,d) = luDecomp m then (p*m = l*u) && (detLaplace p = d)"+ $ \(Sq m) -> (detLaplace m /= 0) ==>+ (let (u,l,p,d) = luDecompUnsafe m in p*m == l*u && detLaplace p == d)+ , QC.testProperty "detLaplace m = detLU m"+ $ \(Sq m) -> detLaplace m == detLU m+ , QC.testProperty "if (u,l,p,q,d,e) = luDecomp' m then (p*m*q = l*u) && (detLU p = d) && (detLU q = e)"+ $ \(Sq m) -> (detLU m /= 0) ==>+ (let (u,l,p,q,d,e) = luDecompUnsafe' m in p*m*q == l*u && detLU p == d && detLU q == e)+ , QC.testProperty "detLU (scaleRow k i m) = k * detLU m"+ $ \(Sq m) k -> forAll (choose (1,nrows m))+ $ \i -> detLU (scaleRow k i m) == k * detLU m+ , QC.testProperty "let n = nrows m in detLU (switchRows i j m) = detLU (permMatrix n i j) * detLU m"+ $ \(Sq m) -> let n = nrows m in forAll (choose (1,n))+ $ \i -> forAll (choose (1,n))+ $ \j -> detLU (switchRows i j m) == detLU (permMatrix n i j) * detLU m+ , QC.testProperty "switchCols i j = transpose . switchRows i j . transpose"+ $ \m -> forAll (choose (1,ncols m))+ $ \i -> forAll (choose (1,ncols m))+ $ \j -> switchCols i j (m :: Matrix R) == (transpose $ switchRows i j $ transpose m)+ ]