More tests

sparse-linear-algebra.cabal

@@ -71,7 +71,7 @@ library
                        Control.Iterative
                        Numeric.LinearAlgebra.LinearSolvers.Experimental
                        Numeric.LinearAlgebra.EigenSolvers.Experimental
-  build-depends:       base >= 4.7 && < 5
+  build-depends:       base >= 4.8 && < 5
                      , primitive
                      -- , deepseq
                      , containers
@@ -105,10 +105,10 @@ test-suite spec
   ghc-options:         -Wall
   type:                exitcode-stdio-1.0
   hs-source-dirs:      test
-  other-modules:       LibSpec
+  other-modules:       LibSpec, Laws
   main-is:             Spec.hs
   build-depends:       QuickCheck >= 2.8.2
-                     , base
+                     , base >= 4.8 && < 5
                      , containers
                      -- , criterion == 1.1.4.0
                      , hspec

src/Data/Sparse/Internal/IntM.hs

@@ -2,10 +2,12 @@
 module Data.Sparse.Internal.IntM where
 
 import Data.Sparse.Utils
+import Numeric.Eps
 import Numeric.LinearAlgebra.Class
 
 import GHC.Exts
 import Data.Complex
+import Data.Monoid (All(..))
 
 import qualified Data.IntMap.Strict as IM
 
@@ -100,6 +102,12 @@ instance (Normed a, Magnitude a ~ RealScalar a, RealScalar a ~ Scalar a) => Norm
   norm2  c = sqrt (norm2Sq c)
   norm2' c = sqrt (norm2Sq c)
 
+instance Epsilon a => Epsilon (IntM a) where
+  nearZero = getAll . foldMap (All . nearZero)
+  -- TODO: rewrite with generic merge combinator?
+  near (IntM x) (IntM y) = nearZero (IntM (IM.difference x y))
+                        && nearZero (IntM (IM.difference y x))
+                        && getAll (foldMap All (IM.intersectionWith (\a b -> a `near` b) x y))
 
 
 -- -- | list to IntMap

src/Data/Sparse/PPrint.hs

@@ -47,15 +47,15 @@ prdefC = PPOpts 4 2 16  -- complex values
 
 
 -- | Pretty print an array of real numbers
-printDN :: (PrintfArg a, Epsilon a, Ord a) =>
+printDN :: (PrintfArg a, Epsilon a, Ord a, Floating a) =>
      Int -> Int -> PPrintOptions -> [a] -> String
 printDN l n = printNpad l n f where
   f o x | isNz x = printf (prepD o x) x
         | otherwise = printf "_"
 
 
 -- | Pretty print an array of complex numbers
-printCN :: (PrintfArg a, Epsilon a, Epsilon (Complex a), Ord a) =>
+printCN :: (PrintfArg a, Floating a, Epsilon a, Epsilon (Complex a), Ord a) =>
      Int -> Int -> PPrintOptions -> [Complex a] -> String
 printCN l n = printNpad l n f where
   f o x | nearZero (re x) && isNz (imagPart x) =
@@ -99,7 +99,7 @@ printNpad llen nmax f o@PPOpts{..} xxl = commas [h,l] where
 
 
 -- | printf format string
-prepD :: (Ord t, Epsilon t) => PPrintOptions -> t -> String
+prepD :: (Ord t, Floating t, Epsilon t) => PPrintOptions -> t -> String
 prepD PPOpts{..} x = s where
   s | nearZero x  = "_"
     | abs x > magHi || abs x < magLo = s0 ++ "e"
@@ -115,7 +115,7 @@ prepD PPOpts{..} x = s where
 
 
 -- | printf format string for a Complex
