diff --git a/bed-and-breakfast.cabal b/bed-and-breakfast.cabal --- a/bed-and-breakfast.cabal +++ b/bed-and-breakfast.cabal @@ -1,5 +1,5 @@ Name: bed-and-breakfast -Version: 0.1.4 +Version: 0.2 Synopsis: Efficient Matrix operations in 100% Haskell. Description: Efficient Matrix operations in 100% Haskell. . @@ -18,10 +18,16 @@ invertible. Also @Matrix@ derives 'Data.Typeable' now. . - [@v0.1.4@] Added @scale@, @<->@, and @<|>@. Corrected - a bug in @isUnit@ reported by Charles Durham - (would have returned True for any matrix for which + [@v0.1.4@] Added @scale@, and methods for joining + matrices vertically and horizontally. Corrected + a bug in @isUnit@ reported by Charles Durham. + @isUnit@ returned True for any matrix for which @all (== 1) . trace@ would have, which is wrong). + . + [@v0.2@] A little bit more documentation. Also moved some + functions (@isXXX@) away from the type class @MatrixElement@. + Properly flagged the package as experimental (was + improperly copied as @stable@, copied form a template). License: MIT License-File: LICENSE @@ -30,7 +36,7 @@ Build-Type: Simple Cabal-Version: >= 1.8 Category: Numeric, Math, Linear Algebra -Stability: stable +Stability: experimental Homepage: http://hub.darcs.net/scravy/bed-and-breakfast Source-Repository head diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -3,57 +3,98 @@ , FlexibleContexts , Trustworthy #-} ---{-# OPTIONS -Wall -fno-warn-name-shadowing #-} +{-# OPTIONS -Wall -fno-warn-name-shadowing #-} +-- | Efficient matrix operations in 100% pure Haskell. +-- +-- This package uses miscellaneous implementations, +-- depending on the type of its components. Typically unboxed +-- arrays will perform best, while unboxed arrays give you +-- certain features such as 'Rational' or 'Complex' components. +-- +-- The following component types are supported by 'Matrix': +-- +-- [@Int@] Uses unboxed arrays internally. 'inv' will always +-- return 'Nothing'. +-- +-- [@Integer@] Uses boxed arrays internally. 'inv' will always +-- return 'Nothing'. +-- +-- [@Double@ and @Float@] Uses unboxed arrays internally. +-- All matrix operations will work as expected. +-- @Matrix Double@ will probably yield the best peformance. +-- +-- [@Rational@] Best choice if precision is what you aim for. +-- Uses boxed arrays internally. All matrix operations will +-- work as expected. +-- +-- [@Complex@] Experimental. Uses boxed arrays internally. +-- The current implementation of 'inv' requires an instance +-- of 'Ord' for the component type, therefor it is currently +-- not possible to calculate the inverse of a complex matrix +-- (on my to do list). module Numeric.Matrix ( + Matrix, - MatrixElement ( - matrix, - fromList, - unit, - zero, - diag, - empty, + MatrixElement (..), + + -- * Matrix property and utility functions. - at, - row, - col, + -- | Joins two matrices horizontally. + -- + -- > 1 2 3 1 0 0 1 2 3 1 0 0 + -- > 3 4 5 <|> 2 1 0 -> 3 4 5 2 1 0 + -- > 5 6 7 3 2 1 5 6 7 3 2 1 + (<|>), - select, - toList, + -- | Joins two matrices vertically. + -- + -- > 1 2 3 1 0 0 1 2 3 + -- > 3 4 5 <-> 2 1 0 -> 3 4 5 + -- > 5 6 7 3 2 1 5 6 7 + -- > 1 0 0 + -- > 2 1 0 + -- > 3 2 1 + (<->), - dimensions, - numRows, - numCols, + -- | Scales a matrix by the given factor. + -- + -- > scale s == map (*s) + scale, - isUnit, - isZero, - isDiagonal, - isEmpty, - isSquare, - - det, - rank, - transpose, - trace, + -- * Matrix properties - minus, - plus, - times, - inv, + -- | Check whether the matrix is an identity matrix. + -- + -- > 1 0 0 + -- > 0 1 0 + -- > 0 0 1 (True) + isUnit, - map, - all, - any, + -- | Check whether the matrix consists of all zeros. + -- + -- > isZero == all (== 0) + isZero, - mapWithIndex, - allWithIndex, - anyWithIndex - ), - (<|>), - (<->), - scale + -- | Checks whether the matrix is a diagonal matrix. + -- + -- > 4 0 0 0 + -- > 0 7 0 0 + -- > 0 0 3 0 + -- > 0 0 0 9 (True) + isDiagonal, + + -- | Checks whether the matrix is empty. + -- + -- > isEmpty m = numCols == 0 || numRows == 0 + isEmpty, + + -- | Checks whether the matrix is a square matrix. + -- + -- > isSquare == uncurry (==) . dimensions + isSquare + ) where @@ -66,8 +107,8 @@ import Data.Ratio import Data.Complex import Data.Maybe -import Data.Foldable (Foldable) -import qualified Data.Foldable as F +--import Data.Foldable (Foldable) +--import qualified Data.Foldable as F import qualified Data.List as L import Data.Array.IArray @@ -129,6 +170,7 @@ rnf matrix = matrix `deepseq` () +(<|>) :: MatrixElement e => Matrix e -> Matrix e -> Matrix e m1 <|> m2 = let m = numCols m1 n1 = numRows m1 n2 = numRows m2 @@ -137,6 +179,7 @@ then (if i > n2 then 0 else m2 `at` (i,j-m)) else (if i > n1 then 0 else m1 `at` (i,j)) +(<->) :: MatrixElement e => Matrix e -> Matrix e -> Matrix e m1 <-> m2 = let m = numRows m1 n1 = numCols m1 n2 = numCols m2 @@ -149,6 +192,17 @@ scale m s = map (*s) m +isUnit, isDiagonal, isZero, isEmpty, isSquare :: MatrixElement e => Matrix e -> Bool + +isZero = all (== 0) +isUnit m = isSquare m && allWithIndex (uncurry check) m + where check = \i j e -> if i == j then e == 1 else e == 0 +isEmpty m = numRows m == 0 || numCols m == 0 +isDiagonal m = isSquare m && allWithIndex (uncurry check) m + where check = \i j e -> if i /= j then e == 0 else True +isSquare m = let (a, b) = dimensions m in a == b + + class Division e where divide :: e -> e -> e @@ -178,12 +232,22 @@ numRows :: Matrix e -> Int numCols :: Matrix e -> Int + + -- Builds a matrix from a list of lists. + -- + -- > fromList [[1,2,3],[2,1,3],[3,2,1]] :: Matrix Rational fromList :: [[e]] -> Matrix e toList :: Matrix e -> [[e]] + -- An identity square matrix of the given size. unit :: Int -> Matrix e + + -- A square matrix of the given size consisting of all zeros. zero :: Int -> Matrix e + diag :: [e] -> Matrix e + + -- > dimensions empty == (0, 0) empty :: Matrix e minus :: Matrix e -> Matrix e -> Matrix e @@ -193,25 +257,30 @@ -- adjugate :: Matrix e -> Matrix e -- cofactors :: Matrix e -> Matrix e ; cofactors = undefined + + -- Applies Bareiss multistep integer-preserving + -- algorithm for finding the determinant of a matrix. + -- Returns 0 if the matrix is not a square matrix. det :: Matrix e -> e + + -- Flips rows and columns. + -- + -- 1 8 9 1 2 3 + -- 2 1 8 --transpose-> 8 1 2 + -- 3 2 1 9 8 1 transpose :: Matrix e -> Matrix e rank :: Matrix e -> e trace :: Matrix e -> [e] - isUnit :: Matrix e -> Bool - isDiagonal :: Matrix e -> Bool - isZero :: Matrix e -> Bool - isEmpty :: Matrix e -> Bool - isSquare :: Matrix e -> Bool - - select p m = [ at m (i,j) | i <- [1..numRows m] - , j <- [1..numCols m] - , p (i,j) ] - - at m (i, j) = ((!! j) . (!! i) . toList) m - + -- Applies a function on every component in the matrix. map :: MatrixElement f => (e -> f) -> Matrix e -> Matrix f + + -- Applies a predicate on every component in the matrix + -- and returns True if all components satisfy it. all :: (e -> Bool) -> Matrix e -> Bool + + -- Applies a predicate on every component in the matrix + -- and returns True if one or more components satisfy it. any :: (e -> Bool) -> Matrix e -> Bool mapWithIndex :: MatrixElement f => ((Int, Int) -> e -> f) -> Matrix e -> Matrix f @@ -223,6 +292,12 @@ empty = fromList [] diag xs = matrix (n,n) (\(i,j) -> if i == j then xs !! (i-1) else 0) where n = length xs + + select p m = [ at m (i,j) | i <- [1..numRows m] + , j <- [1..numCols m] + , p (i,j) ] + + at m (i, j) = ((!! j) . (!! i) . toList) m row i m = ((!! (i-1)) . toList) m col i m = (row i . transpose) m @@ -237,14 +312,6 @@ trace = select (uncurry (==)) inv _ = Nothing - isZero = all (== 0) - isUnit m = isSquare m && allWithIndex (uncurry check) m - where check = \i j e -> if i == j then e == 1 else e == 0 - isEmpty m = numRows m == 0 || numCols m == 0 - isDiagonal m = isSquare m && allWithIndex (uncurry check) m - where check = \i j e -> if i /= j then e == 0 else True - isSquare m = let (a, b) = dimensions m in a == b - map f = mapWithIndex (const f) all f = allWithIndex (const f) any f = anyWithIndex (const f) @@ -348,7 +415,7 @@ row i (ComplexMatrix _ _ arr) = _row i arr col j (ComplexMatrix _ _ arr) = _col j arr toList (ComplexMatrix _ _ arr) = _toList arr - inv (ComplexMatrix m n arr) = Nothing + inv (ComplexMatrix _ _ _) = Nothing --if m /= n then Nothing else -- Just $ ComplexMatrix m n $ runST (_inv boxedST arr) det (ComplexMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr)