bed-and-breakfast 0.4.3 → 0.5
raw patch · 6 files changed
+201/−54 lines, 6 filesdep +cpphsdep ~base
Dependencies added: cpphs
Dependency ranges changed: base
Files
- CHANGES.md +10/−1
- README.md +2/−2
- ROADMAP.md +1/−1
- bed-and-breakfast.cabal +9/−5
- src/Numeric/Matrix.hs +118/−45
- src/Numeric/Vector.hs +61/−0
CHANGES.md view
@@ -100,9 +100,18 @@ <dt>v0.4.3</dt> <dd>- Fixed a bug in @transpose@ that prevented it from+ Fixed a bug in <code>transpose</code> that prevented it from working correctly with non-square matrices. Thanks to @owst@ from @hub.darcs.net@.+ </dd>++ <dt>v0.5</dt>+ <dd>+ <code>inv</code> works now for complex matrices too.+ Added conversion functions <code>toDoubleMatrix</code>,+ <code>toComplexMatrix</code>, <code>toRationalMatrix</code>.+ Changed signature of <code>minor</code> and <code>minorMatrix</code>+ to a more natural format. </dd> </dl>
README.md view
@@ -1,7 +1,7 @@-bed-and-breakfast+bed-and-breakfast [](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
ROADMAP.md view
@@ -22,4 +22,4 @@ * Vector product * Outer product -+* Make it compatible with GHC 7.4 - 7.10
bed-and-breakfast.cabal view
@@ -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
src/Numeric/Matrix.hs view
@@ -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
+ src/Numeric/Vector.hs view
@@ -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++