sparse-tensor 0.2 → 0.2.1
raw patch · 8 files changed
+117/−76 lines, 8 filesbuild-type:Customsetup-changed
Files
- CHANGELOG.md +7/−2
- Setup.hs +31/−0
- cabal.project +0/−4
- sparse-tensor.cabal +8/−5
- src/Math/Tensor.hs +16/−14
- src/Math/Tensor/Internal/LinearAlgebra.hs +45/−43
- stack.yaml +0/−6
- test/LinearAlgebra.hs +10/−2
CHANGELOG.md view
@@ -1,8 +1,10 @@ # Changelog -## [0.1.0.0] - 2019-08-16-Initial release. +## [0.2.1] - 2019-08-22+- fixed haddock build for ghc 8.6+- fixed numerical instability in gaussian elimination algorithm+ ## [0.2] - 2019-08-21 - removed singletons from dependencies - removed eigen from dependencies@@ -17,3 +19,6 @@ error when using haddock with ghc 8.6 - ensured backwards compatibility down to ghc 8.0.2 (stackage lts-9.21)++## [0.1.0.0] - 2019-08-16+Initial release.
+ Setup.hs view
@@ -0,0 +1,31 @@+{-+Disable some errors and warnings during the haddock pass+ (caused by compiler plugins and hs-boot)+ -}+{-# OPTIONS_GHC -Wall #-}+module Main (main) where++import Distribution.Simple+import Distribution.Simple.Setup++main :: IO ()+main = defaultMainWithHooks simpleUserHooks+ { confHook = \a -> confHook simpleUserHooks a . tweakFlags }++tweakFlags :: ConfigFlags -> ConfigFlags+tweakFlags flags = flags+ { configProgramArgs = addHaddockArgs (configProgramArgs flags) }++addHaddockArgs :: [(String, [String])] -> [(String, [String])]+addHaddockArgs []+ = [("haddock", newHaddockGhcArgs)]+addHaddockArgs (("haddock", args):otherProgsArgs)+ = ("haddock", args ++ newHaddockGhcArgs) : otherProgsArgs+addHaddockArgs (progArgs:otherProgsArgs)+ = progArgs : addHaddockArgs otherProgsArgs++newHaddockGhcArgs :: [String]+newHaddockGhcArgs =+ [ "--optghc=-fdefer-type-errors"+ , "--optghc=-fno-warn-deferred-type-errors"+ ]
− cabal.project
@@ -1,4 +0,0 @@-packages: ./--package sparse-tensor- haddock-options: --optghc=-fdefer-type-errors
sparse-tensor.cabal view
@@ -1,5 +1,5 @@ name: sparse-tensor-version: 0.2+version: 0.2.1 synopsis: typesafe tensor algebra library description: .@@ -17,13 +17,16 @@ maintainer: tobi.reinhart@fau.de, nils.alex@fau.de copyright: 2019 Tobias Reinhart and Nils Alex category: Data, Math, Algebra-build-type: Simple-cabal-version: >=1.10+build-type: Custom+cabal-version: 1.24 extra-source-files: README.md CHANGELOG.md- stack.yaml- cabal.project++custom-setup+ setup-depends:+ base >= 4.9 && < 5,+ Cabal >= 1.24 && < 3.0 source-repository head type: git
src/Math/Tensor.hs view
@@ -387,24 +387,26 @@ isZero :: forall n. KnownNat n => Proxy n -> IsZero n isZero _ = case sameNat (Proxy @n) (Proxy @0)- of Nothing -> unsafeCoerce (NonZero @1)- Just Refl -> Zero+ of Nothing -> unsafeCoerce (NonZero @1)+ Just Refl -> Zero -fromList' :: forall n. KnownNat n => Proxy n -> forall (a :: *). [a] -> Maybe (IndList n a)-fromList' (n :: Proxy n') xs = case isZero n- of Zero -> case xs- of [] -> Just Empty- _ -> Nothing- NonZero -> case xs- of [] -> Nothing- x:xs' -> case fromList' (undefined :: (Proxy (n' - 1))) xs'- of Just v -> Just (x `Append` v)- Nothing -> Nothing+isZero' :: forall n. KnownNat n => IsZero n+isZero' = case sameNat (Proxy @n) (Proxy @0)+ of Nothing -> unsafeCoerce (NonZero @1)+ Just Refl -> Zero -- | Construction of a length typed @'IndList'@ from an untyped list. -fromList :: forall n.KnownNat n => forall (a :: *). [a] -> Maybe (IndList n a)-fromList = fromList' (undefined :: Proxy n)+fromList :: forall n. KnownNat n => forall (a :: *). [a] -> Maybe (IndList n a)+fromList xs = case isZero' @n+ of Zero -> case xs+ of [] -> Just Empty+ _ -> Nothing+ NonZero -> case xs+ of [] -> Nothing+ x:xs' -> case fromList xs'+ of Just v -> Just (x `Append` v)+ Nothing -> Nothing -- | Construction of a length typed @'IndList'@ (partial function).
src/Math/Tensor/Internal/LinearAlgebra.hs view
@@ -47,82 +47,81 @@ pivotsU mat = go (0,0) where go (i,j)- = case findPivot mat (i,j) of+ = case findPivot mat e (i,j) of Nothing -> [] Just (i', j') -> j' : go (i'+1, j'+1)+ maxAbs = maximum $ map (maximum . map abs) $ toLists mat+ e = eps * maxAbs --- naive check for numerical zero--eps :: Double -> Bool-eps = (< 1e-12) . abs+eps :: Double+eps = 1e-12 -- find next pivot in upper triangular matrix -findPivot :: Matrix Double -> (Int, Int) -> Maybe (Int, Int)-findPivot mat (i, j)+findPivot :: Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int)+findPivot mat e (i, j) | n == j = Nothing | m == i = Nothing | otherwise = case nonZeros of [] -> if n == j+1 then Nothing- else findPivot mat (i, j+1)+ else findPivot mat e (i, j+1) (pi, pj):_ -> Just (pi, pj+j) where m = rows mat n = cols mat col = mat ¿ [j]- nonZeros = filter (\(i', _) -> i' >= i) $ find (not . eps) col+ nonZeros = filter (\(i', _) -> i' >= i) $ find (not . (< e) . abs) col --- | Find pivot element below position (i, j) with greatest absolute value.+-- | Find pivot element below position (i, j) with greatest absolute value in the ST monad. -findPivotMax :: Matrix Double -> Int -> Int -> Maybe (Int, Int)-findPivotMax mat i j- | n == j = Nothing- | m == i = Nothing- | otherwise = case nonZeros of- [] -> if n == j+1- then Nothing- else findPivotMax mat i (j+1)- _ -> Just (pi, pj+j)- where- m = rows mat- n = cols mat- col = mat ¿ [j]- nonZeros = filter (\(i', _) -> i' >= i) $ find (not . eps) col- (pi, pj) = maximumBy (\ix jx -> abs (col `atIndex` ix)- `compare`- abs (col `atIndex` jx))- nonZeros+findPivotMax :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s (Maybe (Int, Int))+findPivotMax m n i j mat+ | n == j = return Nothing+ | m == i = return Nothing+ | otherwise =+ do+ col <- mapM (\i' -> do+ x <- readMatrix mat i' j+ return (i', abs x))+ [i..m-1]+ let nonZeros = filter (not . (<eps) . abs . snd) col+ let (pi, _) = maximumBy (\(_, x) (_, y) -> x `compare` y) nonZeros+ case nonZeros of+ [] -> if n == j+1+ then return Nothing+ else findPivotMax m n i (j+1) mat+ _ -> return $ Just (pi, j) -- gaussian elimination of sub matrix below position (i, j) -gaussian' :: Int -> Int -> STMatrix s Double -> ST s ()-gaussian' i j mat = do- m <- liftSTMatrix rows mat- n <- liftSTMatrix cols mat- iPivot' <- liftSTMatrix (\x -> findPivotMax x i j) mat+gaussian' :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()+gaussian' m n i j mat = do+ iPivot' <- findPivotMax m n i j mat case iPivot' of Nothing -> return () Just (r, p) -> do rowOper (SWAP i r (FromCol j)) mat- pv <- liftSTMatrix (`atIndex` (i, p)) mat+ pv <- readMatrix mat i p mapM_ (reduce pv p) [i+1 .. m-1]- gaussian' (i+1) (p+1) mat+ gaussian' m n (i+1) (p+1) mat where reduce pv p r = do- rv <- liftSTMatrix (`atIndex` (r, p)) mat- if eps rv+ rv <- readMatrix mat r p+ if abs rv < eps then return () else let frac = -rv / pv op = AXPY frac i r (FromCol p)- in rowOper op mat+ in do+ rowOper op mat+ mapM_ (\j' -> modifyMatrix mat r j' (\x -> if abs x < eps then 0 else x)) [p..n-1] -- | Gaussian elimination perfomed in-place in the @'ST'@ monad. -gaussianST :: STMatrix s Double -> ST s ()-gaussianST = gaussian' 0 0+gaussianST :: Int -> Int -> STMatrix s Double -> ST s ()+gaussianST m n = gaussian' m n 0 0 -- | Gaussian elimination as pure function. Involves a copy of the input matrix.@@ -144,9 +143,12 @@ gaussian :: Matrix Double -> Matrix Double gaussian mat = runST $ do- m <- thawMatrix mat- gaussianST m- freezeMatrix m+ matST <- thawMatrix mat+ gaussianST m n matST+ freezeMatrix matST+ where+ m = rows mat+ n = cols mat -- | Returns the indices of a maximal linearly independent subset of the columns -- in the matrix.
− stack.yaml
@@ -1,6 +0,0 @@-resolver: lts-14.1--build:- haddock-arguments:- haddock-args:- - --optghc=-fdefer-type-errors
test/LinearAlgebra.hs view
@@ -13,7 +13,7 @@ import Numeric.LinearAlgebra (rank) import qualified Numeric.LinearAlgebra.Data as Matrix -import Math.Tensor.Internal.LinearAlgebra (independentColumnsMat)+import Math.Tensor.Internal.LinearAlgebra (independentColumns, independentColumnsMat) data SmallInt = S0 | S1 deriving (Show, Ord, Eq, Enum, Bounded) @@ -59,8 +59,16 @@ mat = (rows Matrix.>< cols) xs mat' = independentColumnsMat mat +prop_consec :: Positive Int -> Int -> Bool+prop_consec (Positive dim') start =+ independentColumns mat == [0,1]+ where+ dim = dim' + 100+ mat = (dim Matrix.>< dim) $ map fromIntegral [start..]+ testCase1 = testProperty "prop_smallValues" prop_smallValues testCase2 = testProperty "prop_ints" prop_ints testCase3 = testProperty "prop_doubles" prop_doubles+testCase4 = testProperty "prop_consec" prop_consec -linearAlgebraTest = testGroup "LinearAlgebraTest" [testCase1, testCase2, testCase3]+linearAlgebraTest = testGroup "LinearAlgebraTest" [testCase1, testCase2, testCase3, testCase4]