more tests

3 files changed, 163 insertions, 4 deletions

diff --git a/examples/tests.hs b/examples/tests.hs
new file mode 100644
index 0000000..f167b92
--- /dev/null
+++ b/examples/tests.hs
@@ -0,0 +1,149 @@
+--
+-- QuickCheck tests
+--
+
+-----------------------------------------------------------------------------
+
+import Data.Packed.Internal.Vector
+import Data.Packed.Internal.Matrix
+import LAPACK
+import Test.QuickCheck
+import Complex
+
+{-
+-- Bravo por quickCheck!    
+
+pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m  ~~ ident (cols m)
+    where infix 2 ~~
+          (~~) = approxEqual tol 
+
+pinvProp2 tol m = 0 < r && r <= c ==> (r==c) `trivial` (m <> pinv m <> m  ~~ m)
+    where r = rank m
+          c = cols m
+          infix 2 ~~
+          (~~) = approxEqual tol 
+        
+nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros
+    where nt    = trans (nullspace m)
+          cr    = corank m
+          r     = rows m
+          zeros = create [r,cr] $ replicate (r*cr) 0  
+
+-}
+
+r >< c = f where
+    f l | dim v == r*c = matrixFromVector RowMajor c v
+        | otherwise    = error $ "inconsistent list size = "
+                                 ++show (dim v) ++"in ("++show r++"><"++show c++")"
+        where v = fromList l
+
+r >|< c = f where
+    f l | dim v == r*c = matrixFromVector ColumnMajor c v
+        | otherwise    = error $ "inconsistent list size = "
+                                 ++show (dim v) ++"in ("++show r++"><"++show c++")"
+        where v = fromList l
+
+ac = (2><3) [1 .. 6::Double]
+bc = (3><4) [7 .. 18::Double]
+
+mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)]
+
+af = (2>|<3) [1,4,2,5,3,6::Double]
+bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]
+
+a |=| b = rows a == rows b &&
+          cols a == cols b &&
+          toList (cdat a) == toList (cdat b) &&
+          toList (fdat a) == toList (fdat b)
+
+aprox fun a b = rows a == rows b &&
+          cols a == cols b &&
+          eps > aproxL fun (toList (t a)) (toList (t b))
+    where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat
+
+aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length v1)
+
+(|~|) = aprox abs
+(|~~|) = aprox magnitude
+
+eps = 1E-8::Double
+
+asFortran m = (rows m >|< cols m) $ toList (fdat m)
+asC m = (rows m >< cols m) $ toList (cdat m)
+
+mulC a b = multiply RowMajor a b
+mulF a b = multiply ColumnMajor a b
+
+cc = mulC ac bf
+cf = mulF af bc
+
+r = mulC cc (trans cf)
+
+ident n = diag (constant n 1)
+
+rd = (2><2)
+ [  43492.0,  50572.0
+ , 102550.0, 119242.0 :: Double]
+
+instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
+    arbitrary = do
+        r <- arbitrary
+        i <- arbitrary
+        return (r:+i)
+    coarbitrary = undefined
+
+instance (Field a, Arbitrary a) => Arbitrary (Matrix a) where 
+   arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
+                  m <- choose (1,10)
+                  n <- choose (1,10)
+                  l <- vector (m*n)
+                  ctype <- arbitrary
+                  let h = if ctype then (m><n) else (m>|<n)
+                  trMode <- arbitrary
+                  let tr = if trMode then trans else id
+                  return $ tr (h l)
+   coarbitrary = undefined
+
+data PairM a = PairM (Matrix a) (Matrix a) deriving Show
+instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where
+    arbitrary = do
+        a <- choose (1,10)
+        b <- choose (1,10)
+        c <- choose (1,10)
+        l1 <- vector (a*b)
+        l2 <- vector (b*c)
+        return $ PairM ((a><b) (map fromIntegral (l1::[Int]))) ((b><c) (map fromIntegral (l2::[Int])))
+        --return $ PairM ((a><b) l1) ((b><c) l2)
+    coarbitrary = undefined
+
+type BaseType = Double
+
+
+svdTestR fun prod m = u <> s <> trans v |~| m
+                  && u <> trans u |~| ident (rows m)
+                  && v <> trans v |~| ident (cols m)
+    where (u,s,v) = fun m
+          (<>) = prod
+
+
+svdTestC fun prod m = u <> s' <> (trans v) |~~| m
+                  && u <> (liftMatrix conj) (trans u) |~~| ident (rows m)
+                  && v <> (liftMatrix conj) (trans v) |~~| ident (cols m)
+    where (u,s,v) = fun m
+          (<>) = prod
+          s' = liftMatrix comp s
+
+comp v = toComplex (v,constant (dim v) 0)
+
+main = do
+    quickCheck $ \l -> null l || (toList . fromList) l == (l :: [BaseType])
+    quickCheck $ \m -> m |=| asC (m :: Matrix BaseType)
+    quickCheck $ \m -> m |=| asFortran (m :: Matrix BaseType)
+    quickCheck $ \m -> m |=| (asC . asFortran) (m :: Matrix BaseType)
+    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| mulF m1 (m2 :: Matrix BaseType)
+    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| trans (mulF (trans m2) (trans m1 :: Matrix BaseType))
+    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| multiplyG m1 (m2 :: Matrix BaseType)
+    quickCheck (svdTestR svdR mulC)
+    quickCheck (svdTestR svdR mulF)
+    quickCheck (svdTestC svdC mulC)
+    quickCheck (svdTestC svdC mulF)
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 4383e79..0664e9b 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -179,6 +179,15 @@ joinVert ms = case common cols ms of
 joinHoriz :: Field t => [Matrix t] -> Matrix t
 joinHoriz ms = trans. joinVert . map trans $ ms
 
