diff --git a/Data/Matrix.hs b/Data/Matrix.hs --- a/Data/Matrix.hs +++ b/Data/Matrix.hs @@ -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 - diff --git a/matrix.cabal b/matrix.cabal --- a/matrix.cabal +++ b/matrix.cabal @@ -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 diff --git a/test/Examples.hs b/test/Examples.hs --- a/test/Examples.hs +++ b/test/Examples.hs @@ -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)) ] diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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)) ]