matrix 0.3.4.4 → 0.3.5.0

raw patch · 4 files changed

+166/−4 lines, 4 filesdep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• Data/Matrix.hs +140/−3
• matrix.cabal +1/−1
• test/Examples.hs +19/−0
• test/Main.hs +6/−0

@@ -15,12 +15,14 @@     -- ** Special matrices   , zero   , identity+  , diagonalList+  , diagonal   , permMatrix     -- * List conversions   , fromList , fromLists   , toList   , toLists     -- * Accessing-  , getElem , (!) , unsafeGet , safeGet+  , getElem , (!) , unsafeGet , safeGet, safeSet   , getRow  , getCol   , getDiag   , getMatrixAsVector@@ -28,6 +30,7 @@   , setElem   , unsafeSet   , transpose , setSize , extendTo+  , inverse, rref   , mapRow , mapCol     -- * Submatrices     -- ** Splitting blocks@@ -63,21 +66,26 @@     -- ** Determinants   , detLaplace   , detLU+  , flatten   ) where +import Prelude hiding (foldl1) -- Classes import Control.DeepSeq import Control.Monad (forM_) import Control.Loop (numLoop,numLoopFold)-import Data.Foldable (Foldable, foldMap)+import Data.Foldable (Foldable, foldMap, foldl1)+import Data.Maybe import Data.Monoid import Data.Traversable+import Control.Applicative(Applicative, (<$>), (<*>), pure) -- 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  ------------------------------------------------------- -------------------------------------------------------@@ -152,6 +160,41 @@ ------------------------------------------------------- ------------------------------------------------------- +-------------------------------------------------------+-------------------------------------------------------+---- MONOID INSTANCE++instance Monoid a => Monoid (Matrix a) where+  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'+++-------------------------------------------------------+-------------------------------------------------------+-------------------------------------------------------+-------------------------------------------------------++-------------------------------------------------------+-------------------------------------------------------+---- APPLICATIVE INSTANCE+---- 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)+++-------------------------------------------------------+-------------------------------------------------------++++-- | Flatten a matrix of matrices. All sub matrices must have same dimensions+--   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. --   Example: --@@ -249,6 +292,15 @@ identity :: Num a => Int -> Matrix a identity n = matrix n n $ \(i,j) -> if i == j then 1 else 0 +-- | Similar to 'diagonalList', but using 'V.Vector', which+--   should be more efficient.+diagonal :: a -- ^ Default element+         -> V.Vector a  -- ^ Diagonal vector+         -> Matrix a+diagonal e v = matrix n n $ \(i,j) -> if i == j then V.unsafeIndex v (i - 1) else e+  where+    n = V.length v+ -- | Create a matrix from a non-empty list given the desired size. --   The list must have at least /rows*cols/ elements. --   An example:@@ -283,6 +335,21 @@ toLists :: Matrix a -> [[a]] toLists m = [ [ unsafeGet i j m | j <- [1 .. ncols m] ] | i <- [1 .. nrows m] ] +-- | Diagonal matrix from a non-empty list given the desired size.+--   Non-diagonal elements will be filled with the given default element.+--   The list must have at least /order/ elements.+--+-- > diagonalList n 0 [1..] =+-- >                   n+-- >   1 ( 1 0 ... 0   0 )+-- >   2 ( 0 2 ... 0   0 )+-- >     (     ...       )+-- >     ( 0 0 ... n-1 0 )+-- >   n ( 0 0 ... 0   n )+--+diagonalList :: Int -> a -> [a] -> Matrix a+diagonalList n e xs = matrix n n $ \(i,j) -> if i == j then xs !! (i - 1) else e+ -- | Create a matrix from a non-empty list of non-empty lists. --   /Each list must have at least as many elements as the first list/. --   Examples:@@ -388,6 +455,12 @@  | i > n || j > m || i < 1 || j < 1 = Nothing  | otherwise = Just $ unsafeGet i j a +-- | Variant of 'setElem' that returns Maybe instead of an error.+safeSet:: a -> (Int, Int) -> Matrix a -> Maybe (Matrix a)+safeSet x p@(i,j) a@(M n m _ _ _ _)+  | i > n || j > m || i < 1 || j < 1 = Nothing+  | otherwise = Just $ unsafeSet x p a+ -- | /O(1)/. Get a row of a matrix as a vector. getRow :: Int -> Matrix a -> V.Vector a {-# INLINE getRow #-}@@ -461,6 +534,71 @@ 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.+--   Uses naive Gaussian elimination formula.+inverse :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a)+inverse m+    | ncols m /= nrows m+        = Left+            $ "Inverting non-square matrix with dimensions "+                ++ show (sizeStr (ncols m) (nrows m))+    | otherwise =+        let+            adjoinedWId = m <|> identity (nrows m)+            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+--  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+--  basically a homogenous system of equations, so it will be reduced to+--  identity or an error depending on whether the marix is invertible+--  (this case is allowed for robustness).+rref :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a)+rref m+        | ncols m < nrows m+            = Left $+                "Invalid dimensions "+                    ++ show (sizeStr (ncols m) (nrows m))+                    ++ "; the number of columns must be greater than or equal to the number of rows"+        | otherwise             = rrefRefd (ref m)+    where+    rrefRefd mtx+        | nrows mtx == 1    = Right mtx+        | otherwise =+            let+                resolvedRight = foldr (.) id (map resolveRow [1..col-1]) mtx+                    where+                    col = nrows mtx+                    resolveRow n = combineRows n (-getElem n col mtx) col+                top = submatrix 1 (nrows resolvedRight - 1) 1 (ncols resolvedRight) resolvedRight+                top' = rrefRefd top+                bot = submatrix (nrows resolvedRight) (nrows resolvedRight) 1 (ncols resolvedRight) resolvedRight+            in top' >>= return . (<-> bot)+++ref :: (Fractional a, Eq a) => Matrix a -> Matrix a+ref mtx+        | nrows mtx == 1+            = clearedLeft+        | otherwise =+            let+                (tl, tr, bl, br) = splitBlocks 1 1 clearedLeft+                br' = ref br+            in (tl <|> tr) <-> (bl <|> br')+    where+    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+    clearedLeft = foldr (.) id (map combinator [2..nrows mtx]) normalizedFirstRow+        where+        combinator n = combineRows n (-getElem n 1 normalizedFirstRow) 1+ -- | 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.@@ -1233,4 +1371,3 @@ detLU m = case luDecomp m of   Just (u,_,_,d) -> d * diagProd u   Nothing -> 0-

@@ -1,5 +1,5 @@ Name: matrix
-Version: 0.3.4.4
+Version: 0.3.5.0
 Author: Daniel Díaz
 Category: Math
 Build-type: Simple

@@ -51,6 +51,9 @@   , testEquality "identity"       ( identity 3 , fromList 3 3 [1,0,0 , 0,1,0 , 0,0,1]         )+  , testEquality "diagonalList"+      ( diagonalList 3 0 [1..] , fromList 3 3 [1,0,0 , 0,2,0 , 0,0,3]+        )   , testEquality "transpose"       ( transpose $ fromList 3 3 [1..9]       , fromList 3 3 [1,4,7 , 2,5,8 , 3,6,9]@@ -108,4 +111,20 @@       ( toLists $ fromList 3 3 [1..9]       , [ [1,2,3] , [4,5,6] , [7,8,9] ]         )+  , testEquality "inverse"+      ( 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)"+      ( 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)"+      ( 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))     ]

@@ -125,4 +125,10 @@        $ \m -> fromList (nrows m) (ncols m) (toList m) == (m :: Matrix R)   , QC.testProperty "fromLists . toLists = id"        $ \m -> fromLists (toLists m) == (m :: Matrix R)+  , QC.testProperty "inv m * m = Right identity"+       $ \sq -> let m = fromSq sq in detLU m /= 0 ==> (multStd m <$> inverse m) == Right (identity (nrows m) :: Matrix R)+  , QC.testProperty "inv . inv == id"+       $ \sq -> let m = fromSq sq in detLU m /= 0 ==> (inverse m >>= inverse) == Right (m :: Matrix R)+  , QC.testProperty "rref . fromSquare = const (identity)"+       $ \sq -> let m = fromSq sq :: Matrix R in (detLU m /= 0) ==> rref m == Right (identity (ncols m))     ]