+-- | creates a complex vector from vectors with real and imaginary parts
+toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)
+toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i]
+
+-- | obtains the complex conjugate of a complex vector
+conj :: Vector (Complex Double) -> Vector (Complex Double)
+conj v = asComplex $ cdat $ reshape 2 (asReal v) `mulC` diag (fromList [1,-1])
+    where mulC = multiply RowMajor
+
 ------------------------------------------------------------------------------
 
 -- | Reverse rows 
@@ -191,7 +200,7 @@ fliprl m = fromColumns . reverse . toColumns $ m
 
 -----------------------------------------------------------------
 
-liftMatrix f m = m { dat = f dat, tdat = f tdat } -- check sizes
+liftMatrix f m = m { dat = f (dat m), tdat = f (tdat m) } -- check sizes
 
 ------------------------------------------------------------------
 
@@ -291,7 +300,7 @@ subMatrixG (r0,c0) (rt,ct) x = reshape ct $ fromList $ concat $ map (subList c0
 diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
     m <- createMatrix RowMajor n n
     fun // vec v // mat dat m // check msg [dat m]
-    return m
+    return m {tdat = dat m}
 
 -- | diagonal matrix from a real vector
 diagR :: Vector Double -> Matrix Double
@@ -319,6 +328,6 @@ diagG v = reshape c $ fromList $ [ l!!(i-1) * delta k i | k <- [1..c], i <- [1..
 diagRect s r c
     | dim s < min r c = error "diagRect"
     | r == c    = diag s
-    | r < c     = joinHoriz [diag s , zeros (r,c-r)]
-    | otherwise = joinVert  [diag s , zeros (r-c,c)]
+    | r < c     = trans $ diagRect s c r
+    | r > c     = joinVert  [diag s , zeros (r-c,c)]
     where zeros (r,c) = reshape c $ constant (r*c) 0
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 36d5df7..8f4e6a4 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -145,6 +145,7 @@ asReal v = V { dim = 2*dim v, fptr =  castForeignPtr (fptr v), ptr = castPtr (pt
 asComplex :: Vector Double -> Vector (Complex Double)
 asComplex v = V { dim = dim v `div` 2, fptr =  castForeignPtr (fptr v), ptr = castPtr (ptr v) }
 
+
 constantG n x = fromList (replicate n x)
 
 constantR :: Int -> Double -> Vector Double