bed-and-breakfast 0.2.3 → 0.3

raw patch · 3 files changed

+102/−28 lines, 3 filesdep +QuickCheckdep +bed-and-breakfastPVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck, bed-and-breakfast

API changes (from Hackage documentation)

Files

• bed-and-breakfast.cabal +11/−2
• quick-check-tests.hs +36/−0
• src/Numeric/Matrix.hs +55/−26

@@ -1,5 +1,5 @@ Name:           bed-and-breakfast-Version:        0.2.3+Version:        0.3 Synopsis:       Efficient Matrix operations in 100% Haskell. Description:    Efficient Matrix operations in 100% Haskell.                 .@@ -37,7 +37,10 @@                 .                 [@v0.2.3@] Added 'Read' instance for @Matrix@.                     Improved on documentation.-+                .+                [@v0.3@] Added a QuickCheck test suite, fixed a bug in @det@+                    (det would crash for singular matrices, where it should+                    return 0; it does so no).  License:        MIT License-File:   LICENSE@@ -60,4 +63,10 @@                         array >= 0.4     Hs-Source-Dirs:     src +Test-Suite quickcheck+    type:           exitcode-stdio-1.0+    main-is:        quick-check-tests.hs+    build-depends:  base >= 4.5 && < 5,+                    bed-and-breakfast == 0.3,+                    QuickCheck >= 2.4.2 

@@ -0,0 +1,36 @@+{-# LANGUAGE Haskell2010, TemplateHaskell #-}++import Control.Monad++import Numeric.Matrix++import Test.QuickCheck+import Test.QuickCheck.All++import System.Exit+++dim :: Num a => a+dim = 4++instance (MatrixElement e, Arbitrary e) => Arbitrary (Matrix e) where++    arbitrary = sequence (replicate dim (vector dim)) >>= return . fromList+++prop_det :: Matrix Rational -> Bool+prop_det m1 = let m1' = inv m1 in case m1' of (Just _) -> det m1 /= 0; _ -> det m1 == 0++prop_inv :: Matrix Rational -> Bool+prop_inv m1 = let m1' = inv m1 in case m1' of (Just m1') -> m1' * m1 == unit dim; _ -> True++prop_rank :: Matrix Rational -> Bool+prop_rank m1 = let r = rank m1 in if r < dim then det m1 == 0 else r == dim+++main = do+    success <- $(quickCheckAll)+    +    (if success then exitSuccess else exitFailure)++

@@ -290,22 +290,48 @@      -- | The dimensions of a given matrix.     dimensions :: Matrix e -> (Int, Int)++    -- | The number of rows in the given matrix.     numRows :: Matrix e -> Int++    -- | The number of columns in the given matrix.     numCols :: Matrix e -> Int       -- | Builds a matrix from a list of lists.     --+    -- The innermost lists represent the rows. This function will create a m-n-matrix,+    -- where m is the number of rows, which is the minimum length of the row lists+    -- and n is the number of columns, i.e. the length of the outer list.+    --     -- > 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 4+    -- 1 0 0 0+    -- 0 1 0 0+    -- 0 0 1 0+    -- 0 0 0 1     unit  :: Int -> Matrix e      -- | A square matrix of the given size consisting of all zeros.+    -- +    -- >> zero 2+    -- 0 0+    -- 0 0     zero  :: Int -> Matrix e +    -- | A square matrix which trace is the given list, all other components+    -- set to zero.+    --+    -- >> diag [1,4,7,9]+    -- 1 0 0 0+    -- 0 4 0 0+    -- 0 0 7 0+    -- 0 0 0 9     diag  :: [e] -> Matrix e      -- | Check whether the matrix is the empty matrix.@@ -318,9 +344,6 @@     times :: Matrix e -> Matrix e -> Matrix e     inv   :: Matrix e -> Maybe (Matrix e) ---    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.@@ -376,7 +399,7 @@                                     (x:xs) -> (length xs + 1, length x)      adjugate = transpose . cofactors-    transpose = fromList . L.transpose . toList+    transpose m = matrix (dimensions m) (\(i,j) -> m `at` (j,i))     trace = select (uncurry (==))     inv _ = Nothing @@ -703,34 +726,40 @@     pivotR <- newSTRef 1      flip mapM_ [1..size] $ \k -> do+        sign <- readSTRef signR+        unless (sign == 0) $ do -        prev  <- readSTRef pivotR-        pivot <- read a k k >>= tee (writeSTRef pivotR)+            prev  <- readSTRef pivotR+            pivot <- read a k k >>= tee (writeSTRef pivotR) -        when (pivot == 0) $ do-            s <- flip mapM [(k+1)..size] $ \r -> do-                a_rk <- read a r k-                if a_rk == 0 then return 0 else return r-            let sf = filter (>0) s+            when (pivot == 0) $ do+                s <- flip mapM [(k+1)..size] $ \r -> do+                    a_rk <- read a r k+                    if a_rk == 0 then return 0 else return r+                let sf = filter (>0) s -            when (not $ null sf) $ do-                let sw = head sf+                when (not $ null sf) $ do+                    let sw = head sf -                row <- readArray a sw-                readArray a k >>= writeArray a sw-                writeArray a k row+                    row <- readArray a sw+                    readArray a k >>= writeArray a sw+                    writeArray a k row -                read a k k >>= writeSTRef pivotR-                readSTRef signR >>= writeSTRef signR . negate+                    read a k k >>= writeSTRef pivotR+                    readSTRef signR >>= writeSTRef signR . negate -        pivot' <- readSTRef pivotR-        flip mapM [(k+1)..size] $ \i -> do-            a_i <- readArray a i-            flip mapM [(k+1)..size] $ \j -> do-                a_ij <- readArray a_i j-                a_ik <- readArray a_i k-                a_kj <- read a k j-                writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)+                when (null sf) (writeSTRef signR 0)++            sign' <- readSTRef signR+            unless (sign' == 0) $ do+                pivot' <- readSTRef pivotR+                flip mapM_ [(k+1)..size] $ \i -> do+                    a_i <- readArray a i+                    flip mapM [(k+1)..size] $ \j -> do+                        a_ij <- readArray a_i j+                        a_ik <- readArray a_i k+                        a_kj <- read a k j+                        writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)      liftM2 (*) (readSTRef pivotR) (readSTRef signR)