matrix 0.3.6.1 → 0.3.6.2

raw patch · 5 files changed

+181/−114 lines, 5 filesdep ~basedep ~deepseqPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base, deepseq

API changes (from Hackage documentation)

Files

• Data/Matrix.hs +56/−48
• bench/rref.hs +19/−0
• changelog.md +2/−0
• matrix.cabal +74/−65
• test/Examples.hs +30/−1

@@ -83,7 +83,7 @@ import GHC.Generics (Generic) -- Data import           Control.Monad.Primitive (PrimMonad, PrimState)-import           Data.List               (maximumBy,foldl1')+import           Data.List               (maximumBy,foldl1',find) import           Data.Ord                (comparing) import qualified Data.Vector             as V import qualified Data.Vector.Mutable     as MV@@ -339,7 +339,16 @@          -> [a] -- ^ List of elements          -> Matrix a {-# INLINE fromList #-}-fromList n m = M n m 0 0 m . V.fromListN (n*m)+fromList n m xs+    | n*m > V.length v =+        (error $+            "List size "+            ++ show (V.length v)+            ++ " is inconsistent with matrix size "+            ++ sizeStr n m+            ++ " in fromList")+    | otherwise       = M n m 0 0 m v+    where v = V.fromListN (n*m) xs  -- | Get the elements of a matrix stored in a list. --@@ -584,61 +593,60 @@         let             adjoinedWId = m <|> identity (nrows m)             rref'd = rref adjoinedWId-        in rref'd >>= return . submatrix 1 (nrows m) (ncols m + 1) (ncols m * 2)+            checkInvertible a = if unsafeGet (ncols m) (nrows m) a == 1+                then Right a+                else Left "Attempt to invert a non-invertible matrix"+        in rref'd >>= checkInvertible >>= return . submatrix 1 (nrows m) (ncols m + 1) (ncols m * 2) --- | /O(rows*rows*cols*cols)/. Converts a matrix to reduced row echelon form, thus---  solving a linear system of equations. This requires that (cols > rows)---  if cols < rows, then there are fewer variables than equations and the---  problem cannot be solved consistently. If rows = cols, then it is---  basically a homogenous system of equations, so it will be reduced to---  identity or an error depending on whether the marix is invertible---  (this case is allowed for robustness).++-- | Converts a matrix to reduced row echelon form, thus+--   solving a linear system of equations. This requires that (cols > rows)+--   if cols < rows, then there are fewer variables than equations and the+--   problem cannot be solved consistently. If rows = cols, then it is+--   basically a homogenous system of equations, so it will be reduced to+--   identity or an error depending on whether the marix is invertible+--   (this case is allowed for robustness).+--   This implementation is taken from rosettacode+--   https://rosettacode.org/wiki/Reduced_row_echelon_form#Haskell rref :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a) rref m         | ncols m < nrows m             = Left $                 "Invalid dimensions "                     ++ show (sizeStr (ncols m) (nrows m))-                    ++ "; the number of columns must be greater than or equal to the number of rows"-        | otherwise             = rrefRefd =<< (ref m)-    where-    rrefRefd mtx-        | nrows mtx == 1    = Right mtx-        | otherwise =-            let-                -- this is super-slow: [resolvedRight] is cubic because [combineRows] is quadratic-                resolvedRight = foldr (.) id (map resolveRow [1..col-1]) mtx-                    where-                    col = nrows mtx-                    resolveRow n = combineRows n (-getElem n col mtx) col-                top = submatrix 1 (nrows resolvedRight - 1) 1 (ncols resolvedRight) resolvedRight-                top' = rrefRefd top-                bot = submatrix (nrows resolvedRight) (nrows resolvedRight) 1 (ncols resolvedRight) resolvedRight-            in top' >>= return . (<-> bot)+        | otherwise = Right . fromLists $ f matM 0 [0 .. rows - 1]+  where+    matM = toLists m+    rows = nrows m+    cols = ncols m +    f a _    []           = a+    f a lead (r : rs)+      | isNothing indices = a+      | otherwise         = f a' (lead' + 1) rs+      where+        indices = find p l+        p (col, row) = a !! row !! col /= 0+        l = [(col, row) |+            col <- [lead .. cols - 1],+            row <- [r .. rows - 1]] -ref :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a)-ref mtx-        | nrows mtx == 1-            = clearedLeft-        | otherwise = do -                (tl, tr, bl, br) <- (splitBlocks 1 1 <$> clearedLeft)-                br' <- ref br-                return  ((tl <|> tr) <-> (bl <|> br'))-    where-    sigAtTop = (\row -> switchRows 1 row mtx) <$> goodRow-        where-        significantRow n = getElem n 1 mtx /= 0-        goodRow = case listToMaybe (filter significantRow [1..nrows mtx]) of-            Nothing -> Left "Attempt to invert a non-invertible matrix"-            Just x -> return x-    normalizedFirstRow = (\sigAtTop' -> scaleRow (1 / getElem 1 1 sigAtTop') 1 sigAtTop') <$> sigAtTop-    clearedLeft =  do -            comb <- mapM combinator [2..nrows mtx]-            firstRow <- normalizedFirstRow-            return $ (foldr (.) id comb) firstRow-        where-        combinator n = (\normalizedFirstRow'  ->combineRows n (-getElem n 1 normalizedFirstRow') 1) <$> normalizedFirstRow+        Just (lead', i) = indices+        newRow = map (/ a !! i !! lead') $ a !! i++        a' = zipWith g [0..] $+            replace r newRow $+            replace i (a !! r) a+        g n row+            | n == r    = row+            | otherwise = zipWith h newRow row+              where h = subtract . (* row !! lead')++        replace :: Int -> b -> [b] -> [b]+        {- Replaces the element at the given index. -}+        replace n e t = a ++ e : b+          where (a, _ : b) = splitAt n t+  -- | Extend a matrix to a given size adding a default element. --   If the matrix already has the required size, nothing happens.

@@ -0,0 +1,19 @@++import Criterion.Main+--+import Data.Matrix++mat :: Int -> Matrix Double+mat n = fromList n n $ map (^2) [1..]++testrref :: Int -> Either String (Matrix Double)+testrref n = rref (mat n)+++brref :: Int -> Benchmark+brref n = bgroup ("rref" ++ show n)+ [ bench "rref" $ nf testrref n+ ]++main :: IO ()+main = defaultMain $ fmap brref [10,25,100]

@@ -0,0 +1,2 @@+## 0.3.6.2+* Fixes for `rref` and `fromList`.

@@ -1,65 +1,74 @@-Name: matrix
-Version: 0.3.6.1
-Author: Daniel Díaz
-Category: Math
-Build-type: Simple
-License: BSD3
-License-file: license
-Maintainer: Daniel Díaz (dhelta `dot` diaz `at` gmail `dot` com)
-Stability: In development
-Bug-reports: https://github.com/Daniel-Diaz/matrix/issues
-Synopsis: A native implementation of matrix operations.
-Description:
-  Matrix library. Basic operations and some algorithms.
-  .
-  To get the library update your cabal package list (if needed) with @cabal update@ and
-  then use @cabal install matrix@, assuming that you already have Cabal installed.
-  Usage examples are included in the API reference generated by Haddock.
-  .
-  If you want to use GSL, BLAS and LAPACK, @hmatrix@ (<http://hackage.haskell.org/package/hmatrix>)
-  is the way to go.
-Cabal-version: >= 1.8
-Extra-source-files:
-  readme.md
-  -- Benchmarks
-  bench/mult.hs
-
-Source-repository head
-  type: git
-  location: git://github.com/Daniel-Diaz/matrix.git
-
-Library
-  Build-depends: base >= 4.5 && < 5
-               , vector >= 0.10
-               , deepseq >= 1.3.0.0 && < 1.5
-               , primitive >= 0.5
-               , semigroups >= 0.9
-               , loop >= 0.2
-  Exposed-modules: Data.Matrix
-  GHC-Options: -Wall -O2 -fstatic-argument-transformation
-
-Benchmark matrix-mult
-  type: exitcode-stdio-1.0
-  hs-source-dirs: bench
-  main-is: mult.hs
-  build-depends: base == 4.*
-               , matrix
-               , criterion
-  ghc-options: -O2 -fstatic-argument-transformation
-
-Test-Suite matrix-test
-  type: exitcode-stdio-1.0
-  hs-source-dirs: test
-  main-is: Main.hs
-  build-depends: base == 4.*
-               , matrix
-               , tasty
-               , QuickCheck
-               , tasty-quickcheck
-
-Test-Suite matrix-examples
-  type: exitcode-stdio-1.0
-  hs-source-dirs: test
-  main-is: Examples.hs
-  build-depends: base == 4.*
-               , matrix
+cabal-version: 3.0+name: matrix+version: 0.3.6.2+author: Daniel Casanueva+category: Math+build-type: Simple+license: BSD-3-Clause+license-file: license+maintainer: Daniel Casanueva (daniel.casanueva `at` proton.me)+bug-reports: https://github.com/Daniel-Diaz/matrix/issues+synopsis: A native implementation of matrix operations.+description:+  Matrix library. Basic operations and some algorithms.+  .+  Usage examples are included in the API reference generated by Haddock.+  .+  If you want to use GSL, BLAS and LAPACK, @hmatrix@ (<http://hackage.haskell.org/package/hmatrix>)+  is the way to go.+extra-doc-files: readme.md, changelog.md++Source-repository head+  type: git+  location: git://github.com/Daniel-Diaz/matrix.git++Library+  default-language: Haskell2010+  build-depends: base >= 4.5 && < 5+               , vector >= 0.10+               , deepseq >= 1.3.0.0 && < 1.6+               , primitive >= 0.5+               , semigroups >= 0.9+               , loop >= 0.2+  exposed-modules: Data.Matrix+  ghc-options: -Wall -O2 -fstatic-argument-transformation++Benchmark matrix-mult+  default-language: Haskell2010+  type: exitcode-stdio-1.0+  hs-source-dirs: bench+  main-is: mult.hs+  build-depends: base == 4.*+               , matrix+               , criterion+  ghc-options: -O2 -fstatic-argument-transformation+++Benchmark matrix-rref+  default-language: Haskell2010+  type: exitcode-stdio-1.0+  hs-source-dirs: bench+  main-is: rref.hs+  build-depends: base == 4.*+               , matrix+               , criterion+  ghc-options: -O2 -fstatic-argument-transformation++Test-Suite matrix-test+  default-language: Haskell2010+  type: exitcode-stdio-1.0+  hs-source-dirs: test+  main-is: Main.hs+  build-depends: base == 4.*+               , matrix+               , tasty+               , QuickCheck+               , tasty-quickcheck++Test-Suite matrix-examples+  default-language: Haskell2010+  type: exitcode-stdio-1.0+  hs-source-dirs: test+  main-is: Examples.hs+  build-depends: base == 4.*+               , matrix

@@ -1,5 +1,6 @@  import Data.Matrix+import Data.Ratio import System.Exit (exitFailure) import System.IO (hFlush,stdout) import Data.Either (isLeft)@@ -140,4 +141,32 @@       , Right $ fromList 3 3 [1,0,0, 0,0,1, 0,1,0] :: Either String (Matrix Rational))   , testExample "inverse (6)" $ isLeft $ inverse $ fromList 2 2 [1,1,2,2]   , testEquality "inverse (7)" ( inverse $ fromList 3 3 [0,0,0, 0,0,0, 0,0,0 :: Double], Left "Attempt to invert a non-invertible matrix")-    ]+  , testEquality "rref (bug #42)"+      ( rref (fromList 9 10+                [-1,1,0,0,0,0,0,0,1,1+                ,1,-1,0,0,0,0,1,0,1,1+                ,1,1,-1,0,0,0,0,1,0,1+                ,0,0,1,-1,0,0,0,0,0,1+                ,0,0,0,1,-1,1,0,0,0,1+                ,0,0,0,0,0,-1,1,0,0,1+                ,0,0,0,0,0,0,-1,0,1,1+                ,0,0,0,0,0,0,0,-1,0,1+                ,0,0,0,0,0,0,0,0,-1,1+                ] :: Matrix Rational)+      , Right (fromList 9 10+        [1,0,0,0,(-1)%2,0,0,0,0,0+        ,0,1,0,0,(-1)%2,0,0,0,0,0+        ,0,0,1,0,-1,0,0,0,0,0+        ,0,0,0,1,-1,0,0,0,0,0+        ,0,0,0,0,0,1,0,0,0,0+        ,0,0,0,0,0,0,1,0,0,0+        ,0,0,0,0,0,0,0,1,0,0+        ,0,0,0,0,0,0,0,0,1,0+        ,0,0,0,0,0,0,0,0,0,1+        ] :: Matrix Rational)+      )+  , testEquality "rref (bug #52)"+      ( rref $ fromList 2 3 [1,2,3,2,4,6]+      , Right $ fromList 2 3 [1,2,3,0,0,0]+      )+  ]