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.3.1 +Version: 0.3.2 Synopsis: Efficient Matrix operations in 100% Haskell. Description: Efficient Matrix operations in 100% Haskell. . @@ -40,11 +40,14 @@ . [@v0.3@] Added a QuickCheck test suite, fixed a bug in @det@ (det would crash for singular matrices, where it should - return 0; it does so no). + return 0). . [@v0.3.1@] Added TemplateHaskell syntactic sugar (see "Numeric.Matrix.Sugar"). Rewrote multiplication. @matrix@ function build an array faster now. + . + [@v0.3.2@] Numeric.Matrix.Sugar was not mentioned in the + cabal file. Improved test suite. Improved documentation. License: MIT @@ -62,16 +65,19 @@ location: hub.darcs.net:bed-and-breakfast Library - Exposed-Modules: Numeric.Matrix + Exposed-Modules: Numeric.Matrix, + Numeric.Matrix.Sugar Build-Depends: base >= 4.5 && < 5, deepseq >= 1.3, - array >= 0.4 + array >= 0.4, + template-haskell >= 2.7 Hs-Source-Dirs: src Test-Suite quickcheck type: exitcode-stdio-1.0 - main-is: quick-check-tests.hs + main-is: quickcheck-tests.hs + GHC-Options: -O2 build-depends: base >= 4.5 && < 5, - bed-and-breakfast == 0.3.1, + bed-and-breakfast == 0.3.2, QuickCheck >= 2.4.2 diff --git a/quick-check-tests.hs b/quick-check-tests.hs deleted file mode 100644 --- a/quick-check-tests.hs +++ /dev/null @@ -1,36 +0,0 @@ -{-# LANGUAGE Haskell2010, TemplateHaskell #-} - -import Control.Monad - -import Numeric.Matrix - -import Test.QuickCheck -import Test.QuickCheck.All - -import System.Exit - - -dim :: Num a => a -dim = 5 - -instance (MatrixElement e, Arbitrary e) => Arbitrary (Matrix e) where - - arbitrary = sequence (replicate dim (vector dim)) >>= return . fromList - - -prop_det :: Matrix Rational -> Bool -prop_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 - -prop_rank :: Matrix Rational -> Bool -prop_rank m1 = let r = rank m1 in if r < dim then det m1 == 0 else r == dim - - -main = do - success <- $(quickCheckAll) - - (if success then exitSuccess else exitFailure) - - diff --git a/quickcheck-tests.hs b/quickcheck-tests.hs new file mode 100644 --- /dev/null +++ b/quickcheck-tests.hs @@ -0,0 +1,49 @@ +{-# LANGUAGE Haskell2010, TemplateHaskell #-} + +import Control.Monad + +import Numeric.Matrix + +import Test.QuickCheck +import Test.QuickCheck.All + +import System.Exit + + +dim :: Num a => a +dim = 6 + +instance (MatrixElement e, Arbitrary e) => Arbitrary (Matrix e) where + + arbitrary = sequence (replicate dim (vector dim)) >>= return . fromList + + +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 + +prop_rank :: Matrix Rational -> Bool +prop_rank m1 = let r = rank m1 in if r < dim then det m1 == 0 else r == dim + +prop_trace_select :: Matrix Rational -> Bool +prop_trace_select m = trace m == select (uncurry (==)) m + +prop_transpose_twice :: Matrix Rational -> Bool +prop_transpose_twice m = transpose (transpose m) == m + +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) + + diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -59,6 +59,10 @@ 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 @@ -71,9 +75,12 @@ import Data.Array.MArray import Data.Array.Unboxed import Data.Array.ST -import qualified Data.Array.Unsafe as U +--import Data.Array.IO import Data.STRef +--import System.IO.Unsafe +import qualified Data.Array.Unsafe as U + import Data.Typeable import Prelude hiding (any, all, read, map) @@ -276,7 +283,14 @@ class (Eq e, Num e) => MatrixElement e where + -- | Creates a matrix of the given size using a generator + -- function for the value of each component. matrix :: (Int, Int) -> ((Int, Int) -> e) -> Matrix e + + -- | Builds a list from a matrix for the indices for which + -- the given predicate matches. + -- + -- > trace == select (uncurry (==)) select :: ((Int, Int) -> Bool) -> Matrix e -> [e] -- | Returns the component at the given position in the matrix. @@ -313,7 +327,7 @@ -- | An identity square matrix of the given size. -- - -- >> unit 4 + -- >>> unit 4 -- 1 0 0 0 -- 0 1 0 0 -- 0 0 1 0 @@ -322,7 +336,7 @@ -- | A square matrix of the given size consisting of all zeros. -- - -- >> zero 2 + -- >>> zero 2 -- 0 0 -- 0 0 zero :: Int -> Matrix e @@ -330,7 +344,7 @@ -- | A square matrix which trace is the given list, all other components -- set to zero. -- - -- >> diag [1,4,7,9] + -- >>> diag [1,4,7,9] -- 1 0 0 0 -- 0 4 0 0 -- 0 0 7 0 @@ -492,6 +506,7 @@ 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 @@ -587,6 +602,15 @@ => (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 @@ -796,3 +820,35 @@ U.unsafeFreeze cols >>= writeArray 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 + +"det/pow" + forall a k. det (a ^ k) = (det a) ^ k + + #-} diff --git a/src/Numeric/Matrix/Sugar.hs b/src/Numeric/Matrix/Sugar.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Matrix/Sugar.hs @@ -0,0 +1,31 @@ +{-# LANGUAGE Haskell2010, TemplateHaskell #-} + +module Numeric.Matrix.Sugar ( + iMatrix, dMatrix + ) where + +import Numeric.Matrix +import Language.Haskell.TH +import Language.Haskell.TH.Quote + +import Data.Data +import Data.Complex +import Data.Ratio + +type ReadM a = String -> Matrix a + +iMatrix = QuasiQuoter + (quoter (read :: ReadM Integer)) + undefined undefined undefined + +dMatrix = QuasiQuoter + (quoter (read :: ReadM Double)) + undefined undefined undefined + +quoter :: (Data e, MatrixElement e) + => (String -> Matrix e) -> String -> Q Exp + +quoter read str = do + let qExp = dataToExpQ (const Nothing) (toList (read str)) + [e| fromList $qExp |] +