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.4 +Version: 0.4.1 Synopsis: Efficient Matrix operations in 100% Haskell. Description: Efficient Matrix operations in 100% Haskell. . @@ -52,6 +52,11 @@ [@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. + . + [@v0.4.1@] The unsafe operations used in v0.4 turned out + to fatally fail on certain platforms. Revoked this + change. Thanks to Marcin Fatyga for pointing at the + issue. License: MIT diff --git a/quickcheck-tests.hs b/quickcheck-tests.hs --- a/quickcheck-tests.hs +++ b/quickcheck-tests.hs @@ -1,5 +1,8 @@ {-# LANGUAGE Haskell2010, TemplateHaskell #-} +import Prelude +import qualified Prelude as P + import Control.Monad import Numeric.Matrix @@ -9,7 +12,11 @@ import System.Exit +import Data.Ratio +import Data.Int +import Data.Word + dim :: Num a => a dim = 6 @@ -18,12 +25,40 @@ arbitrary = sequence (replicate dim (vector dim)) >>= return . fromList +basic m = m + m == scale m 2 && m * m * m == m ^ 3 + +prop_zero :: Word8 -> Bool +prop_zero = isZero . (zero :: Int -> Matrix Double) . fromIntegral + +prop_sum_trace :: Word8 -> Bool +prop_sum_trace n = ((== nd) . P.sum . trace . (unit :: Int -> Matrix Double) . fromIntegral) n' + where n' = n + 1 + nd = fromIntegral n' + +prop_basic_integer :: Matrix Integer -> Bool +prop_basic_integer = basic + +prop_basic_int :: Matrix Int -> Bool +prop_basic_int = basic + +prop_basic_float :: Matrix Float -> Bool +prop_basic_float = basic + +prop_basic_double :: Matrix Double -> Bool +prop_basic_double = basic + +prop_basic_rational :: Matrix Rational -> Bool +prop_basic_rational = basic + +prop_zero_det_double :: Matrix Double -> Bool +prop_zero_det_double m1 = let m1' = inv m1 in case m1' of (Just _) -> det m1 /= 0; _ -> det m1 == 0 + +prop_zero_det_rational :: Matrix Rational -> Bool +prop_zero_det_rational m1 = let m1' = inv m1 in case m1' of (Just _) -> det m1 /= 0; _ -> det m1 == 0 + 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 @@ -39,11 +74,8 @@ 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 @@ -71,7 +71,6 @@ import Data.Array.MArray import Data.Array.Unboxed import Data.Array.ST -import Data.Array.Base (unsafeRead, unsafeWrite) import Data.STRef import Data.Binary @@ -136,9 +135,16 @@ -- See @encode@ and @decode@. data family Matrix e + data instance Matrix Int = IntMatrix !Int !Int (Array Int (UArray Int Int)) +data instance Matrix Int32 + = Int32Matrix !Int !Int (Array Int (UArray Int Int32)) + +data instance Matrix Int64 + = Int64Matrix !Int !Int (Array Int (UArray Int Int64)) + data instance Matrix Float = FloatMatrix !Int !Int (Array Int (UArray Int Float)) @@ -189,8 +195,8 @@ = allWithIndex (\ix e -> m `at` ix == e) n | otherwise = False -instance (MatrixElement e) => NFData (Matrix e) where - rnf matrix = matrix `deepseq` () +instance (MatrixElement e, NFData e) => NFData (Matrix e) where + rnf matrix = map (\x -> x `deepseq` x) matrix `deepseq` () instance (MatrixElement e, Binary e) => Binary (Matrix e) where @@ -378,9 +384,23 @@ -- > dimensions empty == (0, 0) empty :: Matrix e + -- | Subtract a matrix from another. minus :: Matrix e -> Matrix e -> Matrix e + + -- | Add two matrices. + -- + -- You may also use the 'Num' instance for matrices, + -- i.e. @matrix1 + matrix2@ will work, too. plus :: Matrix e -> Matrix e -> Matrix e + + -- | Multiply two matrices /O(n^3)/. + -- + -- You may also use the 'Num' instance for matrices, + -- i.e. @matrix1 * matrix2@ will work, too. times :: Matrix e -> Matrix e -> Matrix e + + -- | Compute the inverse of a matrix. Returns @Nothing@ + -- if the matrix is not invertible. inv :: Matrix e -> Maybe (Matrix e) -- | Applies Bareiss multistep integer-preserving @@ -486,6 +506,30 @@ det (IntMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr) rank (IntMatrix _ _ arr) = runST (_rank thawsBoxed arr) +instance MatrixElement Int32 where + matrix d g = runST (_matrix Int32Matrix arrayST arraySTU d g) + fromList = _fromList Int32Matrix + + at (Int32Matrix _ _ arr) = _at arr + dimensions (Int32Matrix m n _) = (m, n) + row i (Int32Matrix _ _ arr) = _row i arr + col j (Int32Matrix _ _ arr) = _col j arr + toList (Int32Matrix _ _ arr) = _toList arr + det (Int32Matrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr) + rank (Int32Matrix _ _ arr) = runST (_rank thawsBoxed arr) + +instance MatrixElement Int64 where + matrix d g = runST (_matrix Int64Matrix arrayST arraySTU d g) + fromList = _fromList Int64Matrix + + at (Int64Matrix _ _ arr) = _at arr + dimensions (Int64Matrix m n _) = (m, n) + row i (Int64Matrix _ _ arr) = _row i arr + col j (Int64Matrix _ _ arr) = _col j arr + toList (Int64Matrix _ _ arr) = _toList arr + det (Int64Matrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr) + rank (Int64Matrix _ _ arr) = runST (_rank thawsBoxed arr) + instance MatrixElement Integer where matrix d g = runST (_matrix IntegerMatrix arrayST arrayST d g) fromList = _fromList IntegerMatrix @@ -530,7 +574,6 @@ 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 @@ -632,7 +675,7 @@ read :: (MArray a1 b m, MArray a (a1 Int b) m) => a Int (a1 Int b) -> Int -> Int -> m b -read a i j = unsafeRead a i >>= flip unsafeRead j +read a i j = readArray a i >>= flip readArray j _inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Ord e, Show e) @@ -644,9 +687,9 @@ n = 2*m swap a i j = do - tmp <- unsafeRead a i - unsafeRead a j >>= unsafeWrite a i - unsafeWrite a j tmp + tmp <- readArray a i + readArray a j >>= writeArray a i + writeArray a j tmp okay <- newSTRef True @@ -662,37 +705,37 @@ swap a iPivot k forM_ [k+1..m] $ \i -> do - a_i <- unsafeRead a i - a_k <- unsafeRead a k + a_i <- readArray a i + a_k <- readArray 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 + 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 invertible <- readSTRef okay if invertible then do forM_ [ m - v | v <- [0..m-1] ] $ \i -> do - a_i <- unsafeRead a i - p <- unsafeRead a_i i - unsafeWrite a_i i 1 + a_i <- readArray a i + p <- readArray a_i i + writeArray a_i i 1 forM_ [i+1..n] $ \j -> do - unsafeRead a_i j >>= unsafeWrite a_i j . (/ p) + readArray a_i j >>= writeArray a_i j . (/ p) unless (i == m) $ do forM_ [i+1..m] $ \k -> do - a_k <- unsafeRead a k - p <- unsafeRead a_i k + a_k <- readArray a k + p <- readArray a_i k 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) + a_ij <- readArray a_i j + a_kj <- readArray a_k j + writeArray a_i j (a_ij - p * a_kj) - mapM (\i -> unsafeRead a i >>= getElems + mapM (\i -> readArray a i >>= getElems >>= return . listArray (1, m) . drop m) [1..m] >>= return . Just . listArray (1, m) @@ -707,9 +750,9 @@ n = snd $ bounds (mat ! 1) swap a i j = do - tmp <- unsafeRead a i - unsafeRead a j >>= unsafeWrite a i - unsafeWrite a j tmp + tmp <- readArray a i + readArray a j >>= writeArray a i + writeArray a j tmp a <- thaws mat >>= arrays @@ -726,17 +769,17 @@ let ix = fromJust switchRow + pivotRow when (pivotRow /= ix) (swap a pivotRow ix) - a_p <- unsafeRead a k - pivot <- unsafeRead a_p k + a_p <- readArray a k + pivot <- readArray a_p k prev <- readSTRef prevR forM_ [pivotRow+1..m] $ \i -> do - a_i <- unsafeRead a i + a_i <- readArray 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) + 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) `divide` prev) writeSTRef ixPivot (pivotRow + 1) @@ -774,9 +817,9 @@ when (not $ null sf) $ do let sw = head sf - row <- unsafeRead a sw - unsafeRead a k >>= unsafeWrite a sw - unsafeWrite a k row + row <- readArray a sw + readArray a k >>= writeArray a sw + writeArray a k row read a k k >>= writeSTRef pivotR readSTRef signR >>= writeSTRef signR . negate @@ -787,21 +830,23 @@ unless (sign' == 0) $ do pivot' <- readSTRef pivotR forM_ [(k+1)..size] $ \i -> do - a_i <- unsafeRead a i + a_i <- readArray a i forM [(k+1)..size] $ \j -> do - a_ij <- unsafeRead a_i j - a_ik <- unsafeRead a_i k + a_ij <- readArray a_i j + a_ik <- readArray a_i k a_kj <- read a k j - unsafeWrite a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev) + writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev) liftM2 (*) (readSTRef pivotR) (readSTRef signR) _mult :: MatrixElement e => Matrix e -> Matrix e -> Matrix e -_mult a b = let rowsA = numRows a - rowsB = numRows b - colsB = numCols b - in matrix (rowsA, colsB) (\(i,j) -> L.foldl' (+) 0 [a `at` (i, k) * b `at` (k, j) | k <- [1..rowsB]]) +_mult a b = + let rowsA = numRows a + rowsB = numRows b + colsB = numCols b + gen i j = L.foldl' (+) 0 [a `at` (i, k) * b `at` (k, j) | k <- [1..rowsB]] + in matrix (rowsA, colsB) (uncurry gen) _matrix :: (IArray a1 (u Int e), IArray u e, @@ -818,8 +863,8 @@ forM_ [1..m] $ \i -> do cols <- newArrayU (1, n) 0 forM_ [1..n] $ \j -> do - unsafeWrite cols j (g (i,j)) - U.unsafeFreeze cols >>= unsafeWrite rows i + writeArray cols j (g (i,j)) + U.unsafeFreeze cols >>= writeArray rows i U.unsafeFreeze rows >>= return . c m n