-prepC :: (Epsilon t, Ord t) => PPrintOptions -> Complex t -> String
+prepC :: (Epsilon t, Floating t, Ord t) => PPrintOptions -> Complex t -> String
 prepC opts (r :+ i) = prepD opts r ++ oi where
   oi = concat [s, prepD opts i', "i"]
   s | signum i >= 0 = " + "

src/Data/Sparse/SpMatrix.hs

@@ -419,7 +419,7 @@ isUpperTriSM m = m == lm where
   lm = ifilterSM (\i j _ -> i <= j) m
 
 -- |Is the matrix orthogonal? i.e. Q^t ## Q == I
-isOrthogonalSM :: (Eq a, Epsilon a, MatrixRing (SpMatrix a)) => SpMatrix a -> Bool
+isOrthogonalSM :: (Eq a, Num a, Epsilon a, MatrixRing (SpMatrix a)) => SpMatrix a -> Bool
 isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where
   rsm = roundZeroOneSM $ transpose sm ## sm
 
@@ -650,7 +650,7 @@ sparsifySM (SM d im) = SM d $ sparsifyIM2 im
 
 -- ** Value rounding
 -- | Round almost-0 and almost-1 to 0 and 1 respectively
-roundZeroOneSM :: Epsilon a => SpMatrix a -> SpMatrix a
+roundZeroOneSM :: (Num a, Epsilon a) => SpMatrix a -> SpMatrix a
 roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im  
 
 

src/Data/Sparse/SpVector.hs

@@ -130,6 +130,9 @@ instance (Normed a, Magnitude a ~ RealScalar a, RealScalar a ~ Scalar a) => Norm
   norm2  c = sqrt (norm2Sq c)
   norm2' c = sqrt (norm2Sq c)
 
+instance Epsilon a => Epsilon (SpVector a) where
+  nearZero v = nearZero (svData v)
+  near v w = near (svData v) (svData w)
 
 dotS :: InnerSpace t => SpVector t -> SpVector t -> Scalar (IntM t)
 (SV m a) `dotS` (SV n b)

src/Numeric/Eps.hs

@@ -1,4 +1,4 @@
-{-# language FlexibleInstances #-}
+{-# language FlexibleInstances, DefaultSignatures #-}
 -----------------------------------------------------------------------------
 -- |
 -- Copyright   :  (C) 2016 Marco Zocca, 2012-2015 Edward Kmett
@@ -29,9 +29,13 @@ import Foreign.C.Types (CFloat, CDouble)
 --
 -- >>> nearZero (1e-7 :: Float)
 -- True
-class (Floating a, Num a) => Epsilon a where
+class Epsilon a where
   -- | Determine if a quantity is near zero.
   nearZero :: a -> Bool
+  -- | Determine if two quantities are near.
+  near :: a -> a -> Bool
+  default near :: Num a => a -> a -> Bool
+  near x y = nearZero (x - y)
 
 -- | @'abs' a '<=' 1e-6@
 instance Epsilon Float where
@@ -49,6 +53,8 @@ instance Epsilon CFloat where
 instance Epsilon CDouble where
   nearZero a = abs a <= 1e-12
 
+instance Epsilon Rational where
+  nearZero a = a == 0
 
 -- | @'magnitude' a '<=' 1e-6@
 instance Epsilon (Complex Float) where
@@ -72,8 +78,8 @@ instance Epsilon (Complex CDouble) where
 
 
 -- | Is this quantity close to 1 ?
-nearOne :: Epsilon a => a -> Bool
-nearOne x = nearZero (1 - x)
+nearOne :: (Epsilon a, Num a) => a -> Bool
+nearOne = near 1
 
 -- | Is this quantity distinguishable from 0 ?
 isNz :: Epsilon a => a -> Bool
@@ -83,7 +89,7 @@ withDefault :: (t -> Bool) -> t -> t -> t
 withDefault q d x | q x = d
                   | otherwise = x
 
-roundZero, roundOne, roundZeroOne :: Epsilon a => a -> a
+roundZero, roundOne, roundZeroOne :: (Epsilon a, Num a) => a -> a
 roundZero = withDefault nearZero (fromIntegral (0 :: Int))
 roundOne = withDefault nearOne (fromIntegral (1 :: Int))
 

src/Numeric/LinearAlgebra/Class.hs

@@ -121,7 +121,7 @@ hilbertDistSq x y = t <.> t where
 
 -- * Normed vector spaces
 
-class (InnerSpace v, Num (RealScalar v), Eq (RealScalar v), Epsilon (Magnitude v), Show (Magnitude v), Ord (Magnitude v)) => Normed v where
+class (InnerSpace v, Num (RealScalar v), Eq (RealScalar v), Floating (RealScalar v), Floating (Magnitude v), Epsilon (Magnitude v), Show (Magnitude v), Ord (Magnitude v)) => Normed v where
   type Magnitude v :: *
   type RealScalar v :: *
   -- | L1 norm
@@ -391,6 +391,7 @@ instance VectorSpace (t) where {type Scalar (t) = t; (.*) = (*) };
 -- ScalarType(Integer)
 ScalarType(Float)
 ScalarType(Double)
+ScalarType(Rational)
 ScalarType(Complex Float)
 ScalarType(Complex Double)
 -- ScalarType(CSChar)
@@ -407,6 +408,7 @@ ScalarType(Complex Double)
 
 instance InnerSpace Float  where {(<.>) = (*)}
 instance InnerSpace Double where {(<.>) = (*)}
+instance InnerSpace Rational where {(<.>) = (*)}
 instance InnerSpace (Complex Float)  where {x <.> y = x * conjugate y}
 instance InnerSpace (Complex Double) where {x <.> y = x * conjugate y}
 

src/Numeric/LinearAlgebra/Sparse.hs

@@ -598,7 +598,7 @@ tmc6 = fromListDenseSM 2 $ zipWith (:+) [1,2,3,4] [4,3,2,1]
 -- | Given a matrix A, a vector b and a positive integer `n`, this procedure finds the basis of an order `n` Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q.
 -- At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector.
 arnoldi :: (MatrixType (SpVector a) ~ SpMatrix a, V (SpVector a) ,
-            Scalar (SpVector a) ~ a, Epsilon a, MonadThrow m) =>
+            Scalar (SpVector a) ~ a, Floating a, Epsilon a, MonadThrow m) =>
      SpMatrix a                    -- ^ System matrix
      -> SpVector a                 -- ^ r.h.s.
      -> Int                        -- ^ Max. # of iterations

stack.yaml

@@ -17,7 +17,7 @@
 #  location: "./custom-snapshot.yaml"
 
 # resolver: lts-7.2
-resolver: lts-9.6
+resolver: lts-11.0
 
 # User packages to be built.
 # Various formats can be used as shown in the example below.

test/Laws.hs

@@ -0,0 +1,45 @@
+{-# LANGUAGE PartialTypeSignatures, FlexibleContexts #-}
+{-# GHC_OPTIONS -Wno-partial-type-signatures #-}
+module Laws where
+
+import Test.Hspec
+import Numeric.LinearAlgebra.Class
+
+-- All the axioms of an Abelian group
+associative :: (Eq a, AdditiveGroup a, _) => a -> a -> a -> _
+associative a b c = (a ^+^ b) ^+^ c `shouldBe` a ^+^ (b ^+^ c)
+
+cancellative :: (Eq a, AdditiveGroup a, _) => a -> _
+cancellative a = a ^-^ a `shouldBe` zeroV
+
+commutative :: (Eq a, AdditiveGroup a, _) => a -> a -> _
+commutative a b = a ^+^ b `shouldBe` b ^+^ a
+
+neutralZero :: (Eq a, AdditiveGroup a, _) => a -> _
+neutralZero a = a ^+^ zeroV `shouldBe` a
+
+-- All the axioms of a module (action is associative and bilinear)
+associativeScalar :: (Eq v, VectorSpace v, _) => Scalar v -> Scalar v -> v -> _
+associativeScalar a b v = (a * b) .* v `shouldBe` a .* (b .* v)
+
+neutralScalar :: (Eq v, VectorSpace v, _) => v -> _
+neutralScalar v = 1 .* v `shouldBe` v
+
+distributiveScalar :: (Eq v, VectorSpace v, _) => Scalar v -> Scalar v -> v -> _
+distributiveScalar a b v = (a + b) .* v `shouldBe` a .* v ^+^ b .* v
+
+distributiveScalar2 :: (Eq v, VectorSpace v, _) => Scalar v -> v -> v -> _
+distributiveScalar2 a v w = a .* (v ^+^ w) `shouldBe` a .* v ^+^ a .* w
+
+-- Inner product should be bilinear, commutative in the real case and positive
+innerProductBilinear :: (Eq (Scalar v), InnerSpace v, _) => v -> v -> v -> _
+innerProductBilinear a b c = (a ^+^ b) <.> c `shouldBe` (a <.> c) + (b <.> c)
+
+innerProductBilinear2 :: (Eq (Scalar v), InnerSpace v, _) => Scalar v -> v -> v -> _
+innerProductBilinear2 l b c = (l .* b) <.> c `shouldBe` l * (b <.> c)
+
+innerProductCommutative :: (Eq (Scalar v), InnerSpace v, _) => v -> v -> _
+innerProductCommutative v w = v <.> w `shouldBe` w <.> v
+
+innerProductPositive :: (Ord (Scalar v), InnerSpace v, _) => v -> _
+innerProductPositive v = (v <.> v) `shouldSatisfy` (>= 0)

test/LibSpec.hs

@@ -1,6 +1,7 @@
 {-# LANGUAGE FlexibleContexts, TypeFamilies #-}
 {-# language ScopedTypeVariables, FlexibleInstances #-}
-{-# OPTIONS_GHC -Wno-missing-signatures #-}
+{-# LANGUAGE TypeApplications #-}
+{-# OPTIONS_GHC -Wno-missing-signatures -Wno-orphans #-}
 -----------------------------------------------------------------------------
 -- |
 -- Copyright   :  (C) 2016 Marco Zocca
@@ -26,15 +27,39 @@ import Test.Hspec
 import Test.Hspec.QuickCheck
 import Test.QuickCheck
 
+import Laws
 
 main :: IO ()
 main = hspec spec
 
 spec :: Spec
 spec = do
+  describe "Additive group instance for SpVector Rational" $ do
+    prop "Addition is commutative"     $ commutative @(SpVector Rational)
+    prop "Subtraction is cancellative" $ cancellative @(SpVector Rational)
+    prop "Zero is neutral"             $ neutralZero @(SpVector Rational)
+    prop "Addition is associative"     $ associative @(SpVector Rational)
+  describe "Vector space instance for SpVector Rational" $ do
+    prop "Scalar multiplication is associative"    $ associativeScalar @(SpVector Rational)
+    prop "Scalar multiplication is unital"         $ neutralScalar @(SpVector Rational)
+    prop "Scalar multiplication is distributive"   $ distributiveScalar @(SpVector Rational)
+    prop "Scalar multiplication is distributive 2" $ distributiveScalar2 @(SpVector Rational)
+  describe "Inner space instance for SpVector Rational" $ do
+    prop "Inner product is linear in first arg"  $ innerProductBilinear @(SpVector Rational)
+    prop "Inner product is linear in first arg2" $ innerProductBilinear2 @(SpVector Rational)
+    prop "Inner product is commutative"          $ innerProductCommutative @(SpVector Rational)
+    prop "Inner product is semidefinite"         $ innerProductPositive @(SpVector Rational)
+  describe "Additive group instance for SpMatrix Rational" $ do
+    prop "Addition is commutative"     $ commutative @(SpMatrix Rational)
+    prop "Subtraction is cancellative" $ cancellative @(SpMatrix Rational)
+    prop "Zero is neutral"             $ neutralZero @(SpMatrix Rational)
+    prop "Addition is associative"     $ associative @(SpMatrix Rational)
+  describe "Vector space instance for SpMatrix Rational" $ do
+    prop "Scalar multiplication is associative"    $ associativeScalar @(SpMatrix Rational)
+    prop "Scalar multiplication is unital"         $ neutralScalar @(SpMatrix Rational)
+    prop "Scalar multiplication is distributive"   $ distributiveScalar @(SpMatrix Rational)
+    prop "Scalar multiplication is distributive 2" $ distributiveScalar2 @(SpMatrix Rational)
   describe "Numeric.LinearAlgebra.Sparse : Library" $ do
-    prop "Subtraction is cancellative" $ \(x :: SpVector Double) ->
-      norm2Sq (x ^-^ x) `shouldBe` zeroV
     it "<.> : inner product (Real)" $
       tv0 <.> tv0 `shouldBe` 61
     it "<.> : inner product (Complex)" $
@@ -375,7 +400,7 @@ checkTriLowerSolve lmat rhs = do
 
 checkArnoldi :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,
       Normed (SpVector t), MatrixRing (SpMatrix t),
-      LinearVectorSpace (SpVector t), Epsilon t, MonadThrow m) =>
+      LinearVectorSpace (SpVector t), Floating t, Epsilon t, MonadThrow m) =>
      SpMatrix t -> Int -> m Bool
 checkArnoldi aa kn = do -- nearZero (normFrobenius $ lhs ^-^ rhs) where
   let b = onesSV (nrows aa)
@@ -522,9 +547,13 @@ genSpVDense n = do
 
 
 -- | An Arbitrary SpVector such that at least one entry is nonzero
-instance Arbitrary (SpVector Double) where
+instance (Arbitrary a, Epsilon a) => Arbitrary (SpVector a) where
   arbitrary = sized genSpV `suchThat` any isNz 
 
+instance (Arbitrary a, Epsilon a) => Arbitrary (SpMatrix a) where
+  arbitrary = sized2 $ \n m -> do
+                d <- choose (0, n*m)
+                genSpM0 n m d
 
 -- | An arbitrary square SpMatrix
 newtype PropMat0 a = PropMat0 (SpMatrix a) deriving (Eq, Show)