diff --git a/bed-and-breakfast.cabal b/bed-and-breakfast.cabal --- a/bed-and-breakfast.cabal +++ b/bed-and-breakfast.cabal @@ -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 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 diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -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])