diff --git a/matrix-sized.cabal b/matrix-sized.cabal --- a/matrix-sized.cabal +++ b/matrix-sized.cabal @@ -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 diff --git a/src/Data/Matrix/Static/Dense.hs b/src/Data/Matrix/Static/Dense.hs --- a/src/Data/Matrix/Static/Dense.hs +++ b/src/Data/Matrix/Static/Dense.hs @@ -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 #-} diff --git a/src/Data/Matrix/Static/Dense/Mutable.hs b/src/Data/Matrix/Static/Dense/Mutable.hs --- a/src/Data/Matrix/Static/Dense/Mutable.hs +++ b/src/Data/Matrix/Static/Dense/Mutable.hs @@ -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) diff --git a/src/Data/Matrix/Static/Generic/Mutable.hs b/src/Data/Matrix/Static/Generic/Mutable.hs --- a/src/Data/Matrix/Static/Generic/Mutable.hs +++ b/src/Data/Matrix/Static/Generic/Mutable.hs @@ -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) diff --git a/src/Data/Matrix/Static/IO.hs b/src/Data/Matrix/Static/IO.hs new file mode 100644 --- /dev/null +++ b/src/Data/Matrix/Static/IO.hs @@ -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 #-} diff --git a/src/Data/Matrix/Static/LinearAlgebra.hs b/src/Data/Matrix/Static/LinearAlgebra.hs --- a/src/Data/Matrix/Static/LinearAlgebra.hs +++ b/src/Data/Matrix/Static/LinearAlgebra.hs @@ -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 diff --git a/src/Data/Matrix/Static/LinearAlgebra/Internal.hs b/src/Data/Matrix/Static/LinearAlgebra/Internal.hs --- a/src/Data/Matrix/Static/LinearAlgebra/Internal.hs +++ b/src/Data/Matrix/Static/LinearAlgebra/Internal.hs @@ -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) diff --git a/src/Data/Matrix/Static/Sparse.hs b/src/Data/Matrix/Static/Sparse.hs --- a/src/Data/Matrix/Static/Sparse.hs +++ b/src/Data/Matrix/Static/Sparse.hs @@ -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 diff --git a/tests/Test/Base.hs b/tests/Test/Base.hs --- a/tests/Test/Base.hs +++ b/tests/Test/Base.hs @@ -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 diff --git a/tests/Test/LinearAlgebra.hs b/tests/Test/LinearAlgebra.hs --- a/tests/Test/LinearAlgebra.hs +++ b/tests/Test/LinearAlgebra.hs @@ -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 --} diff --git a/tests/Test/Utils.hs b/tests/Test/Utils.hs --- a/tests/Test/Utils.hs +++ b/tests/Test/Utils.hs @@ -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)