matrix-sized 0.0.4 → 0.1.0
raw patch · 11 files changed
+367/−35 lines, 11 filesdep +bytestringdep +bytestring-lexingdep +conduit-extradep ~vector
Dependencies added: bytestring, bytestring-lexing, conduit-extra, store
Dependency ranges changed: vector
Files
- matrix-sized.cabal +9/−4
- src/Data/Matrix/Static/Dense.hs +29/−0
- src/Data/Matrix/Static/Dense/Mutable.hs +9/−0
- src/Data/Matrix/Static/Generic/Mutable.hs +16/−1
- src/Data/Matrix/Static/IO.hs +125/−0
- src/Data/Matrix/Static/LinearAlgebra.hs +33/−4
- src/Data/Matrix/Static/LinearAlgebra/Internal.hs +0/−1
- src/Data/Matrix/Static/Sparse.hs +90/−10
- tests/Test/Base.hs +33/−5
- tests/Test/LinearAlgebra.hs +21/−9
- tests/Test/Utils.hs +2/−1
matrix-sized.cabal view
@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: a94165a01f28a1a01026b95b46341900f92789d9a5b17aec1fe9ab414b560a26+-- hash: 2b1169a59b262a2ee39daffd72cfe59b2b4f716e82ef7df496b186c155e1b003 name: matrix-sized-version: 0.0.4+version: 0.1.0 synopsis: Haskell matrix library with interface to C++ linear algebra libraries. description: A Haskell implementation of matrices with statically known sizes. The library also comes with the bindings to high performance C++ linear algebra libraries such as Eigen and Spectra. category: Math@@ -390,6 +390,7 @@ Data.Matrix.Static.LinearAlgebra.Types Data.Matrix.Static.Generic Data.Matrix.Static.Generic.Mutable+ Data.Matrix.Static.IO other-modules: Data.Matrix.Static.Internal Data.Matrix.Static.LinearAlgebra.Internal@@ -409,9 +410,13 @@ stdc++ build-depends: base >=4.10 && <5+ , bytestring+ , bytestring-lexing , conduit+ , conduit-extra , primitive >=0.6.4.0 , singletons+ , store , vector >=0.11 if flag(parallel) cc-options: -fopenmp@@ -435,9 +440,9 @@ , data-ordlist , ieee754 , matrix-sized- , primitive >=0.6.4.0 , singletons+ , store , tasty , tasty-quickcheck- , vector >=0.11+ , vector default-language: Haskell2010
src/Data/Matrix/Static/Dense.hs view
@@ -105,6 +105,8 @@ import Data.Tuple (swap) import qualified Data.List as L import Text.Printf (printf)+import Data.Store (Store(..), Size(..), decodeExWith)+import Foreign.Storable (sizeOf) import Data.Matrix.Static.Dense.Mutable (MMatrix (..)) import qualified Data.Matrix.Static.Dense.Mutable as DM@@ -116,6 +118,25 @@ data Matrix :: C.MatrixKind where Matrix :: (SingI r, SingI c) => v a -> Matrix r c v a +instance (G.Vector v a, Store (v a), SingI r, SingI c) =>+ Store (Matrix r c v a) where+ size = VarSize $ \(Matrix vec) -> case size of+ VarSize f ->+ 2 * sizeOf (0 :: Int) + f vec++ poke mat@(Matrix vec) = poke r >> poke c >> poke vec+ where+ (r,c) = C.dim mat+ peek = do+ r' <- peek+ c' <- peek+ if r' /= r || c' /= c+ then error $ "Dimensions donot match: " <> show (r,c) <> " /= " <> show (r',c')+ else Matrix <$> peek+ where+ r = fromIntegral $ fromSing (sing :: Sing r) :: Int+ c = fromIntegral $ fromSing (sing :: Sing c) :: Int+ instance (G.Vector v a, Show a) => Show (Matrix r c v a) where show mat = printf "(%d x %d)\n%s" r c vals where@@ -169,6 +190,14 @@ -- | Create a vector by concatenating columns. flatten (Matrix vec) = vec {-# INLINE flatten #-}++ transpose mat@(Matrix vec)+ | r == 1 || c == 1 = Matrix vec+ | otherwise = Matrix $ G.generate (r*c) $ \x ->+ C.unsafeIndex mat $ x `divMod` c+ where+ (r, c) = C.dim mat+ {-# INLINE transpose #-} thaw (Matrix v) = MMatrix <$> G.thaw v {-# INLINE thaw #-}
src/Data/Matrix/Static/Dense/Mutable.hs view
@@ -49,6 +49,15 @@ idx = i + j * r {-# INLINE unsafeWrite #-} + unsafeModify mat@(MMatrix v) f (i,j) = GM.unsafeModify v f idx+ where + (r, _) = C.dim mat+ idx = i + j * r+ {-# INLINE unsafeModify #-}++ fill (MMatrix v) x = GM.set v x+ {-# INLINE fill #-}+ new :: forall r c s. (SingI r, SingI c, PrimMonad s) => s (MMatrix r c v (PrimState s) a) new = MMatrix <$> GM.new (r*c)
src/Data/Matrix/Static/Generic/Mutable.hs view
@@ -5,6 +5,7 @@ module Data.Matrix.Static.Generic.Mutable ( MMatrix(..) , MMatrixKind+ , fillDiagonal ) where import Control.Monad.Primitive (PrimMonad, PrimState)@@ -23,12 +24,26 @@ unsafeWrite :: PrimMonad s => mat r c v (PrimState s) a -> (Int, Int) -> a -> s () + unsafeModify :: PrimMonad s+ => mat r c v (PrimState s) a+ -> (a -> a)+ -> (Int, Int) -> s ()++ fill :: PrimMonad s => mat r c v (PrimState s) a -> a -> s ()+ -- | Create a mutable matrix without initialization new :: (SingI r, SingI c, PrimMonad s) => s (mat r c v (PrimState s) a) replicate :: (SingI r, SingI c, PrimMonad s) => a -> s (mat r c v (PrimState s) a) - {-# MINIMAL dim, unsafeRead, unsafeWrite, new, replicate #-}+fillDiagonal :: (PrimMonad s, MMatrix mat v a)+ => mat r c v (PrimState s) a -> a -> s ()+fillDiagonal mat x = mapM_ (\i -> unsafeWrite mat (i,i) x) [0..n]+ where+ n = min r c - 1+ (r,c) = dim mat+{-# INLINE fillDiagonal #-}+ {- write :: (PrimMonad s, MMatrix m v a)
+ src/Data/Matrix/Static/IO.hs view
@@ -0,0 +1,125 @@+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE StrictData #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}++module Data.Matrix.Static.IO+ ( fromMM+ , withMM+ ) where++import qualified Data.ByteString.Char8 as B+import Conduit+import Control.Monad (when)+import qualified Data.Vector.Generic as G+import Data.ByteString.Lex.Fractional (readExponential, readSigned)+import Data.ByteString.Lex.Integral (readDecimal, readDecimal_)+import Data.Singletons+import Data.Maybe+import Data.Singletons.TypeLits++import qualified Data.Matrix.Static.Sparse as S++data MMElem = MMReal+ | MMComplex+ | MMInteger+ | MMPattern+ deriving (Eq)++class IOElement a where+ decodeElem :: B.ByteString -> a+ elemType :: Proxy a -> MMElem++instance IOElement Int where+ decodeElem x = fst . fromMaybe errMsg . readSigned readDecimal $ x+ where+ errMsg = error $ "readInt: Fail to cast ByteString to Int:" ++ show x+ elemType _ = MMInteger++instance IOElement Double where+ decodeElem x = fst . fromMaybe errMsg . readSigned readExponential $ x+ where+ errMsg = error $ "readDouble: Fail to cast ByteString to Double:" ++ show x+ elemType _ = MMReal++withMM :: forall m v a b. (Monad m, G.Vector v a, IOElement a)+ => ConduitT () B.ByteString m ()+ -> (forall r c. S.SparseMatrix r c v a -> m b)+ -> m b+withMM conduit f = do+ (ty, (r,c,nnz)) <- parseHeader conduit+ when (elemType (Proxy :: Proxy a) /= ty) $ error "Element types do not match"+ withSomeSing (fromIntegral (r :: Int)) $ \(SNat :: Sing r) ->+ withSomeSing (fromIntegral (c :: Int)) $ \(SNat :: Sing c) -> do+ mat@(S.SparseMatrix v _ _) <- S.fromTripletC triplets+ let n = G.length v+ when (n /= nnz) $ error $+ "number of non-zeros do not match: " <> show nnz <> "/=" <> show n+ f (mat :: S.SparseMatrix r c v a)+ where+ triplets = conduit .| linesUnboundedAsciiC .| filterC (not . (=='%') . B.head) .|+ (dropC 1 >> streamTriplet)++fromMM :: forall m r c v a. (Monad m, SingI r, SingI c, G.Vector v a, IOElement a)+ => ConduitT () B.ByteString m () -> m (S.SparseMatrix r c v a)+fromMM conduit = do+ (ty, (r,c,nnz)) <- parseHeader conduit+ mat@(S.SparseMatrix v _ _) <- case () of+ _ | elemType (Proxy :: Proxy a) /= ty -> error "Element types do not match"+ | (r, c) /= (nrow, ncol) -> error $ "Dimensions do not match: " <>+ show (r,c) <> "/=" <> show (nrow,ncol)+ | otherwise -> S.fromTripletC triplets+ let n = G.length v+ if n /= nnz+ then error $ "number of non-zeros do not match: " <> show nnz <> "/=" <> show n+ else return mat+ where+ triplets = conduit .| linesUnboundedAsciiC .| filterC (not . (=='%') . B.head) .|+ (dropC 1 >> streamTriplet)+ nrow = fromIntegral $ fromSing (sing :: Sing r) :: Int+ ncol = fromIntegral $ fromSing (sing :: Sing c) :: Int++parseHeader :: Monad m => ConduitT () B.ByteString m () -> m (MMElem, (Int, Int, Int))+parseHeader conduit = runConduit $+ conduit .| linesUnboundedAsciiC .| headerParser+ where+ parse x+ | "%%MatrixMarket" `B.isPrefixOf` x = case B.words x of+ [_, _, format, ty, form] -> + let ty' = case ty of+ "real" -> MMReal+ "complex" -> MMComplex+ "integer" -> MMInteger+ "pattern" -> MMPattern+ t -> error $ "Unknown type: " <> show t+ in ty'+ _ -> error $ "Cannot parse header: " <> show x+ | otherwise = error $ "Cannot parse header: " <> show x+ headerParser = do+ ty <- headC >>= \case+ Nothing -> error "Empty file"+ Just header -> return $ parse header+ line <- filterC (not . (=='%') . B.head) .| headC+ case line of+ Nothing -> error "Empty file"+ Just x ->+ let [r, c, nnz] = map decodeElem $ B.words x+ in return (ty, (r, c, nnz))+{-# INLINE parseHeader #-}++streamTriplet :: (Monad m, IOElement a) => ConduitT B.ByteString (Int, Int, a) m ()+streamTriplet = mapC (f . B.words)+ where+ f [i,j,x] = (readDecimal_ i, readDecimal_ j, decodeElem x)+ f x = error $ "Formatting error: " <> show x+{-# INLINE streamTriplet #-}
src/Data/Matrix/Static/LinearAlgebra.hs view
@@ -18,6 +18,8 @@ , LinearAlgebra(..) -- * Dense matrix operation+ , zeros+ , ones , inverse , eig , svd@@ -28,9 +30,10 @@ import Data.Vector.Storable (Vector) import System.IO.Unsafe (unsafePerformIO) import Data.Complex (Complex)-import Data.Singletons.Prelude hiding ((@@), type (==))+import Data.Singletons.Prelude hiding ((@@), type (==), type (-), type (<=)) import Data.Type.Bool (If)-import Data.Type.Equality (type (==))+import Data.Type.Equality+import GHC.TypeLits import qualified Data.Matrix.Static.Dense as D import qualified Data.Matrix.Static.Sparse as S@@ -100,6 +103,24 @@ class LinearAlgebra (mat :: C.MatrixKind) where ident :: (Numeric a, SingI n) => mat n n Vector a + colSum :: (Numeric a, SingI n, C.Matrix mat Vector a)+ => mat m n Vector a+ -> Matrix 1 n a+ colSum mat = D.create $ do+ m <- CM.replicate 0+ flip C.imapM_ mat $ \(_,j) v -> CM.unsafeModify m (+v) (0, j)+ return m+ {-# INLINE colSum #-}++ rowSum :: (Numeric a, SingI m, C.Matrix mat Vector a)+ => mat m n Vector a+ -> Matrix m 1 a+ rowSum mat = D.create $ do+ m <- CM.replicate 0+ flip C.imapM_ mat $ \(i,_) x -> CM.unsafeModify m (+x) (i, 0)+ return m+ {-# INLINE rowSum #-}+ instance LinearAlgebra D.Matrix where ident = D.diag 0 $ D.replicate 1 @@ -109,14 +130,14 @@ class Factorization mat where -- | Eigenvalues (from largest to smallest) and -- eigenvectors (as columns) of a general square matrix.- eigS :: (SingI k, SingI n, (k <= n - 2) ~ 'True)+ eigS :: (SingI k, SingI n, k <= n - 2) => Sing k -> mat n n Vector Double -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double)) -- | Eigenvalues (from largest to smallest) and -- eigenvectors (as columns) of a symmetric square matrix.- eigSH :: (SingI k, SingI n, (k <= n - 1) ~ 'True)+ eigSH :: (SingI k, SingI n, k <= n - 1) => Sing k -> mat n n Vector Double -> (Matrix k 1 Double, Matrix n k Double)@@ -192,6 +213,14 @@ R Double = Double R (Complex Double) = Double R (Complex Float) = Float++zeros :: (SingI m, SingI n) => Matrix m n Double+zeros = D.replicate 0+{-# INLINE zeros #-}++ones :: (SingI m, SingI n) => Matrix m n Double+ones = D.replicate 1+{-# INLINE ones #-} -- | The inverse of a dense matrix. inverse :: (SingI n, Numeric a) => Matrix n n a -> Matrix n n a
src/Data/Matrix/Static/LinearAlgebra/Internal.hs view
@@ -16,7 +16,6 @@ , unsafeWithS ) where -import Data.Vector.Storable (Storable) import qualified Data.Vector.Storable as VS import qualified Data.Vector.Storable.Mutable as VSM import System.IO.Unsafe (unsafePerformIO)
src/Data/Matrix/Static/Sparse.hs view
@@ -31,12 +31,14 @@ , C.unsafeIndex , C.unsafeTakeRow , C.unsafeTakeColumn+ , unsafeTakeColumnC -- * Construction , C.empty , fromTriplet , fromTripletC , toTriplet+ , withDecodedMatrix , C.fromVector , C.fromList , C.unsafeFromVector@@ -44,6 +46,7 @@ , diagRect -- * Conversions+ , toDense , C.flatten , C.toList @@ -51,6 +54,7 @@ , C.convertAny ) where +import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Storable as S@@ -65,8 +69,13 @@ import GHC.TypeLits (type (<=)) import Foreign.C.Types import Data.Complex+import Data.Store (Store(..), Size(..), decodeExWith)+import Foreign.Storable (sizeOf)+import Data.ByteString (ByteString)+import Data.Singletons.TypeLits import qualified Data.Matrix.Static.Dense as D+import qualified Data.Matrix.Static.Dense.Mutable as DM import qualified Data.Matrix.Static.Generic as C import Data.Matrix.Static.Sparse.Mutable @@ -108,6 +117,38 @@ -- non-zero in the previous two arrays. -> SparseMatrix r c v a +instance (G.Vector v a, Zero a, Store (v a), SingI r, SingI c) =>+ Store (SparseMatrix r c v a) where+ size = VarSize $ \(SparseMatrix nnz inner outer) ->+ case (size, size) of+ (VarSize f, VarSize g) ->+ 2 * sizeOf (0 :: Int) + f nnz + g inner + g outer++ poke mat@(SparseMatrix nnz inner outer) = poke r >> poke c >> poke nnz >>+ poke inner >> poke outer+ where+ (r,c) = C.dim mat+ peek = do+ r' <- peek+ c' <- peek+ if r' /= r || c' /= c+ then error $ "Dimensions donot match: " <> show (r,c) <> " /= " <> show (r',c')+ else SparseMatrix <$> peek <*> peek <*> peek+ where+ r = fromIntegral $ fromSing (sing :: Sing r) :: Int+ c = fromIntegral $ fromSing (sing :: Sing c) :: Int++withDecodedMatrix :: forall v a b. (G.Vector v a, Store (v a))+ => ByteString -> (forall r c. SparseMatrix r c v a -> b) -> b+withDecodedMatrix bs f = withSomeSing (fromIntegral (r :: Int)) $ \(SNat :: Sing r) ->+ withSomeSing (fromIntegral (c :: Int)) $ \(SNat :: Sing c) ->+ f (SparseMatrix nnz inner outer :: SparseMatrix r c v a)+ where+ (r,c,nnz,inner,outer) = decodeExWith + ((,,,,) <$> peek <*> peek <*> peek <*> peek <*> peek)+ bs+{-# INLINE withDecodedMatrix #-}+ instance (G.Vector v a, Eq (v a)) => Eq (SparseMatrix r c v a) where (==) (SparseMatrix a b c) (SparseMatrix a' b' c') = a == a' && b == b' && c == c'@@ -137,12 +178,19 @@ r1 = fromIntegral $ outer `S.unsafeIndex` (j+1) - 1 {-# INLINE unsafeIndex #-} + unsafeTakeColumn mat i = G.create $ do+ vec <- GM.replicate (C.rows mat) zero+ let f (r,_,v) = GM.unsafeWrite vec r v+ runConduit $ unsafeTakeColumnC mat i .| mapM_C f+ return vec+ {-# INLINE unsafeTakeColumn #-}+ -- | O(1) Create matrix from vector containing columns. unsafeFromVector :: forall r c. (G.Vector v a, SingI r, SingI c) => v a -> SparseMatrix r c v a unsafeFromVector vec = fromTriplet vec' where- vec' = map (\((a,b),c) -> (a,b,c)) $ filter ((/=zero) . snd) $+ vec' = V.fromList $ map (\((a,b),c) -> (a,b,c)) $ filter ((/=zero) . snd) $ zipWith (\i x -> (toIndex i, x)) [0..] $ G.toList vec toIndex i = swap $ i `divMod` r r = fromIntegral $ fromSing (sing :: Sing r)@@ -173,35 +221,67 @@ vec' = runST $ runConduit $ toTriplet mat .| mapC g .| sinkVector g (i,j,x) = f (i,j) x {-# INLINE imap #-}- + imapM_ f mat@(SparseMatrix _ _ _) = runConduit $ toTriplet mat .| mapM_C g+ where+ g (i,j,x) = f (i,j) x >> return ()+ {-# INLINE imapM_ #-}++ sequence (SparseMatrix vec inner outer) = do+ vec' <- G.sequence vec+ return $ SparseMatrix vec' inner outer+ {-# INLINE sequence #-}++ sequence_ (SparseMatrix vec _ _) = G.sequence_ vec+ {-# INLINE sequence_ #-}++toDense :: (Zero a, G.Vector v a, SingI r, SingI c)+ => SparseMatrix r c v a -> D.Matrix r c v a+toDense mat = D.create $ do+ m <- DM.replicate zero+ flip C.imapM_ mat $ \idx -> DM.unsafeWrite m idx+ return m+{-# INLINE toDense #-}++-- | Stream a column.+unsafeTakeColumnC :: (Monad m, G.Vector v a)+ => SparseMatrix r c v a -> Int -> ConduitT i (Int, Int, a) m ()+unsafeTakeColumnC (SparseMatrix nnz inner outer) i = enumFromToC lo hi .| mapC f+ where+ f idx = (fromIntegral $ inner `S.unsafeIndex` idx, i, nnz `G.unsafeIndex` idx)+ lo = fromIntegral $ outer S.! i+ hi = fromIntegral $ outer S.! (i+1) - 1+{-# INLINE unsafeTakeColumnC #-}+ -- | O(n) Create matrix from triplet. row and column indices *are not* assumed to be ordered -- duplicate entries are carried over to the CSR represention-fromTriplet :: forall t r c v a. (Traversable t, G.Vector v a, SingI r, SingI c)- => t (Int, Int, a) -> SparseMatrix r c v a+fromTriplet :: forall u r c v a. (G.Vector u (Int, Int, a), G.Vector v a, SingI r, SingI c)+ => u (Int, Int, a) -> SparseMatrix r c v a fromTriplet triplets = SparseMatrix val inner outer where outer = S.scanl (+) 0 $ S.create $ do vec <- SM.replicate c 0- _ <- flip mapM triplets $ \(_, j, _) -> + G.forM_ triplets $ \(_, j, _) -> SM.unsafeModify vec (+1) j return vec (val, inner) = runST $ do outer' <- S.thaw outer val' <- GM.new nnz inner' <- SM.new nnz- _ <- flip mapM triplets $ \(i, j, v) -> do+ G.forM_ triplets $ \(i, j, v) -> do idx <- fromIntegral <$> SM.unsafeRead outer' j GM.unsafeWrite val' idx v SM.unsafeWrite inner' idx $ fromIntegral i SM.unsafeModify outer' (+1) j (,) <$> G.unsafeFreeze val' <*> S.unsafeFreeze inner'- nnz = length triplets+ nnz = G.length triplets c = fromIntegral $ fromSing (sing :: Sing c) {-# INLINE fromTriplet #-} --- | O(n) Create matrix from triplet. row and column indices *are not* assumed to be ordered--- duplicate entries are carried over to the CSR represention+-- | O(n) Create matrix from triplet. Row and column indices *are not* assumed+-- to be ordered. Duplicate entries are carried over to the CSR represention.+-- NOTE: The Conduit will be consumed twice. Use `fromTriplet` if generating+-- the Conduit is expensive. fromTripletC :: forall m r c v a. (Monad m, G.Vector v a, SingI r, SingI c) => ConduitT () (Int, Int, a) m () -> m (SparseMatrix r c v a)@@ -236,6 +316,7 @@ clone x = S.create $ S.thaw x {-# INLINE fromTripletC #-} +-- | Convert sparse matrix to triplets in column order. toTriplet :: (Monad m, G.Vector v a, SingI r, SingI c) => SparseMatrix r c v a -> ConduitT i (Int, Int, a) m () toTriplet (SparseMatrix val inner outer) =@@ -277,7 +358,6 @@ k = (u+l) `shiftR` 1 x' = vec `S.unsafeIndex` k {-# INLINE binarySearchByBounds #-}- ------------------------------------------------------------------------------- -- Helper
tests/Test/Base.hs view
@@ -16,9 +16,11 @@ import Data.Singletons hiding ((@@)) import Data.Singletons.Prelude (Min) import Data.Vector (Vector)+import qualified Data.Vector as V import Control.Monad.IO.Class (liftIO) import Test.Tasty.QuickCheck import Conduit+import Data.Store import Test.Utils @@ -26,13 +28,30 @@ base = testGroup "Base" [ pConversion , pTranspose+ , pSerialization ] +pSerialization :: TestTree+pSerialization = testGroup "Serialization"+ [ testProperty "Dense: id == decode . encode" tStoreD+ , testProperty "Sparse: id == decode . encode" tStoreS+ , testProperty "Sparse: id == decode . encode" tStoreS' ]+ where+ tStoreD :: D.Matrix 80 60 Vector Double -> Bool+ tStoreD mat = mat == decodeEx (encode mat)+ tStoreS :: S.SparseMatrix 80 60 Vector Double -> Bool+ tStoreS mat = mat == decodeEx (encode mat)+ tStoreS' :: S.SparseMatrix 80 60 Vector Double -> Bool+ tStoreS' mat = G.flatten mat == S.withDecodedMatrix (encode mat) G.flatten++pConversion :: TestTree pConversion = testGroup "Conversion"- [ testProperty "Dense: fromVector . flatten" t1- , testProperty "Sparse: fromVector . flatten" t2- , testProperty "Sparse -- fromTriplet . toTriplet" tTri- , testProperty "Sparse -- fromTripletC . toTriplet" tTriC+ [ testProperty "Dense: id == fromVector . flatten" t1+ , testProperty "Sparse: id == fromVector . flatten" t2+ , testProperty "Sparse -- id == fromTriplet . toTriplet" tTri+ , testProperty "Sparse -- id == fromTripletC . toTriplet" tTriC+ , testProperty "Sparse -- id == fromColumn . toColumn" tCol+ , testProperty "Sparse -- toTripletC == mapM_ takeColumnC" tCol' , testProperty "Sparse -- dense" t3 ] where t1 :: D.Matrix 80 60 Vector Int -> Bool@@ -46,10 +65,19 @@ tTri :: S.SparseMatrix 80 60 Vector Int -> Bool tTri mat = S.fromTriplet xs == mat where- xs = runIdentity $ runConduit $ S.toTriplet mat .| sinkList+ xs = V.fromList $ runIdentity $ runConduit $ S.toTriplet mat .| sinkList tTriC :: S.SparseMatrix 80 60 Vector Int -> Bool tTriC mat = mat == runIdentity (S.fromTripletC (S.toTriplet mat))+ tCol :: S.SparseMatrix 80 60 Vector Int -> Bool+ tCol mat = mat == G.fromColumns (map (G.unsafeTakeColumn mat) [0..G.cols mat -1])+ tCol' :: S.SparseMatrix 80 60 Vector Int -> Bool+ tCol' mat = a == b+ where+ a = runIdentity $ runConduit $+ mapM_ (S.unsafeTakeColumnC mat) [0..G.cols mat -1] .| sinkList+ b = runIdentity $ runConduit $ S.toTriplet mat .| sinkList +pTranspose :: TestTree pTranspose = testGroup "Transpose" [ testProperty "Dense" t1 , testProperty "Sparse" tSp
tests/Test/LinearAlgebra.hs view
@@ -18,6 +18,7 @@ import Data.Singletons.Prelude (Min) import Data.Complex import Data.Vector.Storable (Storable)+import qualified Data.Vector.Storable as V import GHC.TypeNats (KnownNat) import Control.Monad.IO.Class (liftIO) import Test.Tasty.QuickCheck@@ -26,10 +27,28 @@ linearAlgebra :: TestTree linearAlgebra = testGroup "Linear algebra"- [ svdTest+ [ basicTest+ , svdTest , eigenTest ] +basicTest :: TestTree+basicTest = testGroup "Basic"+ [ testProperty "Row sum (Dense)" drs+ , testProperty "Column sum (Dense)" dcs+ , testProperty "Row sum (Sparse)" srs+ , testProperty "Column sum (Sparse)" scs ]+ where+ drs :: Matrix 50 30 Double -> Bool+ drs m = D.toList (rowSum m) == map V.sum (G.toRows m)+ dcs :: Matrix 50 30 Double -> Bool+ dcs m = D.toList (colSum m) == map V.sum (G.toColumns m)+ srs :: SparseMatrix 50 30 Double -> Bool+ srs m = D.toList (rowSum m) == map V.sum (G.toRows m)+ scs :: SparseMatrix 50 30 Double -> Bool+ scs m = D.toList (colSum m) == map V.sum (G.toColumns m)++svdTest :: TestTree svdTest = testGroup "SVD" [ testProperty "SVD (Float)" svd1 , testProperty "SVD (Double)" svd2@@ -46,6 +65,7 @@ m' = u @@ S.diag d @@ D.transpose v (u,d,v) = svd m +eigenTest :: TestTree eigenTest = testGroup "Eigendecomposition" [ testProperty "Full" eigen1 , testProperty "Partial dense" eigen2@@ -73,11 +93,3 @@ where m = raw %+% D.transpose raw (d, v)= eigSH (sing :: Sing 99) m----{--propTranspose :: Matrix 50 100 Double -> Bool-propTranspose m = D.transpose (D.transpose m) == m && - D.convertAny (S.transpose $ S.transpose (D.convertAny m)) == m--}
tests/Test/Utils.hs view
@@ -17,6 +17,7 @@ import qualified Data.Matrix.Static.Sparse as S import Data.Vector.Storable (Storable) import Data.List.Ordered+import qualified Data.Vector as V import Data.Ord import qualified Data.Vector.Generic as G import Data.AEq@@ -57,7 +58,7 @@ i <- choose (0, m-1) j <- choose (0, n-1) return ((i,j),v)- return $ S.fromTriplet $ map (\((a,b),c) -> (a,b,c)) xs+ return $ S.fromTriplet $ V.fromList $ map (\((a,b),c) -> (a,b,c)) xs where p = (m * n) `div` 10 m = fromIntegral $ fromSing (sing :: Sing m)