bed-and-breakfast 0.1 → 0.1.1

raw patch · 2 files changed

+53/−57 lines, 2 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• bed-and-breakfast.cabal +10/−2
• src/Numeric/Matrix.hs +43/−55

@@ -1,8 +1,11 @@ Name:           bed-and-breakfast-Version:        0.1+Version:        0.1.1 Synopsis:       Efficient Matrix operations in 100% Haskell. Description:    Efficient Matrix operations in 100% Haskell.-                +                .+                [@v0.1.1@] Fixed wrong algorithm for computing the+                    inverse of a Matrix.+ License:        MIT License-File:   LICENSE Author:         Julian Fleischer <julian.fleischer@fu-berlin.de>@@ -11,6 +14,11 @@ Cabal-Version:  >= 1.8 Category:       Data Stability:      stable+Homepage:       http://hub.darcs.net/scravy/bed-and-breakfast++Source-Repository head+    type: darcs+    location: hub.darcs.net:bed-and-breakfast  Library     Exposed-Modules:    Numeric.Matrix

@@ -3,7 +3,7 @@     , FlexibleContexts     , Trustworthy  #-}-{-# OPTIONS -Wall -fno-warn-name-shadowing #-}+--{-# OPTIONS -Wall -fno-warn-name-shadowing #-}  module Numeric.Matrix (     Matrix,@@ -58,6 +58,7 @@ import Control.Monad import Control.Monad.ST +import Data.Function import Data.Ratio import Data.Complex import qualified Data.List as L@@ -70,6 +71,7 @@ import Prelude hiding (any, all, read) import qualified Prelude as P +import qualified Debug.Trace as D  data family Matrix e @@ -258,7 +260,7 @@     det        (DoubleMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)     rank       (DoubleMatrix _ _ arr) = runST (_rank thawsBoxed arr) -instance Integral a => MatrixElement (Ratio a) where+instance (Show a, Integral a) => MatrixElement (Ratio a) where     matrix   = _matrix RatioMatrix     fromList = _fromList RatioMatrix @@ -272,7 +274,7 @@     det        (RatioMatrix m n arr) = if m /= n then 0 else  runST (_det thawsBoxed arr)     rank       (RatioMatrix _ _ arr) = runST (_rank thawsBoxed arr) -instance RealFloat a => MatrixElement (Complex a) where+instance (Show a, RealFloat a) => MatrixElement (Complex a) where     matrix   = _matrix ComplexMatrix     fromList = _fromList ComplexMatrix @@ -281,8 +283,9 @@     row i      (ComplexMatrix _ _ arr) = _row i arr     col j      (ComplexMatrix _ _ arr) = _col j arr     toList     (ComplexMatrix _ _ arr) = _toList arr-    inv        (ComplexMatrix m n arr) = if m /= n then Nothing else-                                          Just $ ComplexMatrix m n $ runST (_inv boxedST arr)+    inv        (ComplexMatrix m n arr) = 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) @@ -348,7 +351,7 @@         row (a,i) = newListArray (1, 2*n)                                  [ if j > n then (if j == i + n then 1 else 0)                                             else a ! j-                                 | j <- [1..(2*n)] ]+                                 | j <- [1..2*n] ]          mapM row (zip (elems arr) [1..]) >>= newListArray (1, n) @@ -369,7 +372,7 @@ read a i j = readArray a i >>= flip readArray j  -_inv :: (IArray a e, MArray (u s) e (ST s), Fractional e)+_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))@@ -377,43 +380,54 @@     let m = snd $ bounds mat         n = 2*m +        swap a i j = do+            tmp <- readArray a i+            readArray a j >>= writeArray a i+            writeArray a j tmp+     a <- augment mkArrayST mat      flip mapM_ [1..m] $ \k -> do-        flip mapM_ [(k+1)..m] $ \i -> do+        iPivot <- zip [k..m] <$> mapM (\i -> abs <$> read a i k) [k..m]+                    >>= 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+            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-                a_kk <- readArray a_k k-                writeArray a_i j (a_ij - a_kj * (a_ik / a_kk))+                writeArray a_i j (a_ij - a_kj * (a_ik / p))             writeArray a_i k 0 -    flip mapM_ [ m - k | k <- [1..(m-1)] ] $ \i -> do-        r1 <- readArray a i-        r2 <- readArray a (i+1)+    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) -        p <- readArray r2 (i+1) >>= return . (1 /)+        unless (i == m) $ do+            flip mapM_ [i+1..m] $ \k -> do+                a_k <- readArray a k+                p   <- readArray a_i k -        flip mapM_ [(i+1)..2*m] $ \j -> do-            c1 <- readArray r1 j-            c2 <- readArray r2 j-            writeArray r1 j (c2 * p + c1)+                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) -    result <- flip mapM [1..m] $ \i -> do-        r <- readArray a i-        p <- readArray r i-        -        mapM (\j -> (/ p) <$> readArray r j) [(m+1)..(2*m)]-            >>= return . listArray (1, m)-    -    return $ listArray (1, m) result+    mapM (\i -> readArray a i >>= getElems+                    >>= return . listArray (1, m) . drop m) [1..m]+        >>= return . listArray (1, m)  -_rank :: (IArray a e, MArray (u s) e (ST s),-           Fractional e, Eq e)+_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])       -> Array Int (a Int e)       -> ST s e@@ -484,30 +498,4 @@      liftM2 (*) (readSTRef pivotR) (readSTRef signR) --mat n = fromList [ [ j | j <- take n [i,(i^i)..] ] | i <- take n [1..] ]---{---- | The 'findIndex' function takes a predicate and a list and returns--- the index of the first element in the list satisfying the predicate,--- or 'Nothing' if there is no such element.-findIndex       :: (a -> Bool) -> [a] -> Maybe Int-findIndex p     = listToMaybe . findIndices p---- | The 'findIndices' function extends 'findIndex', by returning the--- indices of all elements satisfying the predicate, in ascending order.-findIndices      :: (a -> Bool) -> [a] -> [Int]--#if defined(USE_REPORT_PRELUDE) || !defined(__GLASGOW_HASKELL__)-findIndices p xs = [ i | (x,i) <- zip xs [0..], p x]-#else--- Efficient definition-findIndices p ls = loop 0# ls-                 where-                   loop _ [] = []-                   loop n (x:xs) | p x       = I# n : loop (n +# 1#) xs-                                 | otherwise = loop (n +# 1#) xs-#endif  /* USE_REPORT_PRELUDE */--}