diff --git a/sparse-linear-algebra.cabal b/sparse-linear-algebra.cabal index 707e2e9..37d1115 100644 --- a/sparse-linear-algebra.cabal +++ b/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 diff --git a/src/Data/Sparse/Internal/IntM.hs b/src/Data/Sparse/Internal/IntM.hs index 0dc2112..27e051e 100644 --- a/src/Data/Sparse/Internal/IntM.hs +++ b/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 diff --git a/src/Data/Sparse/PPrint.hs b/src/Data/Sparse/PPrint.hs index e49ef90..d576bef 100644 --- a/src/Data/Sparse/PPrint.hs +++ b/src/Data/Sparse/PPrint.hs @@ -47,7 +47,7 @@ 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 @@ -55,7 +55,7 @@ printDN l n = printNpad l n f where -- | 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 = " + " diff --git a/src/Data/Sparse/SpMatrix.hs b/src/Data/Sparse/SpMatrix.hs index 2d31653..53b4234 100644 --- a/src/Data/Sparse/SpMatrix.hs +++ b/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 diff --git a/src/Data/Sparse/SpVector.hs b/src/Data/Sparse/SpVector.hs index ccb05fc..a43f4c9 100644 --- a/src/Data/Sparse/SpVector.hs +++ b/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) diff --git a/src/Numeric/Eps.hs b/src/Numeric/Eps.hs index a957e21..679aedf 100644 --- a/src/Numeric/Eps.hs +++ 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)) diff --git a/src/Numeric/LinearAlgebra/Class.hs b/src/Numeric/LinearAlgebra/Class.hs index 57d9c82..0cbc04e 100644 --- a/src/Numeric/LinearAlgebra/Class.hs +++ b/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} diff --git a/src/Numeric/LinearAlgebra/Sparse.hs b/src/Numeric/LinearAlgebra/Sparse.hs index bf4b5f5..fd0553f 100644 --- a/src/Numeric/LinearAlgebra/Sparse.hs +++ b/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 diff --git a/stack.yaml b/stack.yaml index b454d36..596416c 100644 --- a/stack.yaml +++ b/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. diff --git a/test/Laws.hs b/test/Laws.hs new file mode 100644 index 0000000..250523c --- /dev/null +++ b/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) diff --git a/test/LibSpec.hs b/test/LibSpec.hs index 505e1db..ce39288 100644 --- a/test/LibSpec.hs +++ b/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)