sparse-linear-algebra 0.2.0.0 → 0.2.0.1
raw patch · 3 files changed
+182/−177 lines, 3 filesdep ~basenew-component:exe:sparse-linear-algebraPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: base
API changes (from Hackage documentation)
Files
- app/Main.hs +4/−0
- sparse-linear-algebra.cabal +8/−12
- src/Numeric/LinearAlgebra/Sparse.hs +170/−165
+ app/Main.hs view
@@ -0,0 +1,4 @@+module Main where++main :: IO ()+main = putStrLn $ "2 + 3 = " ++ show (2 + 3)
sparse-linear-algebra.cabal view
@@ -1,5 +1,5 @@ name: sparse-linear-algebra-version: 0.2.0.0+version: 0.2.0.1 synopsis: Numerical computation in native Haskell description: Currently it provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details homepage: https://github.com/ocramz/sparse-linear-algebra@@ -29,17 +29,13 @@ , mwc-random --- executable sparse-linear-algebra--- default-language: Haskell2010--- ghc-options: -threaded -rtsopts -with-rtsopts=-N--- hs-source-dirs: app--- main-is: Main.hs--- build-depends: base--- , mtl >= 2.2.1--- , mwc-random--- , primitive >= 0.6.1.0--- , sparse-linear-algebra--- , transformers >= 0.5.2.0+executable sparse-linear-algebra+ default-language: Haskell2010+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ hs-source-dirs: app+ main-is: Main.hs+ build-depends: base+ test-suite spec default-language: Haskell2010
src/Numeric/LinearAlgebra/Sparse.hs view
@@ -82,31 +82,31 @@ -- ** Norms and related results --- *** squared 2-norm+-- | Squared 2-norm normSq :: (Hilbert f, Num a) => f a -> a normSq v = v `dot` v --- *** L1 norm+-- |L1 norm norm1 :: (Foldable t, Num a, Functor t) => t a -> a norm1 v = sum (fmap abs v) --- *** Euclidean norm+-- |Euclidean norm norm2 :: (Hilbert f, Floating a) => f a -> a norm2 v = sqrt (normSq v) --- *** Lp norm (p > 0)+-- |Lp norm (p > 0) normP :: (Foldable t, Functor t, Floating a) => a -> t a -> a normP p v = sum u**(1/p) where u = fmap (**p) v --- *** infinity-norm+-- |Infinity-norm normInfty :: (Foldable t, Ord a) => t a -> a normInfty = maximum --- *** normalize w.r.t. p-norm (p finite)+-- |Normalize w.r.t. p-norm (p finite) normalize :: (Normed f, Floating a, Eq a) => a -> f a -> f a normalize p v = (1 / norm p v) .* v @@ -116,19 +116,19 @@ --- *** Lp inner product (p > 0)+-- |Lp inner product (p > 0) dotLp :: (Set t, Foldable t, Floating a) => a -> t a -> t a -> a dotLp p v1 v2 = sum u**(1/p) where f a b = (a*b)**p u = liftI2 f v1 v2 --- *** reciprocal+-- |Reciprocal reciprocal :: (Functor f, Fractional b) => f b -> f b reciprocal = fmap recip --- *** scale+-- |Scale scale :: (Num b, Functor f) => b -> f b -> f b scale n = fmap (* n) @@ -138,19 +138,13 @@ --- ** FiniteDim : finite-dimensional objects+-- * FiniteDim : finite-dimensional objects class Additive f => FiniteDim f where type FDSize f :: * dim :: f a -> FDSize f ------ -- | unary dimension-checking bracket withDim :: (FiniteDim f, Show e) => f a@@ -179,16 +173,14 @@ --- ** HasData : accessing inner data (do not export)+-- * HasData : accessing inner data (do not export) class Additive f => HasData f a where type HDData f a :: * dat :: f a -> HDData f a ----- ** Sparse : sparse datastructures+-- * Sparse : sparse datastructures class (FiniteDim f, HasData f a) => Sparse f a where spy :: Fractional b => f a -> b@@ -196,16 +188,13 @@ ------ ** Set : things that behave as sets (e.g. of which we can take the union and the intersection)+-- * Set : things that behave as sets class Functor f => Set f where- -- |union binary lift+ -- |union binary lift : apply function on _union_ of two Sets liftU2 :: (a -> a -> a) -> f a -> f a -> f a - -- |intersection binary lift+ -- |intersection binary lift : apply function on _intersection_ of two Sets liftI2 :: (a -> b -> c) -> f a -> f b -> f c @@ -306,7 +295,6 @@ -- | empty sparse vector (size n, no entries)- zeroSV :: Int -> SpVector a zeroSV n = SV n IM.empty @@ -314,7 +302,7 @@ singletonSV :: a -> SpVector a singletonSV x = SV 1 (IM.singleton 0 x) --- ** Create new sparse vector+-- ** Creation -- | create a sparse vector from an association list while discarding all zero entries mkSpVector :: (Num a, Eq a) => Int -> IM.IntMap a -> SpVector a@@ -352,10 +340,12 @@ | otherwise = error "insertSpVector : index out of bounds" ++-- ** fromList fromListSV :: Int -> [(Int, a)] -> SpVector a fromListSV d iix = SV d (IM.fromList (filter (inBounds0 d . fst) iix )) --- |toList+-- ** toList toListSV :: SpVector a -> [(IM.Key, a)] toListSV sv = IM.toList (dat sv) @@ -376,8 +366,9 @@ show (SV d x) = "SV (" ++ show d ++ ") "++ show (IM.toList x) --- | lookup an index in a SpVector (returns 0 if lookup fails)+-- ** Lookup +-- |lookup an index in a SpVector (returns 0 if lookup fails) lookupDenseSV :: Num a => IM.Key -> SpVector a -> a lookupDenseSV i (SV _ im) = IM.findWithDefault 0 i im @@ -386,7 +377,7 @@ -+-- ** Sub-vectors -- | Tail elements tailSV :: SpVector a -> SpVector a tailSV (SV n sv) = SV (n-1) ta where@@ -420,7 +411,7 @@ --- *** Outer vector product+-- ** Outer vector product outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a outerProdSV v1 v2 = fromListSM (m, n) ixy where@@ -471,18 +462,19 @@ spy = spySM --- | TODO : use semilattice properties instead+-- | Componentwise tuple operations+-- TODO : use semilattice properties instead maxTup, minTup :: Ord t => (t, t) -> (t, t) -> (t, t) maxTup (x1,y1) (x2,y2) = (max x1 x2, max y1 y2) minTup (x1,y1) (x2,y2) = (min x1 x2, min y1 y2) --- | empty matrix of size d+-- | Empty matrix of size d emptySpMatrix :: (Int, Int) -> SpMatrix a emptySpMatrix d = SM d IM.empty --- *** multiply matrix by a scalar+-- *** Multiply matrix by a scalar matScale :: Num a => a -> SpMatrix a -> SpMatrix a matScale a = fmap (*a) @@ -497,8 +489,6 @@ --- ** Matrix metadata- -- type synonyms type Rows = Int type Cols = Int@@ -522,7 +512,10 @@ ff irow row = IM.size row == 1 && IM.size (IM.filterWithKey (\j _ -> j == irow) row) == 1 -+-- |is the matrix orthogonal? i.e. Q^t ## Q == I+isOrthogonalSM :: SpMatrix Double -> Bool+isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where+ rsm = roundZeroOneSM $ transposeSM sm ## sm @@ -533,6 +526,12 @@ immSM :: SpMatrix t -> IM.IntMap (IM.IntMap t) immSM (SM _ imm) = imm +++++-- *** Metadata+ -- | (Number of rows, Number of columns) dimSM :: SpMatrix t -> (Rows, Cols) dimSM (SM d _) = d@@ -549,9 +548,6 @@ ncols :: SpMatrix a -> Cols ncols = snd . dim ----- *** SpMatrix information data SMInfo = SMInfo { smNz :: Int, smSpy :: Double} deriving (Eq, Show) @@ -565,6 +561,7 @@ spySM s = fromIntegral (nzSM s) / fromIntegral (nelSM s) + -- *** Non-zero elements in a row nzRow :: SpMatrix a -> IM.Key -> Int@@ -576,7 +573,7 @@ --- *** Mandwidth bounds (min, max)+-- *** Bandwidth bounds (min, max) bwMinSM :: SpMatrix a -> Int bwMinSM = fst . bwBoundsSM@@ -602,21 +599,15 @@ --- ** Sparse matrix builders --- | Zero SpMatrix of size (m, n)-zeroSM :: Int -> Int -> SpMatrix a-zeroSM m n = SM (m,n) IM.empty --- | Insert an element in a preexisting Spmatrix at the specified indices-insertSpMatrix :: IxRow -> IxCol -> a -> SpMatrix a -> SpMatrix a-insertSpMatrix i j x s- | inBounds02 d (i,j) = SM d $ insertIM2 i j x smd - | otherwise = error "insertSpMatrix : index out of bounds" where- smd = immSM s- d = dim s +-- ** Creation +-- | Zero SpMatrix of size (m, n)+zeroSM :: Int -> Int -> SpMatrix a+zeroSM m n = SM (m,n) IM.empty+ -- | Add to existing SpMatrix using data from list (row, col, value) fromListSM' :: Foldable t => t (IxRow, IxCol, a) -> SpMatrix a -> SpMatrix a fromListSM' iix sm = foldl ins sm iix where@@ -631,20 +622,9 @@ fromListDenseSM :: Int -> [a] -> SpMatrix a fromListDenseSM m ll = fromListSM (m, n) $ denseIxArray2 m ll where n = length ll `div` m- ---- |Convert SpMatrix to list and populate missing entries with 0-toDenseListSM :: Num t => SpMatrix t -> [(IxRow, IxCol, t)]-toDenseListSM m =- [(i, j, m @@ (i, j)) | i <- [0 .. nrows m - 1], j <- [0 .. ncols m- 1]]-------- ** Diagonal matrix+-- *** Diagonal matrix mkDiagonal :: Int -> [a] -> SpMatrix a mkDiagonal n = mkSubDiagonal n 0 @@ -654,10 +634,7 @@ - ----- *** Create Super- or sub- diagonal matrix+-- *** Super- or sub- diagonal matrix mkSubDiagonal :: Int -> Int -> [a] -> SpMatrix a mkSubDiagonal n o xx | abs o < n = if o >= 0@@ -678,6 +655,35 @@ +-- ** Element insertion++-- | Insert an element in a preexisting Spmatrix at the specified indices+insertSpMatrix :: IxRow -> IxCol -> a -> SpMatrix a -> SpMatrix a+insertSpMatrix i j x s+ | inBounds02 d (i,j) = SM d $ insertIM2 i j x smd + | otherwise = error "insertSpMatrix : index out of bounds" where+ smd = immSM s+ d = dim s++++ ++-- |Convert SpMatrix to list and populate missing entries with 0+toDenseListSM :: Num t => SpMatrix t -> [(IxRow, IxCol, t)]+toDenseListSM m =+ [(i, j, m @@ (i, j)) | i <- [0 .. nrows m - 1], j <- [0 .. ncols m- 1]]+++++++++++ -- encode :: (Int, Int) -> (Rows, Cols) -> Int -- encode (nr,_) (i,j) = i + (j * nr) @@ -801,7 +807,7 @@ --- *** Misc. SpMatrix operations+-- ** Misc. SpMatrix operations -- | Left fold over SpMatrix foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b@@ -855,7 +861,7 @@ --- ** Rounding operations (!!!)+-- ** Value rounding -- | Round almost-0 and almost-1 to 0 and 1 respectively roundZeroOneSM :: SpMatrix Double -> SpMatrix Double roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im@@ -923,10 +929,6 @@ ---- -- ** Matrix-matrix product matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a@@ -964,7 +966,7 @@ --- ** Sparsified matrix products+-- *** Sparsified matrix products of two matrices -- | A^T B (#^#) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double@@ -980,12 +982,7 @@ --- * Predicates --- |is the matrix orthogonal? i.e. Q^t ## Q == I-isOrthogonalSM :: SpMatrix Double -> Bool-isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where- rsm = roundZeroOneSM $ transposeSM sm ## sm @@ -994,8 +991,7 @@ ---- ** Condition number+-- * Matrix condition number -- |uses the R matrix from the QR factorization conditionNumberSM :: SpMatrix Double -> Double@@ -1013,7 +1009,7 @@ --- ** Householder transformation+-- * Householder transformation hhMat :: Num a => a -> SpVector a -> SpMatrix a hhMat beta x = eye n ^-^ scale beta (x >< x) where@@ -1036,7 +1032,7 @@ --- ** Givens rotation matrix+-- * Givens rotation matrix hypot :: Floating a => a -> a -> a@@ -1097,7 +1093,7 @@ --- ** QR decomposition+-- * QR decomposition -- | Applies Givens rotation iteratively to zero out sub-diagonal elements@@ -1131,9 +1127,9 @@ --- ** Eigenvalue algorithms+-- * Eigenvalue algorithms --- *** All eigenvalues using QR algorithm+-- ** All eigenvalues (QR algorithm) eigsQR :: Int -> SpMatrix Double -> SpVector Double@@ -1149,7 +1145,7 @@ --- *** One eigenvalue and corresponding eigenvector, using Rayleigh iteration +-- ** One eigenvalue and eigenvector (Rayleigh iteration) -- | Cubic-order convergence, but it requires a mildly educated guess on the initial eigenpair rayleighStep ::@@ -1174,7 +1170,7 @@ --- ** Householder vector (G & VL Alg. 5.1.1, function `house`)+-- * Householder vector (G & VL Alg. 5.1.1, function `house`) hhV :: SpVector Double -> (SpVector Double, Double) hhV x = (v, beta) where@@ -1227,7 +1223,7 @@ --- * LINEAR SOLVERS : solve A x = b+-- * Iterative linear solvers -- | numerical tolerance for e.g. solution convergence eps :: Double@@ -1339,7 +1335,7 @@ --- * LINEAR SOLVERS INTERFACE+-- * Linear solver interface data LinSolveMethod = CGS_ | BICGSTAB_ deriving (Eq, Show) @@ -1402,86 +1398,17 @@ --- ** Pretty printing of SpVector and SpMatrix -sizeStr :: SpMatrix a -> String-sizeStr sm =- unwords ["(",show (nrows sm),"rows,",show (ncols sm),"columns ) ,",show nz,"NZ ( sparsity",show sy,")"] where- (SMInfo nz sy) = infoSM sm -showNonZero :: (Show a, Num a, Eq a) => a -> String-showNonZero x = if x == 0 then " " else show x -toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]-toDenseRow (SM (_,ncol) im) irow =- fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1] -toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String-toDenseRowClip sm irow ncomax- | ncols sm > ncomax = unwords (map show h) ++ " ... " ++ show t- | otherwise = show dr- where dr = toDenseRow sm irow- h = take (ncomax - 2) dr- t = last dr -newline :: IO ()-newline = putStrLn "" -printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()-printDenseSM sm = do- newline- putStrLn $ sizeStr sm- newline- printDenseSM' sm 5 5- newline- where - printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()- printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where- rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]- rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]- | otherwise = rr_ --toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String-toDenseListClip sv ncomax- | dim sv > ncomax = unwords (map show h) ++ " ... " ++ show t- | otherwise = show dr- where dr = toDenseListSV sv- h = take (ncomax - 2) dr- t = last dr--printDenseSV :: (Show t, Num t) => SpVector t -> IO ()-printDenseSV sv = do- newline- printDenseSV' sv 5- newline where- printDenseSV' v nco = putStrLn rr_' where- rr_ = toDenseListClip v nco :: String- rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]- | otherwise = rr_---- ** Pretty printer typeclass-class PrintDense a where- prd :: a -> IO ()--instance (Show a, Num a) => PrintDense (SpVector a) where- prd = printDenseSV--instance (Show a, Num a) => PrintDense (SpMatrix a) where- prd = printDenseSM------------ ** Control primitives for bounded iteration with convergence check+-- * Control primitives for bounded iteration with convergence check -- | transform state until a condition is met modifyUntil :: MonadState s m => (s -> Bool) -> (s -> s) -> m s@@ -1579,7 +1506,7 @@ --- ** Rounding operations+-- * Rounding operations -- | Rounding rule almostZero, almostOne :: Double -> Bool@@ -1609,7 +1536,7 @@ --- *** Random matrices and vectors+-- * Random matrices and vectors -- |Dense SpMatrix randMat :: PrimMonad m => Int -> m (SpMatrix Double)@@ -1653,8 +1580,86 @@ +-- * Pretty printing --- *** Misc. utilities++sizeStr :: SpMatrix a -> String+sizeStr sm =+ unwords ["(",show (nrows sm),"rows,",show (ncols sm),"columns ) ,",show nz,"NZ ( sparsity",show sy,")"] where+ (SMInfo nz sy) = infoSM sm +++showNonZero :: (Show a, Num a, Eq a) => a -> String+showNonZero x = if x == 0 then " " else show x+++toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]+toDenseRow (SM (_,ncol) im) irow =+ fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1]++toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String+toDenseRowClip sm irow ncomax+ | ncols sm > ncomax = unwords (map show h) ++ " ... " ++ show t+ | otherwise = show dr+ where dr = toDenseRow sm irow+ h = take (ncomax - 2) dr+ t = last dr++newline :: IO ()+newline = putStrLn ""++printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()+printDenseSM sm = do+ newline+ putStrLn $ sizeStr sm+ newline+ printDenseSM' sm 5 5+ newline+ where + printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()+ printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where+ rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]+ rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]+ | otherwise = rr_+++toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String+toDenseListClip sv ncomax+ | dim sv > ncomax = unwords (map show h) ++ " ... " ++ show t+ | otherwise = show dr+ where dr = toDenseListSV sv+ h = take (ncomax - 2) dr+ t = last dr++printDenseSV :: (Show t, Num t) => SpVector t -> IO ()+printDenseSV sv = do+ newline+ printDenseSV' sv 5+ newline where+ printDenseSV' v nco = putStrLn rr_' where+ rr_ = toDenseListClip v nco :: String+ rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]+ | otherwise = rr_++-- ** Pretty printer typeclass+class PrintDense a where+ prd :: a -> IO ()++instance (Show a, Num a) => PrintDense (SpVector a) where+ prd = printDenseSV++instance (Show a, Num a) => PrintDense (SpMatrix a) where+ prd = printDenseSM++++ ++++++-- * Misc. utilities -- | integer-indexed ziplist denseIxArray :: [b] -> [(Int, b)]