bed-and-breakfast 0.3.2 → 0.4

raw patch · 2 files changed

+114/−131 lines, 2 filesdep +binaryPVP ok

version bump matches the API change (PVP)

Dependencies added: binary

API changes (from Hackage documentation)

Files

• bed-and-breakfast.cabal +7/−1
• src/Numeric/Matrix.hs +107/−130

@@ -1,5 +1,5 @@ Name:           bed-and-breakfast-Version:        0.3.2+Version:        0.4 Synopsis:       Efficient Matrix operations in 100% Haskell. Description:    Efficient Matrix operations in 100% Haskell.                 .@@ -48,6 +48,10 @@                 .                 [@v0.3.2@] Numeric.Matrix.Sugar was not mentioned in the                     cabal file. Improved test suite. Improved documentation.+                .+                [@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.   License:        MIT@@ -70,6 +74,7 @@     Build-Depends:      base >= 4.5 && < 5,                         deepseq >= 1.3,                         array >= 0.4,+                        binary >= 0.5,                         template-haskell >= 2.7     Hs-Source-Dirs:     src @@ -80,4 +85,5 @@     build-depends:  base >= 4.5 && < 5,                     bed-and-breakfast == 0.3.2,                     QuickCheck >= 2.4.2+ 

@@ -59,26 +59,22 @@ import Control.Monad import Control.Monad.ST ---import Control.Concurrent---import Control.Concurrent.MVar---import qualified Control.Monad.Parallel as MPar- import Data.Function (on) import Data.Ratio import Data.Complex import Data.Maybe---import Data.Foldable (Foldable)---import qualified Data.Foldable as F+import Data.Int+import Data.Word  import qualified Data.List as L import Data.Array.IArray import Data.Array.MArray import Data.Array.Unboxed import Data.Array.ST---import Data.Array.IO+import Data.Array.Base (unsafeRead, unsafeWrite) import Data.STRef+import Data.Binary ---import System.IO.Unsafe import qualified Data.Array.Unsafe as U  import Data.Typeable@@ -86,6 +82,7 @@ import Prelude hiding (any, all, read, map) import qualified Prelude as P + -- | Matrices are represented by a type which fits best the component type. -- For example a @Matrix Double@ is represented by unboxed arrays, -- @Matrix Integer@ by boxed arrays.@@ -95,9 +92,9 @@ -- example you can not compute the inverse matrix of a @Matrix Int@. -- -- Every matrix (regardless of the component type) has instances for--- 'Show', 'Read', 'Num', 'Fractional', 'Eq', 'Typeable', and 'NFData'.--- This means that you can use arithmetic operations like '+', '*', and--- '/', as well as functions like 'show', 'read', or 'typeOf'.+-- 'Show', 'Read', 'Num', 'Fractional', 'Eq', 'Typeable', 'Binary',+-- and 'NFData'. This means that you can use arithmetic operations like+-- '+', '*', and '/', as well as functions like 'show', 'read', or 'typeOf'. -- -- [@Show (Matrix e)@] -- Note that a Show instance for the component type @e@ must exist.@@ -132,7 +129,11 @@ -- monad (see the @monad-par@ package for details). -- -- [@Typeable (Matrix e)@]--- Allows you to use matrices as 'Data.Dynamic' values.+-- Allows you to use matrices as 'Dynamic' values.+--+-- [@Binary (Matrix e)@]+-- Serialize and unserialize matrices using the @binary@ package.+-- See @encode@ and @decode@. data family Matrix e  data instance Matrix Int@@ -191,7 +192,25 @@ instance (MatrixElement e) => NFData (Matrix e) where     rnf matrix = matrix `deepseq` () +instance (MatrixElement e, Binary e) => Binary (Matrix e) where+    +    put m = do+        let (rows, cols) = dimensions m+        +        put rows >> put cols +        forM_ [1..rows] $ \i -> do+            forM_ [1..cols] $ \j -> do+                put (m `at` (i,j))++    get = do+        rows <- get :: Get Int+        cols <- get :: Get Int++        forM [1..rows] (const (forM [1..cols] (const get)))+            >>= return . fromList++ (<|>) :: MatrixElement e => Matrix e -> Matrix e -> Matrix e -- ^ Joins two matrices horizontally. --@@ -268,19 +287,22 @@ class Division e where     divide :: e -> e -> e -instance Division Int    where divide = quot--- instance Division Int8   where divide = quot--- instance Division Int16  where divide = quot--- instance Division Int32  where divide = quot--- instance Division Int64  where divide = quot+instance Division Int     where divide = quot+instance Division Int8    where divide = quot+instance Division Int16   where divide = quot+instance Division Int32   where divide = quot+instance Division Int64   where divide = quot+instance Division Word8   where divide = quot+instance Division Word16  where divide = quot+instance Division Word32  where divide = quot+instance Division Word64  where divide = quot instance Division Integer where divide = quot-instance Division Float  where divide = (/)-instance Division Double where divide = (/)+instance Division Float   where divide = (/)+instance Division Double  where divide = (/) instance Integral a => Division (Ratio a) where divide = (/) instance RealFloat a => Division (Complex a) where divide = (/)  - class (Eq e, Num e) => MatrixElement e where      -- | Creates a matrix of the given size using a generator@@ -366,12 +388,14 @@     -- Returns 0 if the matrix is not a square matrix.     det       :: Matrix e -> e -    -- | Flips rows and columns.+    -- | Flip 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++    -- | Compute the rank of a matrix.     rank      :: Matrix e -> e     trace     :: Matrix e -> [e] @@ -380,15 +404,15 @@     adjugate :: MatrixElement e => Matrix e -> Matrix e     minorMatrix :: MatrixElement e => Matrix e -> (Int, Int) -> Matrix e -    -- | Applies a function on every component in the matrix.+    -- | Apply 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.+    -- | Apply a predicate on every component in the matrix+    -- and returns True iff 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.+    -- | Apply a predicate on every component in the matrix+    -- and return 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@@ -551,7 +575,7 @@           => (Int -> Int -> a Int (u Int e) -> matrix e) -> [[e]] -> matrix e _fromList c xs =     let lengths = P.map length xs-        numCols = foldl1 min lengths+        numCols = if null lengths then 0 else foldl1 min lengths         numRows = length lengths              in  c numRows numCols@@ -602,21 +626,13 @@          => (Int, Int) -> e -> ST s ((STUArray s) Int e) arraySTU = newArray ---arrayIO :: MArray IOArray e IO---        => (Int, Int) -> e -> IO (IOArray Int e)---arrayIO = newArray ---arrayIOU :: MArray IOUArray e IO---         => (Int, Int) -> e -> IO (IOUArray Int e)---arrayIOU = newArray-- tee :: Monad m => (b -> m a) -> b -> m b tee f x = f x >> return x  read :: (MArray a1 b m, MArray a (a1 Int b) m) =>                        a Int (a1 Int b) -> Int -> Int -> m b-read a i j = readArray a i >>= flip readArray j+read a i j = unsafeRead a i >>= flip unsafeRead j   _inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Ord e, Show e)@@ -628,15 +644,15 @@         n = 2*m          swap a i j = do-            tmp <- readArray a i-            readArray a j >>= writeArray a i-            writeArray a j tmp+            tmp <- unsafeRead a i+            unsafeRead a j >>= unsafeWrite a i+            unsafeWrite a j tmp      okay <- newSTRef True      a <- augment mkArrayST mat -    flip mapM_ [1..m] $ \k -> do+    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) @@ -645,38 +661,38 @@              swap a iPivot k -            flip mapM_ [k+1..m] $ \i -> do-                a_i <- readArray a i-                a_k <- readArray a k-                flip mapM_ [k+1..n] $ \j -> do-                    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+            forM_ [k+1..m] $ \i -> do+                a_i <- unsafeRead a i+                a_k <- unsafeRead 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      invertible <- readSTRef okay      if invertible then       do-        flip mapM_ [ m - v | v <- [0..m-1] ] $ \i -> do-            a_i <- readArray a i-            p   <- readArray a_i i-            writeArray a_i i 1-            flip mapM_ [i+1..n] $ \j -> do-                readArray a_i j >>= writeArray a_i j . (/ p)+        forM_ [ m - v | v <- [0..m-1] ] $ \i -> do+            a_i <- unsafeRead a i+            p   <- unsafeRead a_i i+            unsafeWrite a_i i 1+            forM_ [i+1..n] $ \j -> do+                unsafeRead a_i j >>= unsafeWrite a_i j . (/ p)              unless (i == m) $ do-                flip mapM_ [i+1..m] $ \k -> do-                    a_k <- readArray a k-                    p   <- readArray a_i k+                forM_ [i+1..m] $ \k -> do+                    a_k <- unsafeRead a k+                    p   <- unsafeRead a_i k -                    flip mapM_ [k..n] $ \j -> do-                        a_ij <- readArray a_i j-                        a_kj <- readArray a_k j-                        writeArray a_i j (a_ij - p * a_kj)+                    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) -        mapM (\i -> readArray a i >>= getElems+        mapM (\i -> unsafeRead a i >>= getElems                         >>= return . listArray (1, m) . drop m) [1..m]             >>= return . Just . listArray (1, m) @@ -691,16 +707,16 @@         n = snd $ bounds (mat ! 1)          swap a i j = do-            tmp <- readArray a i-            readArray a j >>= writeArray a i-            writeArray a j tmp+            tmp <- unsafeRead a i+            unsafeRead a j >>= unsafeWrite a i+            unsafeWrite a j tmp      a <- thaws mat >>= arrays      ixPivot <- newSTRef 1     prevR   <- newSTRef 1 -    flip mapM_ [1..n] $ \k -> do+    forM_ [1..n] $ \k -> do         pivotRow <- readSTRef ixPivot          switchRow <- mapM (\i -> read a i k) [pivotRow .. m]@@ -710,17 +726,17 @@             let ix = fromJust switchRow + pivotRow             when (pivotRow /= ix) (swap a pivotRow ix) -            a_p   <- readArray a k-            pivot <- readArray a_p k+            a_p   <- unsafeRead a k+            pivot <- unsafeRead a_p k             prev  <- readSTRef prevR             -            flip mapM_ [pivotRow+1..m] $ \i -> do-                a_i <- readArray a i-                flip mapM_ [k+1..n] $ \j -> do-                    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)+            forM_ [pivotRow+1..m] $ \i -> do+                a_i <- unsafeRead 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)                                         `divide` prev)              writeSTRef ixPivot (pivotRow + 1)@@ -742,7 +758,7 @@     signR  <- newSTRef 1     pivotR <- newSTRef 1 -    flip mapM_ [1..size] $ \k -> do+    forM_ [1..size] $ \k -> do         sign <- readSTRef signR         unless (sign == 0) $ do @@ -750,7 +766,7 @@             pivot <- read a k k >>= tee (writeSTRef pivotR)              when (pivot == 0) $ do-                s <- flip mapM [(k+1)..size] $ \r -> do+                s <- forM [(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@@ -758,9 +774,9 @@                 when (not $ null sf) $ do                     let sw = head sf -                    row <- readArray a sw-                    readArray a k >>= writeArray a sw-                    writeArray a k row+                    row <- unsafeRead a sw+                    unsafeRead a k >>= unsafeWrite a sw+                    unsafeWrite a k row                      read a k k >>= writeSTRef pivotR                     readSTRef signR >>= writeSTRef signR . negate@@ -770,13 +786,13 @@             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+                forM_ [(k+1)..size] $ \i -> do+                    a_i <- unsafeRead a i+                    forM [(k+1)..size] $ \j -> do+                        a_ij <- unsafeRead a_i j+                        a_ik <- unsafeRead a_i k                         a_kj <- read a k j-                        writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)+                        unsafeWrite a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)      liftM2 (*) (readSTRef pivotR) (readSTRef signR) @@ -787,21 +803,7 @@                 colsB = numCols b             in  matrix (rowsA, colsB) (\(i,j) -> L.foldl' (+) 0 [a `at` (i, k) * b `at` (k, j) | k <- [1..rowsB]]) -{- Old function: TODO benchmark against new one -_matrix_ :: (IArray a (u Int e), IArray u e)-        => (Int -> Int -> (a Int (u Int e)) -> matrix e)-        -> (Int, Int)-        -> ((Int, Int) -> e)-        -> matrix e-_matrix_ c (numRows, numCols) generator =-    c numRows numCols-      $ array (1, numRows)-      $ [ (i, array (1, numCols) [(j, generator (i, j))-                                 | j <- [1..numCols]])-          | i <- [1..numRows] ]--}- _matrix :: (IArray a1 (u Int e), IArray u e,             MArray a2 (u Int e) (ST s), MArray a3 e (ST s),             Num e)@@ -813,38 +815,13 @@         -> ST s matrix _matrix c newArray newArrayU (m, n) g = do     rows <- newArray (1, m) undefined-    flip mapM_ [1..m] $ \i -> do+    forM_ [1..m] $ \i -> do         cols <- newArrayU (1, n) 0-        flip mapM_ [1..n] $ \j -> do-            writeArray cols j (g (i,j))-        U.unsafeFreeze cols >>= writeArray rows i+        forM_ [1..n] $ \j -> do+            unsafeWrite cols j (g (i,j))+        U.unsafeFreeze cols >>= unsafeWrite rows i     U.unsafeFreeze rows >>= return . c m n -{--_matrixIO c newArray newArrayU (m, n) g = do-    rows <- newArray (1, m) undefined--    let half = m `div` 2--    wait <- newEmptyMVar-    forkOS $ do-        flip mapM_ [1..half] $ \i -> do-            cols <- newArrayU (1, n) 0-            flip mapM_ [1..n] $ \j -> do-                writeArray cols j (g (i,j))-            U.unsafeFreeze cols >>= writeArray rows i-        putMVar wait ()--    flip mapM_ [half+1..m] $ \i -> do-        cols <- newArrayU (1, n) 0-        flip mapM_ [1..n] $ \j -> do-            writeArray cols j (g (i,j))-        U.unsafeFreeze cols >>= writeArray rows i--    takeMVar wait--    U.unsafeFreeze rows >>= return . c m n--}  {-# RULES