bed-and-breakfast 0.4 → 0.4.1

raw patch · 3 files changed

+137/−55 lines, 3 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• bed-and-breakfast.cabal +6/−1
• quickcheck-tests.hs +38/−6
• src/Numeric/Matrix.hs +93/−48

@@ -1,5 +1,5 @@ Name:           bed-and-breakfast-Version:        0.4+Version:        0.4.1 Synopsis:       Efficient Matrix operations in 100% Haskell. Description:    Efficient Matrix operations in 100% Haskell.                 .@@ -52,6 +52,11 @@                 [@v0.4@] Fixed a bug regarding @empty@ and @fromList@.                     Use unsafe operations where it is safe for speed.                     Added RULES. Added an instance for binary.+                .+                [@v0.4.1@] The unsafe operations used in v0.4 turned out+                    to fatally fail on certain platforms. Revoked this+                    change. Thanks to Marcin Fatyga for pointing at the+                    issue.   License:        MIT

@@ -1,5 +1,8 @@ {-# LANGUAGE Haskell2010, TemplateHaskell #-} +import Prelude+import qualified Prelude as P+ import Control.Monad  import Numeric.Matrix@@ -9,7 +12,11 @@  import System.Exit +import Data.Ratio+import Data.Int+import Data.Word + dim :: Num a => a dim = 6 @@ -18,12 +25,40 @@     arbitrary = sequence (replicate dim (vector dim)) >>= return . fromList  +basic m = m + m == scale m 2 && m * m * m == m ^ 3++prop_zero :: Word8 -> Bool+prop_zero = isZero . (zero :: Int -> Matrix Double) . fromIntegral++prop_sum_trace :: Word8 -> Bool+prop_sum_trace n = ((== nd) . P.sum . trace . (unit :: Int -> Matrix Double) . fromIntegral) n'+  where n' = n + 1+        nd = fromIntegral n'++prop_basic_integer :: Matrix Integer -> Bool+prop_basic_integer = basic++prop_basic_int :: Matrix Int -> Bool+prop_basic_int = basic++prop_basic_float :: Matrix Float -> Bool+prop_basic_float = basic++prop_basic_double :: Matrix Double -> Bool+prop_basic_double = basic++prop_basic_rational :: Matrix Rational -> Bool+prop_basic_rational = basic++prop_zero_det_double :: Matrix Double -> Bool+prop_zero_det_double m1 = let m1' = inv m1 in case m1' of (Just _) -> det m1 /= 0; _ -> det m1 == 0++prop_zero_det_rational :: Matrix Rational -> Bool+prop_zero_det_rational m1 = let m1' = inv m1 in case m1' of (Just _) -> det m1 /= 0; _ -> det m1 == 0+ prop_plus_commutative :: Matrix Double -> Matrix Double -> Bool prop_plus_commutative m1 m2 = m1 + m2 == m2 + m1 -prop_zero_det :: Matrix Rational -> Bool-prop_zero_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 @@ -39,11 +74,8 @@ prop_inv_twice :: Matrix Rational -> Bool prop_inv_twice m1 = let m1' = inv m1 in case m1' of (Just m1') -> inv m1' == Just m1; _ -> True - main = do-     success <- $(quickCheckAll)-         (if success then exitSuccess else exitFailure)  

@@ -71,7 +71,6 @@ import Data.Array.MArray import Data.Array.Unboxed import Data.Array.ST-import Data.Array.Base (unsafeRead, unsafeWrite) import Data.STRef import Data.Binary @@ -136,9 +135,16 @@ -- See @encode@ and @decode@. data family Matrix e + data instance Matrix Int     = IntMatrix !Int !Int (Array Int (UArray Int Int)) +data instance Matrix Int32+    = Int32Matrix !Int !Int (Array Int (UArray Int Int32))++data instance Matrix Int64+    = Int64Matrix !Int !Int (Array Int (UArray Int Int64))+ data instance Matrix Float     = FloatMatrix !Int !Int (Array Int (UArray Int Float)) @@ -189,8 +195,8 @@             = allWithIndex (\ix e -> m `at` ix == e) n         | otherwise = False -instance (MatrixElement e) => NFData (Matrix e) where-    rnf matrix = matrix `deepseq` ()+instance (MatrixElement e, NFData e) => NFData (Matrix e) where+    rnf matrix = map (\x -> x `deepseq` x) matrix `deepseq` ()  instance (MatrixElement e, Binary e) => Binary (Matrix e) where     @@ -378,9 +384,23 @@     -- > dimensions empty == (0, 0)     empty :: Matrix e +    -- | Subtract a matrix from another.     minus :: Matrix e -> Matrix e -> Matrix e++    -- | Add two matrices.+    --+    -- You may also use the 'Num' instance for matrices,+    -- i.e. @matrix1 + matrix2@ will work, too.     plus  :: Matrix e -> Matrix e -> Matrix e++    -- | Multiply two matrices /O(n^3)/.+    --+    -- You may also use the 'Num' instance for matrices,+    -- i.e. @matrix1 * matrix2@ will work, too.     times :: Matrix e -> Matrix e -> Matrix e++    -- | Compute the inverse of a matrix. Returns @Nothing@+    -- if the matrix is not invertible.     inv   :: Matrix e -> Maybe (Matrix e)      -- | Applies Bareiss multistep integer-preserving@@ -486,6 +506,30 @@     det        (IntMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)     rank       (IntMatrix _ _ arr) = runST (_rank thawsBoxed arr) +instance MatrixElement Int32 where+    matrix d g = runST (_matrix Int32Matrix arrayST arraySTU d g)+    fromList = _fromList Int32Matrix++    at         (Int32Matrix _ _ arr) = _at arr+    dimensions (Int32Matrix m n _) = (m, n)+    row i      (Int32Matrix _ _ arr) = _row i arr+    col j      (Int32Matrix _ _ arr) = _col j arr+    toList     (Int32Matrix _ _ arr) = _toList arr+    det        (Int32Matrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)+    rank       (Int32Matrix _ _ arr) = runST (_rank thawsBoxed arr)++instance MatrixElement Int64 where+    matrix d g = runST (_matrix Int64Matrix arrayST arraySTU d g)+    fromList = _fromList Int64Matrix++    at         (Int64Matrix _ _ arr) = _at arr+    dimensions (Int64Matrix m n _) = (m, n)+    row i      (Int64Matrix _ _ arr) = _row i arr+    col j      (Int64Matrix _ _ arr) = _col j arr+    toList     (Int64Matrix _ _ arr) = _toList arr+    det        (Int64Matrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)+    rank       (Int64Matrix _ _ arr) = runST (_rank thawsBoxed arr)+ instance MatrixElement Integer where     matrix d g = runST (_matrix IntegerMatrix arrayST arrayST d g)     fromList   = _fromList IntegerMatrix@@ -530,7 +574,6 @@  instance (Show a, Integral a) => MatrixElement (Ratio a) where     matrix d g = runST (_matrix RatioMatrix arrayST arrayST d g)---    matrix d g = unsafePerformIO (_matrixIO RatioMatrix arrayIO arrayIO d g)     fromList   = _fromList RatioMatrix      at         (RatioMatrix _ _ arr) = _at arr@@ -632,7 +675,7 @@  read :: (MArray a1 b m, MArray a (a1 Int b) m) =>                        a Int (a1 Int b) -> Int -> Int -> m b-read a i j = unsafeRead a i >>= flip unsafeRead j+read a i j = readArray a i >>= flip readArray j   _inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Ord e, Show e)@@ -644,9 +687,9 @@         n = 2*m          swap a i j = do-            tmp <- unsafeRead a i-            unsafeRead a j >>= unsafeWrite a i-            unsafeWrite a j tmp+            tmp <- readArray a i+            readArray a j >>= writeArray a i+            writeArray a j tmp      okay <- newSTRef True @@ -662,37 +705,37 @@             swap a iPivot k              forM_ [k+1..m] $ \i -> do-                a_i <- unsafeRead a i-                a_k <- unsafeRead a k+                a_i <- readArray a i+                a_k <- readArray a k                 forM_ [k+1..n] $ \j -> do-                    a_ij <- unsafeRead a_i j-                    a_kj <- unsafeRead a_k j-                    a_ik <- unsafeRead a_i k-                    unsafeWrite a_i j (a_ij - a_kj * (a_ik / p))-                unsafeWrite a_i k 0+                    a_ij <- readArray a_i j+                    a_kj <- readArray a_k j+                    a_ik <- readArray a_i k+                    writeArray a_i j (a_ij - a_kj * (a_ik / p))+                writeArray a_i k 0      invertible <- readSTRef okay      if invertible then       do         forM_ [ m - v | v <- [0..m-1] ] $ \i -> do-            a_i <- unsafeRead a i-            p   <- unsafeRead a_i i-            unsafeWrite a_i i 1+            a_i <- readArray a i+            p   <- readArray a_i i+            writeArray a_i i 1             forM_ [i+1..n] $ \j -> do-                unsafeRead a_i j >>= unsafeWrite a_i j . (/ p)+                readArray a_i j >>= writeArray a_i j . (/ p)              unless (i == m) $ do                 forM_ [i+1..m] $ \k -> do-                    a_k <- unsafeRead a k-                    p   <- unsafeRead a_i k+                    a_k <- readArray a k+                    p   <- readArray a_i k                      forM_ [k..n] $ \j -> do-                        a_ij <- unsafeRead a_i j-                        a_kj <- unsafeRead a_k j-                        unsafeWrite a_i j (a_ij - p * a_kj)+                        a_ij <- readArray a_i j+                        a_kj <- readArray a_k j+                        writeArray a_i j (a_ij - p * a_kj) -        mapM (\i -> unsafeRead a i >>= getElems+        mapM (\i -> readArray a i >>= getElems                         >>= return . listArray (1, m) . drop m) [1..m]             >>= return . Just . listArray (1, m) @@ -707,9 +750,9 @@         n = snd $ bounds (mat ! 1)          swap a i j = do-            tmp <- unsafeRead a i-            unsafeRead a j >>= unsafeWrite a i-            unsafeWrite a j tmp+            tmp <- readArray a i+            readArray a j >>= writeArray a i+            writeArray a j tmp      a <- thaws mat >>= arrays @@ -726,17 +769,17 @@             let ix = fromJust switchRow + pivotRow             when (pivotRow /= ix) (swap a pivotRow ix) -            a_p   <- unsafeRead a k-            pivot <- unsafeRead a_p k+            a_p   <- readArray a k+            pivot <- readArray a_p k             prev  <- readSTRef prevR                          forM_ [pivotRow+1..m] $ \i -> do-                a_i <- unsafeRead a i+                a_i <- readArray a i                 forM_ [k+1..n] $ \j -> do-                    a_ij <- unsafeRead a_i j-                    a_ik <- unsafeRead a_i k-                    a_pj <- unsafeRead a_p j-                    unsafeWrite a_i j ((pivot * a_ij - a_ik * a_pj)+                    a_ij <- readArray a_i j+                    a_ik <- readArray a_i k+                    a_pj <- readArray a_p j+                    writeArray a_i j ((pivot * a_ij - a_ik * a_pj)                                         `divide` prev)              writeSTRef ixPivot (pivotRow + 1)@@ -774,9 +817,9 @@                 when (not $ null sf) $ do                     let sw = head sf -                    row <- unsafeRead a sw-                    unsafeRead a k >>= unsafeWrite a sw-                    unsafeWrite 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@@ -787,21 +830,23 @@             unless (sign' == 0) $ do                 pivot' <- readSTRef pivotR                 forM_ [(k+1)..size] $ \i -> do-                    a_i <- unsafeRead a i+                    a_i <- readArray a i                     forM [(k+1)..size] $ \j -> do-                        a_ij <- unsafeRead a_i j-                        a_ik <- unsafeRead a_i k+                        a_ij <- readArray a_i j+                        a_ik <- readArray a_i k                         a_kj <- read a k j-                        unsafeWrite a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)+                        writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)      liftM2 (*) (readSTRef pivotR) (readSTRef signR)   _mult :: MatrixElement e => Matrix e -> Matrix e -> Matrix e-_mult a b = let rowsA = numRows a-                rowsB = numRows b-                colsB = numCols b-            in  matrix (rowsA, colsB) (\(i,j) -> L.foldl' (+) 0 [a `at` (i, k) * b `at` (k, j) | k <- [1..rowsB]])+_mult a b =+    let rowsA = numRows a+        rowsB = numRows b+        colsB = numCols b+        gen i j = L.foldl' (+) 0 [a `at` (i, k) * b `at` (k, j) | k <- [1..rowsB]]+    in  matrix (rowsA, colsB) (uncurry gen)   _matrix :: (IArray a1 (u Int e), IArray u e,@@ -818,8 +863,8 @@     forM_ [1..m] $ \i -> do         cols <- newArrayU (1, n) 0         forM_ [1..n] $ \j -> do-            unsafeWrite cols j (g (i,j))-        U.unsafeFreeze cols >>= unsafeWrite rows i+            writeArray cols j (g (i,j))+        U.unsafeFreeze cols >>= writeArray rows i     U.unsafeFreeze rows >>= return . c m n