bed-and-breakfast (empty) → 0.1
raw patch · 4 files changed
+556/−0 lines, 4 filesdep +arraydep +basesetup-changed
Dependencies added: array, base
Files
- LICENSE +18/−0
- Setup.hs +4/−0
- bed-and-breakfast.cabal +21/−0
- src/Numeric/Matrix.hs +513/−0
+ LICENSE view
@@ -0,0 +1,18 @@+Permission is hereby granted, free of charge, to any person obtaining+a copy of this software and associated documentation files (the "Software"),+to deal in the Software without restriction, including without limitation+the rights to use, copy, modify, merge, publish, distribute, sublicense,+and/or sell copies of the Software, and to permit persons to whom the+Software is furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included+in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS+OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL+THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING+FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER+DEALINGS IN THE SOFTWARE.+
+ Setup.hs view
@@ -0,0 +1,4 @@+import Distribution.Simple++main = defaultMain+
+ bed-and-breakfast.cabal view
@@ -0,0 +1,21 @@+Name: bed-and-breakfast+Version: 0.1+Synopsis: Efficient Matrix operations in 100% Haskell.+Description: Efficient Matrix operations in 100% Haskell.+ +License: MIT+License-File: LICENSE+Author: Julian Fleischer <julian.fleischer@fu-berlin.de>+Maintainer: Julian Fleischer <julian.fleischer@fu-berlin.de>+Build-Type: Simple+Cabal-Version: >= 1.8+Category: Data+Stability: stable++Library+ Exposed-Modules: Numeric.Matrix+ Build-Depends: base >= 4.5 && < 5,+ array >= 0.4+ Hs-Source-Dirs: src++
+ src/Numeric/Matrix.hs view
@@ -0,0 +1,513 @@+{-# LANGUAGE Haskell2010+ , TypeFamilies+ , FlexibleContexts+ , Trustworthy+ #-}+{-# OPTIONS -Wall -fno-warn-name-shadowing #-}++module Numeric.Matrix (+ Matrix,+ MatrixElement (+ matrix,+ fromList,++ unit,+ zero,+ diag,+ empty,++ at,+ row,+ col,++ select,+ toList,++ dimensions,+ numRows,+ numCols,++ isUnit,+ isZero,+ isDiagonal,+ isEmpty,+ isSquare,+ + det,+ rank,+ transpose,+ trace,++ minus,+ plus,+ times,+ inv,++ map,+ all,+ any,++ mapWithIndex,+ allWithIndex,+ anyWithIndex+ )+) where+++import Control.Applicative+import Control.Monad+import Control.Monad.ST++import Data.Ratio+import Data.Complex+import qualified Data.List as L+import Data.Array.IArray+import Data.Array.MArray+import Data.Array.Unboxed+import Data.Array.ST+import Data.STRef++import Prelude hiding (any, all, read)+import qualified Prelude as P+++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 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 (MatrixElement e, Show e) => Show (Matrix e) where+ show = unlines . P.map showRow . toList+ where+ showRow = unwords . P.map ((' ':) . show)+++class Division e where+ divide :: e -> e -> e++instance Division Int where divide = quot+-- instance Division Int8 where divide = quot+-- instance Division Int16 where divide = quot+-- instance Division Int32 where divide = quot+-- instance Division Int64 where divide = quot+instance Division Integer where divide = quot+instance Division Float where divide = (/)+instance Division Double where divide = (/)+instance Integral a => Division (Ratio a) where divide = (/)+instance RealFloat a => Division (Complex a) where divide = (/)++++class (Eq e, Num e) => MatrixElement e where++ matrix :: (Int, Int) -> ((Int, Int) -> e) -> Matrix e+ select :: ((Int, Int) -> Bool) -> Matrix e -> [e]+ at :: Matrix e -> (Int, Int) -> e++ row :: Int -> Matrix e -> [e]+ col :: Int -> Matrix e -> [e]++ dimensions :: Matrix e -> (Int, Int)+ numRows :: Matrix e -> Int+ numCols :: Matrix e -> Int++ fromList :: [[e]] -> Matrix e+ toList :: Matrix e -> [[e]]++ unit :: Int -> Matrix e+ zero :: Int -> Matrix e+ diag :: [e] -> Matrix e+ empty :: Matrix e++ minus :: Matrix e -> Matrix e -> Matrix e+ plus :: Matrix e -> Matrix e -> Matrix e+ times :: Matrix e -> Matrix e -> Matrix e+ inv :: Matrix e -> Maybe (Matrix e)++-- adjugate :: Matrix e -> Matrix e+-- cofactors :: Matrix e -> Matrix e ; cofactors = undefined+ det :: Matrix e -> e+ transpose :: Matrix e -> Matrix e+ rank :: Matrix e -> e+ trace :: Matrix e -> [e]++ isUnit :: Matrix e -> Bool+ isDiagonal :: Matrix e -> Bool+ isZero :: Matrix e -> Bool+ isEmpty :: Matrix e -> Bool+ isSquare :: Matrix e -> Bool++ select p m = [ at m (i,j) | i <- [1..numRows m]+ , j <- [1..numCols m]+ , p (i,j) ]++ at m (i, j) = ((!! j) . (!! i) . toList) m++ map :: MatrixElement f => (e -> f) -> Matrix e -> Matrix f+ all :: (e -> Bool) -> Matrix e -> Bool+ any :: (e -> Bool) -> Matrix e -> Bool++ mapWithIndex :: MatrixElement f => ((Int, Int) -> e -> f) -> Matrix e -> Matrix f+ allWithIndex :: ((Int, Int) -> e -> Bool) -> Matrix e -> Bool+ anyWithIndex :: ((Int, Int) -> e -> Bool) -> Matrix e -> Bool++ unit n = fromList [[ if i == j then 1 else 0 | j <- [1..n]] | i <- [1..n] ]+ zero n = matrix (n,n) (const 0)+ empty = fromList []+ diag xs = matrix (n,n) (\(i,j) -> if i == j then xs !! (i-1) else 0)+ where n = length xs+ + row i m = ((!! (i-1)) . toList) m+ col i m = (row i . transpose) m++ numRows = fst . dimensions+ numCols = snd . dimensions+ dimensions m = case toList m of [] -> (0, 0)+ (x:xs) -> (length xs + 1, length x)++-- adjugate = transpose . cofactors+ transpose = fromList . L.transpose . toList+ trace = select (uncurry (==))+ inv _ = Nothing++ isZero = all (== 0)+ isUnit m = isSquare m && P.all (== 1) (trace m)+ isEmpty m = numRows m == 0 || numCols m == 0+ isDiagonal = allWithIndex (uncurry $ \x y z -> if x /= y then z == 0 else True)+ isSquare m = let (a, b) = dimensions m in a == b++ map f = mapWithIndex (const f)+ all f = allWithIndex (const f)+ any f = anyWithIndex (const f)++ mapWithIndex f m = matrix (dimensions m) (\x -> f x (m `at` x))+ allWithIndex f m = P.all id [ f (i, j) (m `at` (i,j))+ | i <- [1..numRows m], j <- [1..numCols m]]+ anyWithIndex f m = P.any id [ f (i, j) (m `at` (i,j))+ | i <- [1..numRows m], j <- [1..numCols m]]++ a `plus` b+ | dimensions a /= dimensions b = error "Matrix.plus: dimensions don't match."+ | otherwise = matrix (dimensions a) (\x -> a `at` x + b `at` x)+ a `minus` b+ | dimensions a /= dimensions b = error "Matrix.minus: dimensions don't match."+ | otherwise = matrix (dimensions a) (\x -> a `at` x - b `at` x)+ a `times` b+ | numRows a /= numCols b = error "Matrix.times: `numRows a' and `numCols b' don't match."+ | otherwise = fromList [ [ row i a `dotProd` col j b+ | j <- [1..numCols b] ]+ | i <- [1..numRows a] ]+++instance MatrixElement Int where+ matrix = _matrix IntMatrix+ fromList = _fromList IntMatrix++ at (IntMatrix _ _ arr) = _at arr+ dimensions (IntMatrix m n _) = (m, n)+ 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)+ det (IntMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)+ rank = undefined -- runST (_rank thawsBoxed arr)++instance MatrixElement Integer where+ matrix = _matrix IntegerMatrix+ fromList = _fromList IntegerMatrix++ at (IntegerMatrix _ _ arr) = _at arr+ dimensions (IntegerMatrix m n _) = (m, n)+ 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)+ det (IntegerMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr)+ rank = undefined -- runST (_rank thawsBoxed arr)++instance MatrixElement Float where+ matrix = _matrix FloatMatrix+ fromList = _fromList FloatMatrix++ at (FloatMatrix _ _ arr) = _at arr+ dimensions (FloatMatrix m n _ ) = (m, n)+ row i (FloatMatrix _ _ arr) = _row i arr+ 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)+ det (FloatMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)+ rank (FloatMatrix _ _ arr) = runST (_rank thawsBoxed arr)++instance MatrixElement Double where+ matrix = _matrix DoubleMatrix+ fromList = _fromList DoubleMatrix++ at (DoubleMatrix _ _ arr) = _at arr+ dimensions (DoubleMatrix m n _ ) = (m, n)+ row i (DoubleMatrix _ _ arr) = _row i arr+ 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)+ 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+ matrix = _matrix RatioMatrix+ fromList = _fromList RatioMatrix++ at (RatioMatrix _ _ arr) = _at arr+ dimensions (RatioMatrix m n _ ) = (m, n)+ row i (RatioMatrix _ _ arr) = _row i arr+ 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)+ 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+ matrix = _matrix ComplexMatrix+ fromList = _fromList ComplexMatrix++ at (ComplexMatrix _ _ arr) = _at arr+ dimensions (ComplexMatrix m n _ ) = (m, n)+ 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)+ det (ComplexMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr)+ rank (ComplexMatrix _ _ arr) = runST (_rank thawsBoxed arr)+++_at :: (IArray a (u Int e), IArray u e)+ => a Int (u Int e) -> (Int, Int) -> e+_at arr (i,j) = arr ! i ! j++_row, _col :: (IArray a (u Int e), IArray u e) => Int -> a Int (u Int e) -> [e]+_row i arr = let row = arr ! i in [ row ! j | j <- [1..(snd (bounds arr))] ]+_col j arr = [ arr ! i ! j | i <- [1..(snd (bounds arr))] ]++_matrix :: (IArray a (u Int e), IArray u e)+ => (Int -> Int -> (a Int (u Int e)) -> matrix e)+ -> (Int, Int)+ -> ((Int, Int) -> e)+ -> matrix e+_matrix c (numRows, numCols) generator =+ c numRows numCols+ $ array (1, numRows)+ $ [ (i, array (1, numCols) [(j, generator (i, j))+ | j <- [1..numCols]])+ | i <- [1..numRows] ]++_toList :: (IArray a e) => Array Int (a Int e) -> [[e]]+_toList = P.map elems . elems++_fromList :: (IArray a (u Int e), IArray u e)+ => (Int -> Int -> a Int (u Int e) -> matrix e) -> [[e]] -> matrix e+_fromList c xs =+ let lengths = P.map length xs+ numCols = foldl1 min lengths+ numRows = length lengths+ + in c numRows numCols+ $ array (1, numRows)+ $ zip [1..numRows]+ $ P.map (array (1, numCols) . zip [1..numCols]) xs++dotProd :: Num a => [a] -> [a] -> a+dotProd x = L.foldl' (+) 0 . zipWith (*) x++thawsBoxed :: (IArray a e, MArray (STArray s) e (ST s))+ => Array Int (a Int e)+ -> ST s [STArray s Int e]+thawsBoxed = mapM thaw . elems++thawsUnboxed :: (IArray a e, MArray (STUArray s) e (ST s))+ => Array Int (a Int e)+ -> ST s [STUArray s Int e]+thawsUnboxed = mapM thaw . elems++arrays :: [(u s) Int e]+ -> ST s ((STArray s) Int ((u s) Int e))+arrays list = newListArray (1, length list) list++augment :: (IArray a e, MArray (u s) e (ST s), Num e)+ => ((Int, Int) -> [e] -> ST s ((u s) Int e))+ -> Array Int (a Int e)+ -> ST s (STArray s Int (u s Int e))+augment _ arr = do+ let (_, n) = bounds arr+ 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)] ]+ + mapM row (zip (elems arr) [1..]) >>= newListArray (1, n)++boxedST :: MArray (STArray s) e (ST s)+ => (Int, Int) -> [e] -> ST s ((STArray s) Int e)+boxedST = newListArray++unboxedST :: MArray (STUArray s) e (ST s)+ => (Int, Int) -> [e] -> ST s ((STUArray s) Int e)+unboxedST = newListArray+++tee :: Monad m => (b -> m a) -> b -> m b+tee f x = f x >> return x++read :: (MArray a1 b m, MArray a (a1 Int b) m) =>+ a Int (a1 Int b) -> Int -> Int -> m b+read a i j = readArray a i >>= flip readArray j+++_inv :: (IArray a e, MArray (u s) e (ST s), Fractional e)+ => ((Int, Int) -> [e] -> ST s ((u s) Int e))+ -> Array Int (a Int e)+ -> ST s (Array Int (a Int e))+_inv mkArrayST mat = do+ let m = snd $ bounds mat+ n = 2*m++ a <- augment mkArrayST mat++ flip mapM_ [1..m] $ \k -> do+ 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+ a_kk <- readArray a_k k+ writeArray a_i j (a_ij - a_kj * (a_ik / a_kk))+ writeArray a_i k 0++ flip mapM_ [ m - k | k <- [1..(m-1)] ] $ \i -> do+ r1 <- readArray a i+ r2 <- readArray a (i+1)++ p <- readArray r2 (i+1) >>= return . (1 /)++ flip mapM_ [(i+1)..2*m] $ \j -> do+ c1 <- readArray r1 j+ c2 <- readArray r2 j+ writeArray r1 j (c2 * p + c1)++ 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+++_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+_rank thaws mat = do+ let m = snd $ bounds mat+ n = snd $ bounds (mat ! 1)++ a <- thaws mat >>= arrays++ trace <- flip mapM [1..m] $ \k -> do+ 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+ a_kk <- readArray a_k k+ writeArray a_i j (a_ij - a_kj * (a_ik / a_kk))+ writeArray a_i k 0+ read a k k++ return $ fromIntegral $ length $ filter (/= 0) trace+++_det :: (IArray a e, MArray (u s) e (ST s),+ Num e, Eq e, Division e)+ => (Array Int (a Int e) -> ST s [(u s) Int e])+ -> Array Int (a Int e) -> ST s e+_det thaws mat = do++ let size = snd $ bounds mat++ a <- thaws mat >>= arrays++ signR <- newSTRef 1+ pivotR <- newSTRef 1++ flip mapM_ [1..size] $ \k -> do++ prev <- readSTRef pivotR+ pivot <- read a k k >>= tee (writeSTRef pivotR)++ when (pivot == 0) $ do+ s <- flip mapM [(k+1)..size] $ \r -> do+ a_rk <- read a r k+ if a_rk == 0 then return 0 else return r+ let sf = filter (>0) s++ when (not $ null sf) $ do+ let sw = head sf++ 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++ pivot' <- readSTRef pivotR+ flip mapM [(k+1)..size] $ \i -> do+ a_i <- readArray a i+ flip mapM [(k+1)..size] $ \j -> do+ a_ij <- readArray a_i j+ a_ik <- readArray a_i k+ a_kj <- read a k j+ writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)++ 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 */+-}+