packages feed

dense-linear-algebra (empty) → 0.1.0.0

raw patch · 10 files changed

+567/−0 lines, 10 filesdep +QuickCheckdep +basedep +deepseqsetup-changed

Dependencies added: QuickCheck, base, deepseq, dense-linear-algebra, hspec, math-functions, primitive, vector, vector-algorithms, vector-binary-instances, vector-th-unbox

Files

+ LICENSE view
@@ -0,0 +1,26 @@+Copyright (c) 2009, 2010 Bryan O'Sullivan+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,5 @@+# statistics-dense-linear-algebra++[![Build Status](https://travis-ci.org/githubuser/statistics-dense-linear-algebra.png)](https://travis-ci.org/githubuser/statistics-dense-linear-algebra)++The dense linear algebra functionality and related modules, extracted from `statistics-0.14.0.2`
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ dense-linear-algebra.cabal view
@@ -0,0 +1,50 @@+name:                dense-linear-algebra+version:             0.1.0.0+synopsis:            Simple and incomplete pure haskell implementation of linear algebra+description:+  This library is simply collection of linear-algebra related modules+  split from statistics library.++license:             BSD2+license-file:        LICENSE+author:              Bryan O'Sullivan <bos@serpentine.com>+maintainer:          Alexey Khudaykov <alexey.skladnoy@gmail.com>+copyright:           2018 Author name here+category:            Math, Statistics, Numeric+build-type:          Simple+extra-source-files:  README.md+cabal-version:       >=1.10+tested-with: GHC==7.6.3, GHC==7.8.3, GHC==7.10.3, GHC==8.0.1++library+  default-language:    Haskell2010+  ghc-options:         -Wall+  hs-source-dirs:      src+  exposed-modules:     Statistics.Matrix+                       Statistics.Matrix.Algorithms+                       Statistics.Matrix.Function+                       Statistics.Matrix.Mutable+                       Statistics.Matrix.Types+  build-depends:       base                    >= 4.5 && < 5+                     , deepseq                 >= 1.1.0.2+                     , math-functions          >= 0.1.7                     +                     , primitive               >= 0.3+                     , vector                  >= 0.10+                     , vector-algorithms       >= 0.4+                     , vector-th-unbox+                     , vector-binary-instances >= 0.2.1++test-suite spec+  default-language:    Haskell2010+  ghc-options:         -Wall+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             LibSpec.hs+  build-depends:       base+                     , dense-linear-algebra+                     , hspec+                     , QuickCheck++source-repository head+  type:     git+  location: https://github.com/bos/statistics
+ src/Statistics/Matrix.hs view
@@ -0,0 +1,265 @@+{-# LANGUAGE PatternGuards #-}+-- |+-- Module    : Statistics.Matrix+-- Copyright : 2011 Aleksey Khudyakov, 2014 Bryan O'Sullivan+-- License   : BSD3+--+-- Basic matrix operations.+--+-- There isn't a widely used matrix package for Haskell yet, so+-- we implement the necessary minimum here.++module Statistics.Matrix+    ( -- * Data types+      Matrix(..)+    , Vector+      -- * Conversion from/to lists/vectors+    , fromVector+    , fromList+    , fromRowLists+    , fromRows+    , fromColumns+    , toVector+    , toList+    , toRows+    , toColumns+    , toRowLists+      -- * Other+    , generate+    , generateSym+    , ident+    , diag+    , dimension+    , center+    , multiply+    , multiplyV+    , transpose+    , power+    , norm+    , column+    , row+    , map+    , for+    , unsafeIndex+    , hasNaN+    , bounds+    , unsafeBounds+    ) where++import Prelude hiding (exponent, map)+import Control.Applicative ((<$>))+import Control.Monad.ST+import qualified Data.Vector.Unboxed as U+import           Data.Vector.Unboxed   ((!))+import qualified Data.Vector.Unboxed.Mutable as UM++import Numeric.Sum         (sumVector,kbn)+import Statistics.Matrix.Function+import Statistics.Matrix.Types+import Statistics.Matrix.Mutable  (unsafeNew,unsafeWrite,unsafeFreeze)++++----------------------------------------------------------------+-- Conversion to/from vectors/lists+----------------------------------------------------------------++-- | Convert from a row-major list.+fromList :: Int                 -- ^ Number of rows.+         -> Int                 -- ^ Number of columns.+         -> [Double]            -- ^ Flat list of values, in row-major order.+         -> Matrix+fromList r c = fromVector r c . U.fromList++-- | create a matrix from a list of lists, as rows+fromRowLists :: [[Double]] -> Matrix+fromRowLists = fromRows . fmap U.fromList++-- | Convert from a row-major vector.+fromVector :: Int               -- ^ Number of rows.+           -> Int               -- ^ Number of columns.+           -> U.Vector Double   -- ^ Flat list of values, in row-major order.+           -> Matrix+fromVector r c v+  | r*c /= len = error "input size mismatch"+  | otherwise  = Matrix r c v+  where len    = U.length v++-- | create a matrix from a list of vectors, as rows+fromRows :: [Vector] -> Matrix+fromRows xs+  | [] <- xs        = error "Statistics.Matrix.fromRows: empty list of rows!"+  | any (/=nCol) ns = error "Statistics.Matrix.fromRows: row sizes do not match"+  | nCol == 0       = error "Statistics.Matrix.fromRows: zero columns in matrix"+  | otherwise       = fromVector nRow nCol (U.concat xs)+  where+    nCol:ns = U.length <$> xs+    nRow    = length xs+++-- | create a matrix from a list of vectors, as columns+fromColumns :: [Vector] -> Matrix+fromColumns = transpose . fromRows++-- | Convert to a row-major flat vector.+toVector :: Matrix -> U.Vector Double+toVector (Matrix _ _ v) = v++-- | Convert to a row-major flat list.+toList :: Matrix -> [Double]+toList = U.toList . toVector++-- | Convert to a list of lists, as rows+toRowLists :: Matrix -> [[Double]]+toRowLists (Matrix _ nCol v)+  = chunks $ U.toList v+  where+    chunks [] = []+    chunks xs = case splitAt nCol xs of+      (rowE,rest) -> rowE : chunks rest+++-- | Convert to a list of vectors, as rows+toRows :: Matrix -> [Vector]+toRows (Matrix _ nCol v) = chunks v+  where+    chunks xs+      | U.null xs = []+      | otherwise = case U.splitAt nCol xs of+          (rowE,rest) -> rowE : chunks rest++-- | Convert to a list of vectors, as columns+toColumns :: Matrix -> [Vector]+toColumns = toRows . transpose++++----------------------------------------------------------------+-- Other+----------------------------------------------------------------++-- | Generate matrix using function+generate :: Int                 -- ^ Number of rows+         -> Int                 -- ^ Number of columns+         -> (Int -> Int -> Double)+            -- ^ Function which takes /row/ and /column/ as argument.+         -> Matrix+generate nRow nCol f+  = Matrix nRow nCol $ U.generate (nRow*nCol) $ \i ->+      let (r,c) = i `quotRem` nCol in f r c++-- | Generate symmetric square matrix using function+generateSym+  :: Int                 -- ^ Number of rows and columns+  -> (Int -> Int -> Double)+     -- ^ Function which takes /row/ and /column/ as argument. It must+     --   be symmetric in arguments: @f i j == f j i@+  -> Matrix+generateSym n f = runST $ do+  m <- unsafeNew n n+  for 0 n $ \r -> do+    unsafeWrite m r r (f r r)+    for (r+1) n $ \c -> do+      let x = f r c+      unsafeWrite m r c x+      unsafeWrite m c r x+  unsafeFreeze m+++-- | Create the square identity matrix with given dimensions.+ident :: Int -> Matrix+ident n = diag $ U.replicate n 1.0++-- | Create a square matrix with given diagonal, other entries default to 0+diag :: Vector -> Matrix+diag v+  = Matrix n n $ U.create $ do+      arr <- UM.replicate (n*n) 0+      for 0 n $ \i ->+        UM.unsafeWrite arr (i*n + i) (v ! i)+      return arr+  where+    n = U.length v++-- | Return the dimensions of this matrix, as a (row,column) pair.+dimension :: Matrix -> (Int, Int)+dimension (Matrix r c _) = (r, c)++-- | Matrix-matrix multiplication. Matrices must be of compatible+-- sizes (/note: not checked/).+multiply :: Matrix -> Matrix -> Matrix+multiply m1@(Matrix r1 _ _) m2@(Matrix _ c2 _) =+  Matrix r1 c2 $ U.generate (r1*c2) go+  where+    go t = sumVector kbn $ U.zipWith (*) (row m1 i) (column m2 j)+      where (i,j) = t `quotRem` c2++-- | Matrix-vector multiplication.+multiplyV :: Matrix -> Vector -> Vector+multiplyV m v+  | cols m == c = U.generate (rows m) (sumVector kbn . U.zipWith (*) v . row m)+  | otherwise   = error $ "matrix/vector unconformable " ++ show (cols m,c)+  where c = U.length v++-- | Raise matrix to /n/th power. Power must be positive+-- (/note: not checked).+power :: Matrix -> Int -> Matrix+power mat 1 = mat+power mat n = res+  where+    mat2 = power mat (n `quot` 2)+    pow  = multiply mat2 mat2+    res | odd n     = multiply pow mat+        | otherwise = pow++-- | Element in the center of matrix (not corrected for exponent).+center :: Matrix -> Double+center mat@(Matrix r c _) =+    unsafeBounds U.unsafeIndex mat (r `quot` 2) (c `quot` 2)++-- | Calculate the Euclidean norm of a vector.+norm :: Vector -> Double+norm = sqrt . sumVector kbn . U.map square++-- | Return the given column.+column :: Matrix -> Int -> Vector+column (Matrix r c v) i = U.backpermute v $ U.enumFromStepN i c r+{-# INLINE column #-}++-- | Return the given row.+row :: Matrix -> Int -> Vector+row (Matrix _ c v) i = U.slice (c*i) c v++unsafeIndex :: Matrix+            -> Int              -- ^ Row.+            -> Int              -- ^ Column.+            -> Double+unsafeIndex = unsafeBounds U.unsafeIndex++-- | Apply function to every element of matrix+map :: (Double -> Double) -> Matrix -> Matrix+map f (Matrix r c v) = Matrix r c (U.map f v)++-- | Indicate whether any element of the matrix is @NaN@.+hasNaN :: Matrix -> Bool+hasNaN = U.any isNaN . toVector++-- | Given row and column numbers, calculate the offset into the flat+-- row-major vector.+bounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r+bounds k (Matrix rs cs v) r c+  | r < 0 || r >= rs = error "row out of bounds"+  | c < 0 || c >= cs = error "column out of bounds"+  | otherwise        = k v $! r * cs + c+{-# INLINE bounds #-}++-- | Given row and column numbers, calculate the offset into the flat+-- row-major vector, without checking.+unsafeBounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r+unsafeBounds k (Matrix _ cs v) r c = k v $! r * cs + c+{-# INLINE unsafeBounds #-}++transpose :: Matrix -> Matrix+transpose m@(Matrix r0 c0 _) = Matrix c0 r0 . U.generate (r0*c0) $ \i ->+  let (r,c) = i `quotRem` r0+  in unsafeIndex m c r
+ src/Statistics/Matrix/Algorithms.hs view
@@ -0,0 +1,42 @@+-- |+-- Module    : Statistics.Matrix.Algorithms+-- Copyright : 2014 Bryan O'Sullivan+-- License   : BSD3+--+-- Useful matrix functions.++module Statistics.Matrix.Algorithms+    (+      qr+    ) where++import Control.Applicative ((<$>), (<*>))+import Control.Monad.ST (ST, runST)+import Prelude hiding (replicate)+import Numeric.Sum       (sumVector,kbn)+import Statistics.Matrix (Matrix, column, dimension, for, norm)+import qualified Statistics.Matrix.Mutable as M+import qualified Data.Vector.Unboxed as U++-- | /O(r*c)/ Compute the QR decomposition of a matrix.+-- The result returned is the matrices (/q/,/r/).+qr :: Matrix -> (Matrix, Matrix)+qr mat = runST $ do+  let (m,n) = dimension mat+  r <- M.replicate n n 0+  a <- M.thaw mat+  for 0 n $ \j -> do+    cn <- M.immutably a $ \aa -> norm (column aa j)+    M.unsafeWrite r j j cn+    for 0 m $ \i -> M.unsafeModify a i j (/ cn)+    for (j+1) n $ \jj -> do+      p <- innerProduct a j jj+      M.unsafeWrite r j jj p+      for 0 m $ \i -> do+        aij <- M.unsafeRead a i j+        M.unsafeModify a i jj $ subtract (p * aij)+  (,) <$> M.unsafeFreeze a <*> M.unsafeFreeze r++innerProduct :: M.MMatrix s -> Int -> Int -> ST s Double+innerProduct mmat j k = M.immutably mmat $ \mat ->+  sumVector kbn $ U.zipWith (*) (column mat j) (column mat k)
+ src/Statistics/Matrix/Function.hs view
@@ -0,0 +1,15 @@+{-# LANGUAGE BangPatterns #-}+-- |+module Statistics.Matrix.Function where++-- | Multiply a number by itself.+square :: Double -> Double+square x = x * x++-- | Simple for loop.  Counts from /start/ to /end/-1.+for :: Monad m => Int -> Int -> (Int -> m ()) -> m ()+for n0 !n f = loop n0+  where+    loop i | i == n    = return ()+           | otherwise = f i >> loop (i+1)+{-# INLINE for #-}
+ src/Statistics/Matrix/Mutable.hs view
@@ -0,0 +1,86 @@+-- |+-- Module    : Statistics.Matrix.Mutable+-- Copyright : (c) 2014 Bryan O'Sullivan+-- License   : BSD3+--+-- Basic mutable matrix operations.++module Statistics.Matrix.Mutable+    (+      MMatrix(..)+    , MVector+    , replicate+    , thaw+    , bounds+    , unsafeNew+    , unsafeFreeze+    , unsafeRead+    , unsafeWrite+    , unsafeModify+    , immutably+    , unsafeBounds+    ) where++import Control.Applicative ((<$>))+import Control.DeepSeq (NFData(..))+import Control.Monad.ST (ST)+import Statistics.Matrix.Types (Matrix(..), MMatrix(..), MVector)+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as M+import Prelude hiding (replicate)++replicate :: Int -> Int -> Double -> ST s (MMatrix s)+replicate r c k = MMatrix r c <$> M.replicate (r*c) k++thaw :: Matrix -> ST s (MMatrix s)+thaw (Matrix r c v) = MMatrix r c <$> U.thaw v++unsafeFreeze :: MMatrix s -> ST s Matrix+unsafeFreeze (MMatrix r c mv) = Matrix r c <$> U.unsafeFreeze mv++-- | Allocate new matrix. Matrix content is not initialized hence unsafe.+unsafeNew :: Int                -- ^ Number of row+          -> Int                -- ^ Number of columns+          -> ST s (MMatrix s)+unsafeNew r c+  | r < 0     = error "Statistics.Matrix.Mutable.unsafeNew: negative number of rows"+  | c < 0     = error "Statistics.Matrix.Mutable.unsafeNew: negative number of columns"+  | otherwise = do+      vec <- M.new (r*c)+      return $ MMatrix r c vec++unsafeRead :: MMatrix s -> Int -> Int -> ST s Double+unsafeRead mat r c = unsafeBounds mat r c M.unsafeRead+{-# INLINE unsafeRead #-}++unsafeWrite :: MMatrix s -> Int -> Int -> Double -> ST s ()+unsafeWrite mat row col k = unsafeBounds mat row col $ \v i ->+  M.unsafeWrite v i k+{-# INLINE unsafeWrite #-}++unsafeModify :: MMatrix s -> Int -> Int -> (Double -> Double) -> ST s ()+unsafeModify mat row col f = unsafeBounds mat row col $ \v i -> do+  k <- M.unsafeRead v i+  M.unsafeWrite v i (f k)+{-# INLINE unsafeModify #-}++-- | Given row and column numbers, calculate the offset into the flat+-- row-major vector.+bounds :: MMatrix s -> Int -> Int -> (MVector s -> Int -> r) -> r+bounds (MMatrix rs cs mv) r c k+  | r < 0 || r >= rs = error "row out of bounds"+  | c < 0 || c >= cs = error "column out of bounds"+  | otherwise        = k mv $! r * cs + c+{-# INLINE bounds #-}++-- | Given row and column numbers, calculate the offset into the flat+-- row-major vector, without checking.+unsafeBounds :: MMatrix s -> Int -> Int -> (MVector s -> Int -> r) -> r+unsafeBounds (MMatrix _ cs mv) r c k = k mv $! r * cs + c+{-# INLINE unsafeBounds #-}++immutably :: NFData a => MMatrix s -> (Matrix -> a) -> ST s a+immutably mmat f = do+  k <- f <$> unsafeFreeze mmat+  rnf k `seq` return k+{-# INLINE immutably #-}
+ src/Statistics/Matrix/Types.hs view
@@ -0,0 +1,60 @@+-- |+-- Module    : Statistics.Matrix.Types+-- Copyright : 2014 Bryan O'Sullivan+-- License   : BSD3+--+-- Basic matrix operations.+--+-- There isn't a widely used matrix package for Haskell yet, so+-- we implement the necessary minimum here.++module Statistics.Matrix.Types+    (+      Vector+    , MVector+    , Matrix(..)+    , MMatrix(..)+    , debug+    ) where++import Data.Char (isSpace)+import Numeric (showFFloat)+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as M++type Vector = U.Vector Double+type MVector s = M.MVector s Double++-- | Two-dimensional matrix, stored in row-major order.+data Matrix = Matrix {+      rows     :: {-# UNPACK #-} !Int -- ^ Rows of matrix.+    , cols     :: {-# UNPACK #-} !Int -- ^ Columns of matrix.+    , _vector  :: !Vector             -- ^ Matrix data.+    } deriving (Eq)++-- | Two-dimensional mutable matrix, stored in row-major order.+data MMatrix s = MMatrix+                 {-# UNPACK #-} !Int+                 {-# UNPACK #-} !Int+                 !(MVector s)++-- The Show instance is useful only for debugging.+instance Show Matrix where+    show = debug++debug :: Matrix -> String+debug (Matrix r c vs) = unlines $ zipWith (++) (hdr0 : repeat hdr) rrows+  where+    rrows         = map (cleanEnd . unwords) . split $ zipWith (++) ldone tdone+    hdr0          = show (r,c) ++ " "+    hdr           = replicate (length hdr0) ' '+    pad plus k xs = replicate (k - length xs) ' ' `plus` xs+    ldone         = map (pad (++) (longest lstr)) lstr+    tdone         = map (pad (flip (++)) (longest tstr)) tstr+    (lstr, tstr)  = unzip . map (break (=='.') . render) . U.toList $ vs+    longest       = maximum . map length+    render k      = reverse . dropWhile (=='.') . dropWhile (=='0') . reverse .+                    showFFloat (Just 4) k $ ""+    split []      = []+    split xs      = i : split rest where (i, rest) = splitAt c xs+    cleanEnd      = reverse . dropWhile isSpace . reverse
+ test/LibSpec.hs view
@@ -0,0 +1,16 @@+module Main where++import Test.Hspec+import Test.Hspec.QuickCheck+++main :: IO ()+main = hspec spec++spec :: Spec+spec =+  describe "Lib" $ do+    it "works" $ do+      True `shouldBe` True+    -- prop "ourAdd is commutative" $ \x y ->+    --   ourAdd x y `shouldBe` ourAdd y x