bed-and-breakfast 0.1.2 → 0.1.3

raw patch · 2 files changed

+105/−48 lines, 2 filesdep +deepseqdep ~basePVP ok

version bump matches the API change (PVP)

Dependencies added: deepseq

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• bed-and-breakfast.cabal +9/−3
• src/Numeric/Matrix.hs +96/−45

@@ -1,5 +1,5 @@ Name:           bed-and-breakfast-Version:        0.1.2+Version:        0.1.3 Synopsis:       Efficient Matrix operations in 100% Haskell. Description:    Efficient Matrix operations in 100% Haskell.                 .@@ -8,10 +8,15 @@                     'Float', 'Double', 'Complex', and 'Rational'.                 .                 [@v0.1.1@] Fixed wrong algorithm for computing the-                    inverse of a Matrix.+                    inverse of a 'Matrix'.                 .                 [@v0.1.2@] Added instances for @Num Matrix@,                     @Fractional Matrix@, and @Eq Matrix@.+                .+                [@v0.1.3@] @inv@ is now a total function and will+                    no longer call `error' if a matrix is not+                    invertible. Also @Matrix@ derives 'Data.Typeable'+                    now.  License:        MIT License-File:   LICENSE@@ -19,7 +24,7 @@ Maintainer:     Julian Fleischer <julian.fleischer@fu-berlin.de> Build-Type:     Simple Cabal-Version:  >= 1.8-Category:       Data+Category:       Numeric, Math Stability:      stable Homepage:       http://hub.darcs.net/scravy/bed-and-breakfast @@ -30,6 +35,7 @@ Library     Exposed-Modules:    Numeric.Matrix     Build-Depends:      base >= 4.5 && < 5,+                        deepseq >= 1.3,                         array >= 0.4     Hs-Source-Dirs:     src 

@@ -50,18 +50,24 @@         mapWithIndex,         allWithIndex,         anyWithIndex-    )+    ),+    (<|>),+    (<->),+    scale ) where  -import Control.Applicative+import Control.Applicative ((<$>))+import Control.DeepSeq import Control.Monad import Control.Monad.ST -import Data.Function+import Data.Function (on) import Data.Ratio import Data.Complex import Data.Maybe+import Data.Foldable (Foldable)+import qualified Data.Foldable as F  import qualified Data.List as L import Data.Array.IArray@@ -69,22 +75,33 @@ import Data.Array.Unboxed import Data.Array.ST import Data.STRef+import Data.Typeable -import Prelude hiding (any, all, read)+import Prelude hiding (any, all, read, map) import qualified Prelude as P -import qualified Debug.Trace as D- data family Matrix e -data instance Matrix Int    = IntMatrix    Int Int (Array Int (UArray Int Int))-data instance Matrix Float  = FloatMatrix  Int Int (Array Int (UArray Int Float))-data instance Matrix Double = DoubleMatrix Int Int (Array Int (UArray Int Double))+data instance Matrix Int = IntMatrix !Int !Int (Array Int (UArray Int Int)) -data instance Matrix Integer = IntegerMatrix Int Int (Array Int (Array Int Integer))-data instance Matrix (Ratio a) = RatioMatrix Int Int (Array Int (Array Int (Ratio a)))-data instance Matrix (Complex a) = ComplexMatrix Int Int (Array Int (Array Int (Complex a)))+data instance Matrix Float = FloatMatrix !Int !Int (Array Int (UArray Int Float)) +data instance Matrix Double = DoubleMatrix !Int !Int (Array Int (UArray Int Double))++data instance Matrix Integer = IntegerMatrix !Int !Int (Array Int (Array Int Integer))++data instance Matrix (Ratio a) = RatioMatrix !Int !Int (Array Int (Array Int (Ratio a)))++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@@ -108,7 +125,30 @@             = allWithIndex (\ix e -> m `at` ix == e) n         | otherwise = False +instance (MatrixElement e) => NFData (Matrix e) where+    rnf matrix = matrix `deepseq` () ++m1 <|> m2 = let m = numCols m1+                n1 = numRows m1+                n2 = numRows m2+            in matrix (max n1 n2, m + numCols m2)+              $ \(i,j) -> if j > m+                    then (if i > n2 then 0 else m2 `at` (i,j-m))+                    else (if i > n1 then 0 else m1 `at` (i,j))++m1 <-> m2 = let m = numRows m1+                n1 = numCols m1+                n2 = numCols m2+            in matrix (m + numRows m2, max n1 n2)+              $ \(i,j) -> if i > m+                    then (if j > n2 then 0 else m2 `at` (i-m,j))+                    else (if j > n1 then 0 else m1 `at` (i,j))++scale :: MatrixElement e => Matrix e -> e -> Matrix e+scale m s = map (*s) m++ class Division e where     divide :: e -> e -> e @@ -235,7 +275,7 @@     row i      (IntMatrix _ _ arr) = _row i arr     col j      (IntMatrix _ _ arr) = _col j arr     toList     (IntMatrix _ _ arr) = _toList arr-    inv = undefined -- IntMatrix $ runST (invSTU arr)+    inv _ = Nothing     det        (IntMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)     rank = undefined -- runST (_rank thawsBoxed arr) @@ -248,7 +288,7 @@     row i      (IntegerMatrix _ _ arr) = _row i arr     col j      (IntegerMatrix _ _ arr) = _col j arr     toList     (IntegerMatrix _ _ arr) = _toList arr-    inv = undefined -- IntMatrix $ runST (invSTU arr)+    inv _ = Nothing     det        (IntegerMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr)     rank = undefined -- runST (_rank thawsBoxed arr) @@ -262,7 +302,8 @@     col j      (FloatMatrix _ _ arr) = _col j arr     toList     (FloatMatrix _ _ arr) = _toList arr     inv        (FloatMatrix m n arr) = if m /= n then Nothing else-                                        Just $ FloatMatrix m n $ runST (_inv unboxedST arr)+                                         let x = runST (_inv unboxedST arr)+                                         in maybe Nothing (Just . FloatMatrix m n) x     det        (FloatMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)     rank       (FloatMatrix _ _ arr) = runST (_rank thawsBoxed arr) @@ -276,7 +317,8 @@     col j      (DoubleMatrix _ _ arr) = _col j arr     toList     (DoubleMatrix _ _ arr) = _toList arr     inv        (DoubleMatrix m n arr) = if m /= n then Nothing else-                                         Just $ DoubleMatrix m n $ runST (_inv unboxedST arr)+                                         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) @@ -290,7 +332,8 @@     col j      (RatioMatrix _ _ arr) = _col j arr     toList     (RatioMatrix _ _ arr) = _toList arr     inv        (RatioMatrix m n arr) = if m /= n then Nothing else-                                        Just $ RatioMatrix m n $ runST (_inv boxedST arr)+                                        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) @@ -395,7 +438,7 @@ _inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Ord e, Show e)      => ((Int, Int) -> [e] -> ST s ((u s) Int e))      -> Array Int (a Int e)-     -> ST s (Array Int (a Int e))+     -> ST s (Maybe (Array Int (a Int e))) _inv mkArrayST mat = do     let m = snd $ bounds mat         n = 2*m@@ -405,6 +448,8 @@             readArray a j >>= writeArray a i             writeArray a j tmp +    okay <- newSTRef True+     a <- augment mkArrayST mat      flip mapM_ [1..m] $ \k -> do@@ -412,40 +457,46 @@                     >>= return . fst . L.maximumBy (compare `on` snd)          p <- read a iPivot k-        when (p == 0) (fail "not invertible")-        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+        if p == 0 then writeSTRef okay False else 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)+            swap a iPivot k -        unless (i == m) $ do-            flip mapM_ [i+1..m] $ \k -> do+            flip mapM_ [k+1..m] $ \i -> do+                a_i <- readArray a i                 a_k <- readArray a k-                p   <- readArray a_i k--                flip mapM_ [k..n] $ \j -> do+                flip mapM_ [k+1..n] $ \j -> do                     a_ij <- readArray a_i j                     a_kj <- readArray a_k j-                    writeArray a_i j (a_ij - p * a_kj)+                    a_ik <- readArray a_i k+                    writeArray a_i j (a_ij - a_kj * (a_ik / p))+                writeArray a_i k 0 -    mapM (\i -> readArray a i >>= getElems-                    >>= return . listArray (1, m) . drop m) [1..m]-        >>= return . listArray (1, m)+    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)++            unless (i == m) $ do+                flip mapM_ [i+1..m] $ \k -> do+                    a_k <- readArray a k+                    p   <- readArray 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)++        mapM (\i -> readArray a i >>= getElems+                        >>= return . listArray (1, m) . drop m) [1..m]+            >>= return . Just . listArray (1, m)++      else return Nothing  _rank :: (IArray a e, MArray (u s) e (ST s), Fractional e, Eq e)       => (Array Int (a Int e) -> ST s [(u s) Int e])