new package hmatrix-tests

2 files changed, 523 insertions, 0 deletions

diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
new file mode 100644
index 0000000..647a06c
--- /dev/null
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -0,0 +1,251 @@
+{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-}
+{-# OPTIONS_GHC -fno-warn-unused-imports #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Tests.Instances
+Copyright   :  (c) Alberto Ruiz 2008
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  portable
+
+Arbitrary instances for vectors, matrices.
+
+-}
+
+module Numeric.LinearAlgebra.Tests.Instances(
+    Sq(..),     rSq,cSq,
+    Rot(..),    rRot,cRot,
+    Her(..),    rHer,cHer,
+    WC(..),     rWC,cWC,
+    SqWC(..),   rSqWC, cSqWC,
+    PosDef(..), rPosDef, cPosDef,
+    Consistent(..), rConsist, cConsist,
+    RM,CM, rM,cM,
+    FM,ZM, fM,zM
+) where
+
+import System.Random
+
+import Numeric.LinearAlgebra
+import Control.Monad(replicateM)
+import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
+                      ,sized,classify,Testable,Property
+                      ,quickCheckWith,maxSize,stdArgs,shrink)
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]]
+shrinkListElementwise []     = []
+shrinkListElementwise (x:xs) = [ y:xs | y  <- shrink x                 ]
+                            ++ [ x:ys | ys <- shrinkListElementwise xs ]
+
+shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)]
+shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ]
+#endif
+
+#if MIN_VERSION_QuickCheck(2,1,1)
+#else
+instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
+    arbitrary = do
+        re <- arbitrary
+        im <- arbitrary
+        return (re :+ im)
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+    shrink (re :+ im) = 
+        [ u :+ v | (u,v) <- shrinkPair (re,im) ]
+#else
+    -- this has been moved to the 'Coarbitrary' class in QuickCheck 2
+    coarbitrary = undefined 
+#endif
+
+#endif
+
+chooseDim = sized $ \m -> choose (1,max 1 m)
+
+instance (Field a, Arbitrary a) => Arbitrary (Vector a) where 
+    arbitrary = do m <- chooseDim
+                   l <- vector m
+                   return $ fromList l
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+    -- shrink any one of the components
+    shrink = map fromList . shrinkListElementwise . toList
+
+#else
+    coarbitrary = undefined
+#endif
+
+instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where 
+    arbitrary = do
+        m <- chooseDim
+        n <- chooseDim
+        l <- vector (m*n)
+        return $ (m><n) l
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+    -- shrink any one of the components
+    shrink a = map (rows a >< cols a)
+               . shrinkListElementwise
+               . concat . toLists 
+                     $ a
+#else
+    coarbitrary = undefined
+#endif
+
+
+-- a square matrix
+newtype (Sq a) = Sq (Matrix a) deriving Show
+instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
+    arbitrary = do
+        n <- chooseDim
+        l <- vector (n*n)
+        return $ Sq $ (n><n) l
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+    shrink (Sq a) = [ Sq b | b <- shrink a ]
+#else
+    coarbitrary = undefined
+#endif
+
+
+-- a unitary matrix
+newtype (Rot a) = Rot (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let (q,_) = qr m
+        return (Rot q)
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+#else
+    coarbitrary = undefined
+#endif
+
+
+-- a complex hermitian or real symmetric matrix
+newtype (Her a) = Her (Matrix a) deriving Show
+instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let m' = m/2
+        return $ Her (m' + ctrans m')
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+#else
+    coarbitrary = undefined
+#endif
+
+class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
+instance ArbitraryField Double
+instance ArbitraryField (Complex Double)
+
+
+-- a well-conditioned general matrix (the singular values are between 1 and 100)
+newtype (WC a) = WC (Matrix a) deriving Show
+instance (ArbitraryField a) => Arbitrary (WC a) where
+    arbitrary = do
+        m <- arbitrary
+        let (u,_,v) = svd m
+            r = rows m
+            c = cols m
+            n = min r c
+        sv' <- replicateM n (choose (1,100))
+        let s = diagRect 0 (fromList sv') r c
+        return $ WC (u <> real s <> trans v)
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+#else
+    coarbitrary = undefined
+#endif
+
+
+-- a well-conditioned square matrix (the singular values are between 1 and 100)
+newtype (SqWC a) = SqWC (Matrix a) deriving Show
+instance (ArbitraryField a) => Arbitrary (SqWC a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let (u,_,v) = svd m
+            n = rows m
+        sv' <- replicateM n (choose (1,100))
+        let s = diag (fromList sv')
+        return $ SqWC (u <> real s <> trans v)
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+#else
+    coarbitrary = undefined
+#endif
+
+
+-- a positive definite square matrix (the eigenvalues are between 0 and 100)
+newtype (PosDef a) = PosDef (Matrix a) deriving Show
+instance (ArbitraryField a, Num (Vector a)) 
+    => Arbitrary (PosDef a) where
+    arbitrary = do
+        Her m <- arbitrary
+        let (_,v) = eigSH m
+            n = rows m
+        l <- replicateM n (choose (0,100))
+        let s = diag (fromList l)
+            p = v <> real s <> ctrans v
+        return $ PosDef (0.5 * p + 0.5 * ctrans p)
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+#else
+    coarbitrary = undefined
+#endif
+
+
+-- a pair of matrices that can be multiplied
+newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where
+    arbitrary = do
+        n <- chooseDim
+        k <- chooseDim
+        m <- chooseDim
+        la <- vector (n*k)
+        lb <- vector (k*m)
+        return $ Consistent ((n><k) la, (k><m) lb)
+
+#if MIN_VERSION_QuickCheck(2,0,0)
+    shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ]
+#else
+    coarbitrary = undefined
+#endif
+
+
+
+type RM = Matrix Double
+type CM = Matrix (Complex Double)
+type FM = Matrix Float
+type ZM = Matrix (Complex Float)
+
+
+rM m = m :: RM
+cM m = m :: CM
+fM m = m :: FM
+zM m = m :: ZM
+
+
+rHer (Her m) = m :: RM
+cHer (Her m) = m :: CM
+
+rRot (Rot m) = m :: RM
+cRot (Rot m) = m :: CM
+
+rSq  (Sq m)  = m :: RM
+cSq  (Sq m)  = m :: CM
+
+rWC (WC m) = m :: RM
+cWC (WC m) = m :: CM
+
+rSqWC (SqWC m) = m :: RM
+cSqWC (SqWC m) = m :: CM
+
+rPosDef (PosDef m) = m :: RM
+cPosDef (PosDef m) = m :: CM
+
+rConsist (Consistent (a,b)) = (a,b::RM)
+cConsist (Consistent (a,b)) = (a,b::CM)
+
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
new file mode 100644
index 0000000..c96d3de
--- /dev/null
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -0,0 +1,272 @@
+{-# LANGUAGE CPP, FlexibleContexts #-}
+{-# OPTIONS_GHC -fno-warn-unused-imports #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Tests.Properties
+Copyright   :  (c) Alberto Ruiz 2008
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  portable
+
+Testing properties.
+
+-}
+
+module Numeric.LinearAlgebra.Tests.Properties (
+    dist, (|~|), (~:), Aprox((:~)),
+    zeros, ones,
+    square,
+    unitary,
+    hermitian,
+    wellCond,
+    positiveDefinite,
+    upperTriang,
+    upperHessenberg,
+    luProp,
+    invProp,
+    pinvProp,
+    detProp,
+    nullspaceProp,
+    bugProp,
+    svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
+    svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
+    eigProp, eigSHProp, eigProp2, eigSHProp2,
+    qrProp, rqProp, rqProp1, rqProp2, rqProp3,
+    hessProp,
+    schurProp1, schurProp2,
+    cholProp, exactProp,
+    expmDiagProp,
+    multProp1, multProp2,
+    subProp,
+    linearSolveProp, linearSolveProp2
+) where
+
+import Numeric.LinearAlgebra --hiding (real,complex)
+import Numeric.LinearAlgebra.LAPACK
+import Debug.Trace
+import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
+                      ,sized,classify,Testable,Property
+                      ,quickCheckWith,maxSize,stdArgs,shrink)
+
+trivial :: Testable a => Bool -> a -> Property
+trivial = (`classify` "trivial")
+
+
+-- relative error
+dist :: (Normed c t, Num (c t)) => c t -> c t -> Double
+dist a b = realToFrac r
+    where norm = pnorm Infinity
+          na = norm a
+          nb = norm b
+          nab = norm (a-b)
+          mx = max na nb
+          mn = min na nb
+          r = if mn < peps
+                then mx
+                else nab/mx
+
+infixl 4 |~|
+a |~| b = a :~10~: b
+--a |~| b = dist a b < 10^^(-10)
+
+data Aprox a = (:~) a Int
+-- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
+a :~n~: b = dist a b < 10^^(-n)
+
+------------------------------------------------------
+
+square m = rows m == cols m
+
+-- orthonormal columns
+orthonormal m = ctrans m <> m |~| ident (cols m)
+
+unitary m = square m && orthonormal m
+
+hermitian m = square m && m |~| ctrans m
+
+wellCond m = rcond m > 1/100
+
+positiveDefinite m = minimum (toList e) > 0
+    where (e,_v) = eigSH m
+
+upperTriang m = rows m == 1 || down == z
+    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
+upperHessenberg m = rows m < 3 || down == z
+    where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
+zeros (r,c) = reshape c (constant 0 (r*c))
+
+ones (r,c) = zeros (r,c) + 1
+
+-----------------------------------------------------
+
+luProp m = m |~| p <> l <> u && f (det p) |~| f s
+    where (l,u,p,s) = lu m
+          f x = fromList [x]
+
+invProp m = m <> inv m |~| ident (rows m)
+
+pinvProp m =  m <> p <> m |~| m
+           && p <> m <> p |~| p
+           && hermitian (m<>p)
+           && hermitian (p<>m)
+    where p = pinv m
+
+detProp m = s d1 |~| s d2
+    where d1 = det m
+          d2 = det' * det q
+          det' = product $ toList $ takeDiag r
+          (q,r) = qr m
+          s x = fromList [x]
+
+nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)
+                                     && orthonormal (fromColumns nl))
+    where nl = nullspacePrec 1 m
+          n = fromColumns nl
+          r = rows m
+          c = cols m - rank m
+
+------------------------------------------------------------------
+
+-- testcase for nonempty fpu stack
+-- uncommenting unitary' signature eliminates the problem
+bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v
+    where (u,d,v) = fullSVD m
+          -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool
+          unitary' a = unitary a
+
+------------------------------------------------------------------
+
+-- fullSVD
+svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v
+    where (u,d,v) = fullSVD m
+
+svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where
+    (u,s,v) = svdfun m
+    d = diagRect 0 s (rows m) (cols m)
+
+svdProp1b svdfun m = unitary u && unitary v where
+    (u,_,v) = svdfun m
+
+-- thinSVD
+svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m)
+    where (u,s,v) = thinSVDfun m
+
+-- compactSVD
+svdProp3 m = (m |~| u <> real (diag s) <> trans v
+             && orthonormal u && orthonormal v)
+    where (u,s,v) = compactSVD m
+
+svdProp4 m' = m |~| u <> real (diag s) <> trans v
+           && orthonormal u && orthonormal v
+           && (dim s == r || r == 0 && dim s == 1)
+    where (u,s,v) = compactSVD m
+          m = fromBlocks [[m'],[m']]
+          r = rank m'
+
+svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where
+    s1       = svR  m
+    s2       = svRd m
+    (_,s3,_) = svdR m
+    (_,s4,_) = svdRd m
+    (_,s5,_) = thinSVDR m
+    (_,s6,_) = thinSVDRd m
+
+svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where
+    s1       = svC  m
+    s2       = svCd m
+    (_,s3,_) = svdC m
+    (_,s4,_) = svdCd m
+    (_,s5,_) = thinSVDC m
+    (_,s6,_) = thinSVDCd m
+
+svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
+    where (u,s,v) = svdR m
+          (s',v') = rightSVR m
+          (u',s'') = leftSVR m
+
+svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
+    where (u,s,v) = svdC m
+          (s',v') = rightSVC m
+          (u',s'') = leftSVC m
+
+svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
+    where (u,s,v) = svd m
+          (s',v') = rightSV m
+          (u',_s'') = leftSV m
+          s''' = singularValues m
+
+------------------------------------------------------------------
+
+eigProp m = complex m <> v |~| v <> diag s
+    where (s, v) = eig m
+
+eigSHProp m = m <> v |~| v <> real (diag s)
+              && unitary v
+              && m |~| v <> real (diag s) <> ctrans v
+    where (s, v) = eigSH m
+
+eigProp2 m = fst (eig m) |~| eigenvalues m
+
+eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
+
+------------------------------------------------------------------
+
+qrProp m = q <> r |~| m && unitary q && upperTriang r
+    where (q,r) = qr m
+
+rqProp m = r <> q |~| m && unitary q && upperTriang' r
+    where (r,q) = rq m
+
+rqProp1 m = r <> q |~| m
+    where (r,q) = rq m
+
+rqProp2 m = unitary q
+    where (_r,q) = rq m
+
+rqProp3 m = upperTriang' r
+    where (r,_q) = rq m
+
+upperTriang' r = upptr (rows r) (cols r) * r |~| r
+    where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1
+              where t = f-c
+
+hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
+    where (p,h) = hess m
+
+schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s
+    where (u,s) = schur m
+
+schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme
+    where (u,s) = schur m
+
+cholProp m = m |~| ctrans c <> c && upperTriang c
+    where c = chol m
+
+exactProp m = chol m == chol (m+0)
+
+expmDiagProp m = expm (logm m) :~ 7 ~: complex m
+    where logm = matFunc log
+
+-- reference multiply
+mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
+    where doth u v = sum $ zipWith (*) (toList u) (toList v)
+
+multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)
+
+multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a)
+
+linearSolveProp f m = f m m |~| ident (rows m)
+
+linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
+    where q = min (rows a) (cols a)
+          b = a <> x
+          wc = rank a == q
+
+subProp m = m == (trans . fromColumns . toRows) m
+