diff --git a/CHANGES.md b/CHANGES.md --- a/CHANGES.md +++ b/CHANGES.md @@ -100,9 +100,18 @@
v0.4.3
- Fixed a bug in @transpose@ that prevented it from + Fixed a bug in transpose that prevented it from working correctly with non-square matrices. Thanks to @owst@ from @hub.darcs.net@. +
+ +
v0.5
+
+ inv works now for complex matrices too. + Added conversion functions toDoubleMatrix, + toComplexMatrix, toRationalMatrix. + Changed signature of minor and minorMatrix + to a more natural format.
diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ -bed-and-breakfast +bed-and-breakfast [![Build Status](https://travis-ci.org/scravy/bed-and-breakfast.svg?branch=master)](https://travis-ci.org/scravy/bed-and-breakfast) ================= -Matrix operations in 100% pure Haskell. +Matrix operations in 100% pure awesome Haskell. Bed and Breakfast is a linear algebra library written in Haskell. It provides fast matrix operations like finding the determinant diff --git a/ROADMAP.md b/ROADMAP.md --- a/ROADMAP.md +++ b/ROADMAP.md @@ -22,4 +22,4 @@ * Vector product * Outer product - +* Make it compatible with GHC 7.4 - 7.10 diff --git a/bed-and-breakfast.cabal b/bed-and-breakfast.cabal --- a/bed-and-breakfast.cabal +++ b/bed-and-breakfast.cabal @@ -1,7 +1,7 @@ Name: bed-and-breakfast -Version: 0.4.3 -Synopsis: Efficient Matrix operations in 100% Haskell. -Description: Efficient Matrix operations in 100% Haskell. +Version: 0.5 +Synopsis: Efficient Matrix and Vector operations in 100% Haskell. +Description: Efficient Matrix and Vector operations in 100% Haskell. . This library uses boxed and unboxed arrays in the ST monad, in order to achieve efficiency. @@ -23,13 +23,17 @@ Library Exposed-Modules: Numeric.Matrix, - Numeric.Matrix.Sugar + Numeric.Matrix.Sugar, + Numeric.Vector Build-Depends: base >= 4.5 && < 5, deepseq >= 1.3, array >= 0.4, binary >= 0.5, - template-haskell >= 2.7 + template-haskell >= 2.7, + cpphs >= 1.18 Hs-Source-Dirs: src + GHC-Options: -Wall -pgmPcpphs -optP--cpp + Extensions: CPP Test-Suite quickcheck type: exitcode-stdio-1.0 diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -2,8 +2,11 @@ , TypeFamilies , FlexibleContexts , Trustworthy + , StandaloneDeriving + , DeriveDataTypeable + , ConstrainedClassMethods #-} -{-# OPTIONS -Wall -fno-warn-name-shadowing #-} +{-# OPTIONS_GHC -Wall -fno-warn-name-shadowing #-} -- | Efficient matrix operations in 100% pure Haskell. -- @@ -29,10 +32,9 @@ -- 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). +-- All matrix operations will work as expected, though +-- finding the inverse of a matrix isa tad less numerically +-- stable than with a @Double@ matrix. module Numeric.Matrix ( Matrix, @@ -49,7 +51,12 @@ isZero, isDiagonal, isEmpty, - isSquare + isSquare, + + -- ** Conversions + toDoubleMatrix, + toComplexMatrix, + toRationalMatrix ) where @@ -78,9 +85,18 @@ import Data.Typeable -import Prelude hiding (any, all, read, map) +import Prelude (Show, Read, Num, Fractional, Eq, Bool (..), Integer, Integral, + Float, Double, RealFloat, Ord, Real, + (*), (/), (+), (-), (^), (.), (>=), (==), (/=), ($), (>), (!!), + (&&), (||), + undefined, null, head, zip, abs, flip, length, compare, drop, + negate, not, filter, fromIntegral, fst, snd, foldl1, min, max, + error, fromInteger, signum, lines, words, show, unwords, + unlines, + otherwise, id, const, uncurry, quot, toRational, fromRational) import qualified Prelude as P +import Data.Monoid -- | Matrices are represented by a type which fits best the component type. -- For example a @Matrix Double@ is represented by unboxed arrays, @@ -97,7 +113,7 @@ -- -- [@Show (Matrix e)@] -- Note that a Show instance for the component type @e@ must exist. --- +-- -- [@Read (Matrix e)@] -- You can read a matrix like so: -- @@ -135,6 +151,12 @@ -- See @encode@ and @decode@. data family Matrix e +#if defined(__GLASGOW_HASKELL__) && (__GLASGOW_HASKELL__ >= 707) +deriving instance Typeable Matrix +#else +deriving instance Typeable1 Matrix +#endif + data instance Matrix Int = IntMatrix !Int !Int (Array Int (UArray Int Int)) @@ -153,15 +175,6 @@ data instance Matrix (Complex a) = ComplexMatrix !Int !Int (Array Int (Array Int (Complex a))) - -instance Typeable a => Typeable (Matrix a) where - typeOf x = mkTyConApp (mkTyCon3 "bed-and-breakfast" - "Numeric.Matrix" - "Matrix") [typeOf (unT x)] - where - unT :: Matrix a -> a - unT = undefined - instance (MatrixElement e, Show e) => Show (Matrix e) where show = unlines . P.map showRow . toList where @@ -283,6 +296,16 @@ isSquare m = let (a, b) = dimensions m in a == b +toDoubleMatrix :: (MatrixElement a, Integral a) => Matrix a -> Matrix Double +toDoubleMatrix = map fromIntegral + +toRationalMatrix :: (MatrixElement a, Real a) => Matrix a -> Matrix Rational +toRationalMatrix = map toRational + +toComplexMatrix :: (MatrixElement a, RealFloat a, Show a) => Matrix a -> Matrix (Complex a) +toComplexMatrix = map (:+ 0) + + class Division e where divide :: e -> e -> e @@ -344,6 +367,12 @@ -- -- > fromList [[1,2,3],[2,1,3],[3,2,1]] :: Matrix Rational fromList :: [[e]] -> Matrix e + + -- | Turns a matrix into a list of lists. + -- + -- > (toList . fromList) xs == xs + -- + -- > (fromList . toList) mat == mat toList :: Matrix e -> [[e]] -- | An identity square matrix of the given size. @@ -410,12 +439,31 @@ -- | Compute the rank of a matrix. rank :: Matrix e -> e + + -- | Select the diagonal elements of a matrix as a list. + -- + -- > 1 8 3 + -- > 3 6 5 --trace-> [1, 6, 2] + -- > 7 4 2 trace :: Matrix e -> [e] - minor :: MatrixElement e => Matrix e -> (Int, Int) -> e + -- | Select the minor of a matrix, that is the determinant + -- of the 'minorMatrix'. + -- + -- > minor = det . minorMatrix + minor :: MatrixElement e => (Int, Int) -> Matrix e -> e + + -- | Select the minor matrix of a matrix, a matrix that is obtained + -- by deleting the i-th row and j-th column. + -- + -- > 10 9 95 45 + -- > 8 7 3 27 8 3 27 + -- > 13 17 19 23 --minorMatrix (1,2)-> 13 19 23 + -- > 1 2 5 8 1 5 8 + minorMatrix :: MatrixElement e => (Int, Int) -> Matrix e -> Matrix e + cofactors :: MatrixElement e => Matrix e -> Matrix e adjugate :: MatrixElement e => Matrix e -> Matrix e - minorMatrix :: MatrixElement e => Matrix e -> (Int, Int) -> Matrix e -- | Apply a function on every component in the matrix. map :: MatrixElement f => (e -> f) -> Matrix e -> Matrix f @@ -428,9 +476,16 @@ -- and return True if one or more components satisfy it. any :: (e -> Bool) -> Matrix e -> Bool + -- | Compute the sum of the components of the matrix. + sum :: Matrix e -> e + + -- | Map each component of the matrix to a monoid, and combine the results. + foldMap :: Monoid m => (e -> m) -> Matrix e -> m + mapWithIndex :: MatrixElement f => ((Int, Int) -> e -> f) -> Matrix e -> Matrix f allWithIndex :: ((Int, Int) -> e -> Bool) -> Matrix e -> Bool anyWithIndex :: ((Int, Int) -> e -> Bool) -> Matrix e -> Bool + foldMapWithIndex :: Monoid m => ((Int, Int) -> e -> m) -> Matrix e -> m unit n = fromList [[ if i == j then 1 else 0 | j <- [1..n]] | i <- [1..n] ] zero n = matrix (n,n) (const 0) @@ -442,10 +497,10 @@ , j <- [1..numCols m] , p (i,j) ] - at m (i, j) = ((!! j) . (!! i) . toList) m + at mat (i, j) = ((!! j) . (!! i) . toList) mat - row i m = ((!! (i-1)) . toList) m - col i m = (row i . transpose) m + row i = (!! (i-1)) . toList + col i = row i . transpose numRows = fst . dimensions numCols = snd . dimensions @@ -458,24 +513,28 @@ trace = select (uncurry (==)) inv _ = Nothing - minorMatrix m (i,j) = matrix (numRows m - 1, numCols m - 1) $ - \(i',j') -> m `at` (if i' >= i then i' + 1 else i', - if j' >= j then j' + 1 else j') + minorMatrix (i,j) mat = matrix (numRows mat - 1, numCols mat - 1) $ + \(i',j') -> mat `at` (if i' >= i then i' + 1 else i', + if j' >= j then j' + 1 else j') - minor m = det . minorMatrix m + minor ix = det . minorMatrix ix - cofactors m = matrix (dimensions m) $ - \(i,j) -> fromIntegral ((-1 :: Int)^(i+j)) * minor m (i,j) + cofactors mat = matrix (dimensions mat) $ + \(i,j) -> fromIntegral ((-1 :: Int)^(i+j)) * minor (i,j) mat map f = mapWithIndex (const f) all f = allWithIndex (const f) any f = anyWithIndex (const f) + sum = getSum . foldMap Sum + foldMap f = foldMapWithIndex (const f) mapWithIndex f m = matrix (dimensions m) (\x -> f x (m `at` x)) allWithIndex f m = P.all id [ f (i, j) (m `at` (i,j)) | i <- [1..numRows m], j <- [1..numCols m]] anyWithIndex f m = P.any id [ f (i, j) (m `at` (i,j)) | i <- [1..numRows m], j <- [1..numCols m]] + foldMapWithIndex f m = mconcat [ f (i, j) (m `at` (i,j)) + | i <- [1..numRows m], j <- [1..numCols m]] a `plus` b | dimensions a /= dimensions b = error "Matrix.plus: dimensions don't match." @@ -493,7 +552,7 @@ fromList = _fromList IntMatrix at (IntMatrix _ _ arr) = _at arr - dimensions (IntMatrix m n _) = (m, n) + dimensions (IntMatrix m n _) = (m, n) row i (IntMatrix _ _ arr) = _row i arr col j (IntMatrix _ _ arr) = _col j arr toList (IntMatrix _ _ arr) = _toList arr @@ -505,7 +564,7 @@ fromList = _fromList IntegerMatrix at (IntegerMatrix _ _ arr) = _at arr - dimensions (IntegerMatrix m n _) = (m, n) + dimensions (IntegerMatrix m n _) = (m, n) row i (IntegerMatrix _ _ arr) = _row i arr col j (IntegerMatrix _ _ arr) = _col j arr toList (IntegerMatrix _ _ arr) = _toList arr @@ -524,7 +583,7 @@ det (FloatMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr) rank (FloatMatrix _ _ arr) = runST (_rank thawsBoxed arr) inv (FloatMatrix m n arr) = if m /= n then Nothing else - let x = runST (_inv unboxedST arr) + let x = runST (_inv unboxedST pivotMax arr) in maybe Nothing (Just . FloatMatrix m n) x instance MatrixElement Double where @@ -536,11 +595,11 @@ row i (DoubleMatrix _ _ arr) = _row i arr col j (DoubleMatrix _ _ arr) = _col j arr toList (DoubleMatrix _ _ arr) = _toList arr - inv (DoubleMatrix m n arr) = if m /= n then Nothing else - let x = runST (_inv unboxedST arr) - in maybe Nothing (Just . DoubleMatrix m n) x det (DoubleMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr) rank (DoubleMatrix _ _ arr) = runST (_rank thawsBoxed arr) + inv (DoubleMatrix m n arr) = if m /= n then Nothing else + let x = runST (_inv unboxedST pivotMax arr) + in maybe Nothing (Just . DoubleMatrix m n) x instance (Show a, Integral a) => MatrixElement (Ratio a) where matrix d g = runST (_matrix RatioMatrix arrayST arrayST d g) @@ -551,11 +610,11 @@ row i (RatioMatrix _ _ arr) = _row i arr col j (RatioMatrix _ _ arr) = _col j arr toList (RatioMatrix _ _ arr) = _toList arr - inv (RatioMatrix m n arr) = if m /= n then Nothing else - let x = runST (_inv boxedST arr) - in maybe Nothing (Just . RatioMatrix m n) x det (RatioMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr) rank (RatioMatrix _ _ arr) = runST (_rank thawsBoxed arr) + inv (RatioMatrix m n arr) = if m /= n then Nothing else + let x = runST (_inv boxedST pivotMax arr) + in maybe Nothing (Just . RatioMatrix m n) x instance (Show a, RealFloat a) => MatrixElement (Complex a) where matrix d g = runST (_matrix ComplexMatrix arrayST arrayST d g) @@ -566,11 +625,11 @@ row i (ComplexMatrix _ _ arr) = _row i arr col j (ComplexMatrix _ _ arr) = _col j arr toList (ComplexMatrix _ _ arr) = _toList arr --- 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) rank (ComplexMatrix _ _ arr) = runST (_rank thawsBoxed arr) + inv (ComplexMatrix m n arr) = if m /= n then Nothing else + let x = runST (_inv boxedST pivotNonZero arr) + in maybe Nothing (Just . ComplexMatrix m n) x _at :: (IArray a (u Int e), IArray u e) @@ -647,14 +706,23 @@ a Int (a1 Int b) -> Int -> Int -> m b read a i j = readArray a i >>= flip readArray j +pivotMax :: Ord v => [(i, v)] -> i +pivotMax = fst . L.maximumBy (compare `on` snd) -_inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Ord e, Show e) +pivotNonZero :: (Num v, Eq v) => [(i, v)] -> i +pivotNonZero xs = maybe (fst $ head xs) fst $ L.find ((/= 0) . snd) xs + +_inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Show e, Eq e) => ((Int, Int) -> [e] -> ST s ((u s) Int e)) + -- ^ A function for building a new array + -> ([(Int, e)] -> Int) + -- ^ A function for selecting pivot elements -> Array Int (a Int e) + -- ^ A matrix as arrays or arrays -> ST s (Maybe (Array Int (a Int e))) -_inv mkArrayST mat = do +_inv mkArrayST selectPivot mat = do let m = snd $ bounds mat - n = 2*m + n = 2 * m swap a i j = do tmp <- readArray a i @@ -666,8 +734,8 @@ a <- augment mkArrayST mat forM_ [1..m] $ \k -> do - iPivot <- zip [k..m] <$> mapM (\i -> abs <$> read a i k) [k..m] - >>= return . fst . L.maximumBy (compare `on` snd) + iPivot <- selectPivot <$> zip [k..m] + <$> mapM (\i -> abs <$> read a i k) [k..m] p <- read a iPivot k if p == 0 then writeSTRef okay False else do @@ -713,7 +781,9 @@ _rank :: (IArray a e, MArray (u s) e (ST s), Num e, Division e, Eq e) => (Array Int (a Int e) -> ST s [(u s) Int e]) + -- ^ A function for thawing a boxed array -> Array Int (a Int e) + -- ^ A matrix given as array of arrays -> ST s e _rank thaws mat = do let m = snd $ bounds mat @@ -810,6 +880,9 @@ liftM2 (*) (readSTRef pivotR) (readSTRef signR) +-- TODO: More efficient implementation (decrease the constant factors +-- a little bit by working in the ST monad) +-- [ remark: not sure if this will be faster than lists -> benchmark! ] _mult :: MatrixElement e => Matrix e -> Matrix e -> Matrix e _mult a b = let rowsA = numRows a rowsB = numRows b diff --git a/src/Numeric/Vector.hs b/src/Numeric/Vector.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Vector.hs @@ -0,0 +1,61 @@ +{-# LANGUAGE Haskell2010 + , TypeFamilies + , FlexibleContexts + , Trustworthy + , StandaloneDeriving + , DeriveDataTypeable + , CPP + #-} +{-# OPTIONS -Wall -fno-warn-name-shadowing #-} + +module Numeric.Vector ( + + Vector + +) where + +import Data.Ratio +import Data.Complex +import Data.Int + +import Data.Array.IArray +import Data.Array.Unboxed + +import Data.Typeable + +import qualified Data.Array.Unsafe as U + +data family Vector e + +#if defined(__GLASGOW_HASKELL__) && (__GLASGOW_HASKELL__ >= 707) +deriving instance Typeable Vector +#else +deriving instance Typeable1 Vector +#endif + + +data instance Vector Int + = IntVector !Int (UArray Int Int) + +data instance Vector Float + = FloatVector !Int (UArray Int Float) + +data instance Vector Double + = DoubleVector !Int (UArray Int Double) + +data instance Vector Integer + = IntegerVector !Int (Array Int Integer) + +data instance Vector (Ratio a) + = RatioVector !Int (Array Int (Ratio a)) + +data instance Vector (Complex a) + = ComplexVector !Int (Array Int (Complex a)) + +-- <.> +dotProd = undefined + +-- >< +vectorProd = undefined + +