diff --git a/Data/Matrix.hs b/Data/Matrix.hs --- a/Data/Matrix.hs +++ b/Data/Matrix.hs @@ -1,4 +1,4 @@ - +{-# LANGUAGE DeriveGeneric #-} -- | Matrix datatype and operations. -- -- Every provided example has been tested. @@ -23,7 +23,7 @@ , toList , toLists -- * Accessing , getElem , (!) , unsafeGet , safeGet, safeSet - , getRow , getCol + , getRow , safeGetRow , getCol , safeGetCol , getDiag , getMatrixAsVector -- * Manipulating matrices @@ -31,7 +31,7 @@ , unsafeSet , transpose , setSize , extendTo , inverse, rref - , mapRow , mapCol + , mapRow , mapCol, mapPos -- * Submatrices -- ** Splitting blocks , submatrix @@ -77,15 +77,16 @@ import Data.Foldable (Foldable, foldMap, foldl1) import Data.Maybe import Data.Monoid +import qualified Data.Semigroup as S import Data.Traversable import Control.Applicative(Applicative, (<$>), (<*>), pure) +import GHC.Generics (Generic) -- Data import Control.Monad.Primitive (PrimMonad, PrimState) import Data.List (maximumBy,foldl1') import Data.Ord (comparing) import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV -import Data.Maybe ------------------------------------------------------- ------------------------------------------------------- @@ -114,7 +115,7 @@ , colOffset :: {-# UNPACK #-} !Int , vcols :: {-# UNPACK #-} !Int -- ^ Number of columns of the matrix without offset , mvect :: V.Vector a -- ^ Content of the matrix as a plain vector. - } + } deriving (Generic) instance Eq a => Eq (Matrix a) where m1 == m2 = @@ -129,12 +130,19 @@ -- | Display a matrix as a 'String' using the 'Show' instance of its elements. prettyMatrix :: Show a => Matrix a -> String -prettyMatrix m@(M _ _ _ _ _ v) = unlines - [ "( " <> unwords (fmap (\j -> fill mx $ show $ m ! (i,j)) [1..ncols m]) <> " )" | i <- [1..nrows m] ] +prettyMatrix m = concat + [ "┌ ", unwords (replicate (ncols m) blank), " ┐\n" + , unlines + [ "│ " ++ unwords (fmap (\j -> fill $ strings ! (i,j)) [1..ncols m]) ++ " │" | i <- [1..nrows m] ] + , "└ ", unwords (replicate (ncols m) blank), " ┘" + ] where - mx = V.maximum $ fmap (length . show) v - fill k str = replicate (k - length str) ' ' ++ str + strings@(M _ _ _ _ _ v) = fmap show m + widest = V.maximum $ fmap length v + fill str = replicate (widest - length str) ' ' ++ str + blank = fill "" + instance Show a => Show (Matrix a) where show = prettyMatrix @@ -164,8 +172,11 @@ ------------------------------------------------------- ---- MONOID INSTANCE +instance Monoid a => S.Semigroup (Matrix a) where + (<>) = mappend + instance Monoid a => Monoid (Matrix a) where - mempty = fromList 1 1 [mempty] + mempty = fromList 1 1 [mempty] mappend m m' = matrix (max (nrows m) (nrows m')) (max (ncols m) (ncols m')) $ uncurry zipTogether where zipTogether row column = fromMaybe mempty $ safeGet row column m <> safeGet row column m' @@ -178,11 +189,11 @@ ------------------------------------------------------- ------------------------------------------------------- ---- APPLICATIVE INSTANCE ----- Works like tensor product but applies a function +---- Works like tensor product but applies a function instance Applicative Matrix where - pure x = fromList 1 1 [x] - m <*> m' = flatten $ ((\f -> f <$> m') <$> m) + pure x = fromList 1 1 [x] + m <*> m' = flatten $ (\f -> f <$> m') <$> m ------------------------------------------------------- @@ -191,8 +202,8 @@ -- | Flatten a matrix of matrices. All sub matrices must have same dimensions --- This criteria is not checked. -flatten:: (Matrix (Matrix a)) -> Matrix a +-- This criteria is not checked. +flatten:: Matrix (Matrix a) -> Matrix a flatten m = foldl1 (<->) $ map (foldl1 (<|>) . (\i -> getRow i m)) [1..(nrows m)] -- | /O(rows*cols)/. Map a function over a row. @@ -229,6 +240,20 @@ then f i a else a + +-- | /O(rows*cols)/. Map a function over elements. +-- Example: +-- +-- > ( 1 2 3 ) ( 0 -1 -2 ) +-- > ( 4 5 6 ) ( 1 0 -1 ) +-- > mapPos (\(r,c) a -> r - c) ( 7 8 9 ) = ( 2 1 0 ) +-- +mapPos :: ((Int, Int) -> a -> b) -- ^ Function takes the current Position as additional argument. + -> Matrix a + -> Matrix b +mapPos f m@(M {ncols = cols, mvect = vect})= + m { mvect = V.imap (\i e -> f (decode cols i) e) vect} + ------------------------------------------------------- ------------------------------------------------------- ---- FOLDABLE AND TRAVERSABLE INSTANCES @@ -422,14 +447,15 @@ -> a {-# INLINE getElem #-} getElem i j m = - case safeGet i j m of - Just x -> x - Nothing -> error - $ "getElem: Trying to get the " - ++ show (i,j) - ++ " element from a " - ++ sizeStr (nrows m) (ncols m) - ++ " matrix." + fromMaybe + (error $ + "getElem: Trying to get the " + ++ show (i, j) + ++ " element from a " + ++ sizeStr (nrows m) (ncols m) + ++ " matrix." + ) + (safeGet i j m) -- | /O(1)/. Unsafe variant of 'getElem', without bounds checking. unsafeGet :: Int -- ^ Row @@ -466,11 +492,23 @@ {-# INLINE getRow #-} getRow i (M _ m ro co w v) = V.slice (w*(i-1+ro) + co) m v +-- | Varian of 'getRow' that returns a maybe instead of an error +safeGetRow :: Int -> Matrix a -> Maybe (V.Vector a) +safeGetRow r m + | r > nrows m || r < 1 = Nothing + | otherwise = Just $ getRow r m + -- | /O(rows)/. Get a column of a matrix as a vector. getCol :: Int -> Matrix a -> V.Vector a {-# INLINE getCol #-} getCol j (M n _ ro co w v) = V.generate n $ \i -> v V.! encode w (i+1+ro,j+co) +-- | Varian of 'getColumn' that returns a maybe instead of an error +safeGetCol :: Int -> Matrix a -> Maybe (V.Vector a) +safeGetCol c m + | c > ncols m || c < 1 = Nothing + | otherwise = Just $ getCol c m + -- | /O(min rows cols)/. Diagonal of a /not necessarily square/ matrix. getDiag :: Matrix a -> V.Vector a getDiag m = V.generate k $ \i -> m ! (i+1,i+1) @@ -534,7 +572,7 @@ transpose :: Matrix a -> Matrix a transpose m = matrix (ncols m) (nrows m) $ \(i,j) -> m ! (j,i) --- | /O(rows*rows*rows) = O(cols*cols*cols)/. The inverse of a square matrix. +-- | /O(rows*rows*rows*rows) = O(cols*cols*cols*cols)/. The inverse of a square matrix. -- Uses naive Gaussian elimination formula. inverse :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a) inverse m @@ -548,7 +586,7 @@ rref'd = rref adjoinedWId in rref'd >>= return . submatrix 1 (nrows m) (ncols m + 1) (ncols m * 2) --- | /O(rows*rows*cols)/. Converts a matrix to reduced row echelon form, thus +-- | /O(rows*rows*cols*cols)/. Converts a matrix to reduced row echelon form, thus -- solving a linear system of equations. This requires that (cols > rows) -- if cols < rows, then there are fewer variables than equations and the -- problem cannot be solved consistently. If rows = cols, then it is @@ -568,6 +606,7 @@ | nrows mtx == 1 = Right mtx | otherwise = let + -- this is super-slow: [resolvedRight] is cubic because [combineRows] is quadratic resolvedRight = foldr (.) id (map resolveRow [1..col-1]) mtx where col = nrows mtx @@ -591,10 +630,11 @@ sigAtTop = switchRows 1 goodRow mtx where significantRow n = getElem n 1 mtx /= 0 - goodRow = case listToMaybe (filter significantRow [1..ncols mtx]) of - Nothing -> error "Attempt to invert a non-invertible matrix" - Just x -> x - normalizedFirstRow = scaleRow (1 / getElem 1 1 mtx) 1 sigAtTop + goodRow = fromMaybe + (error "Attempt to invert a non-invertible matrix") + (listToMaybe (filter significantRow [1 .. ncols mtx])) + + normalizedFirstRow = scaleRow (1 / getElem 1 1 sigAtTop) 1 sigAtTop clearedLeft = foldr (.) id (map combinator [2..nrows mtx]) normalizedFirstRow where combinator n = combineRows n (-getElem n 1 normalizedFirstRow) 1 @@ -865,7 +905,7 @@ -- B b11 = b !. (1,1) ; b12 = b !. (1,2) b21 = b !. (2,1) ; b22 = b !. (2,2) - in V.fromList + in V.fromList [ a11*b11 + a12*b21 , a11*b12 + a12*b22 , a21*b11 + a22*b21 , a21*b12 + a22*b22 ] @@ -1191,7 +1231,7 @@ in if i == k then l else M (nrows l) (ncols l) lro lco lw $ - V.modify (\mv -> forM_ [1 .. k-1] $ + V.modify (\mv -> forM_ [1 .. k-1] $ \j -> MV.swap mv (en (i+lro,j+lco)) (en (k+lro,j+lco)) ) $ mvect l diff --git a/matrix.cabal b/matrix.cabal --- a/matrix.cabal +++ b/matrix.cabal @@ -1,5 +1,5 @@ Name: matrix -Version: 0.3.5.0 +Version: 0.3.6.0 Author: Daniel Díaz Category: Math Build-type: Simple @@ -33,6 +33,7 @@ , vector >= 0.10 , deepseq >= 1.3.0.0 && < 1.5 , primitive >= 0.5 + , semigroups >= 0.9 , loop >= 0.2 Exposed-modules: Data.Matrix GHC-Options: -Wall -O2 -fstatic-argument-transformation diff --git a/test/Examples.hs b/test/Examples.hs --- a/test/Examples.hs +++ b/test/Examples.hs @@ -2,6 +2,7 @@ import Data.Matrix import System.Exit (exitFailure) import System.IO (hFlush,stdout) +import Data.Either (isLeft) -- We flush stdout explictly to get output printed -- in real-time in Windows systems. @@ -75,6 +76,10 @@ ( mapCol (\_ x -> x + 1) 2 $ fromList 3 3 [1..9] , fromList 3 3 [1,3,3 , 4,6,6 , 7,9,9] ) + , testEquality "mapPos" + ( mapPos (\(r,c) x -> r + 2 * c) $ fromList 2 2 [1..4] + , fromList 2 2 [3, 5, 4, 6] + ) , testEquality "submatrix" ( submatrix 1 2 2 3 $ fromList 3 3 [1..9] , fromList 2 2 [2,3 , 5,6] @@ -111,20 +116,27 @@ ( toLists $ fromList 3 3 [1..9] , [ [1,2,3] , [4,5,6] , [7,8,9] ] ) - , testEquality "inverse" + , testEquality "inverse (1)" ( inverse $ fromList 2 2 [1,7, 2,4] , Right $ fromList 2 2 [-4/10,7/10, 2/10,-1/10] :: Either String (Matrix Rational) ) - , testEquality "inverse (1)" + , testEquality "inverse (2)" ( inverse $ fromList 3 3 [1,7,-12, 2,4,10, 0,-23,1] , Right $ fromList 3 3 [117/386, 269/772, 59/386, -1/386, 1/772, -17/386, -23/386, 23/772, -5/386] :: Either String (Matrix Rational)) - , testEquality "inverse (2)" + , testEquality "inverse (3)" ( inverse $ fromList 4 4 [1,2345,23,78, 12,34556,123,-1242, 429,-131,0,0, 0,0,0,-1] , Right $ fromList 4 4 [ -5371/72415160, 3013/217245480, 506353/217245480, -41658/1810379, -17589/72415160, 3289/72415160, -51/72415160, -136422/1810379, 617754/9051895, -125767/27155685, -802/27155685, 20050470/1810379, 0, 0, 0, -1] :: Either String (Matrix Rational)) + , testEquality "inverse (4)" + ( inverse $ fromList 2 2 [0,1, 1,0] + , Right $ fromList 2 2 [0,1, 1,0] :: Either String (Matrix Rational)) + , testEquality "inverse (5)" + ( inverse $ fromList 3 3 [1,0,0, 0,0,1, 0,1,0] + , Right $ fromList 3 3 [1,0,0, 0,0,1, 0,1,0] :: Either String (Matrix Rational)) + , testExample "inverse (6)" $ isLeft $ inverse $ fromList 2 2 [1,1,2,2] ]