l.a. algorithms, etc.

1 files changed, 77 insertions, 75 deletions

diff --git a/examples/tests.hs b/examples/tests.hs
index 37800cd..9de9339 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -16,12 +16,26 @@ import GSL.Fourier
 import GSL.Polynomials
 import LAPACK
 import Test.QuickCheck
-import Test.HUnit
+import Test.HUnit hiding ((~:))
 import Complex
+import LinearAlgebra.Algorithms
+import GSL.Matrix
+import Data.Packed.Instances hiding ((<>))
+
+dist :: (Normed t, Num t) => t -> t -> Double
+dist a b = norm (a-b)
+
+infixl 4 |~|
+a |~| b = a :~8~: b
+
+data Aprox a = (:~) a Int
+
+(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
+a :~n~: b = dist a b < 10^^(-n)
 
 
 {-
--- Bravo por quickCheck!    
+-- Bravo por quickCheck!
 
 pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m  ~~ ident (cols m)
     where infix 2 ~~
@@ -32,7 +46,7 @@ pinvProp2 tol m = 0 < r && r <= c ==> (r==c) `trivial` (m <> pinv m <> m  ~~ m)
           c = cols m
           infix 2 ~~
           (~~) = approxEqual tol 
-        
+
 nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros
     where nt    = trans (nullspace m)
           cr    = corank m
@@ -49,31 +63,28 @@ 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))
+          epsTol > 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)
 
 normVR a b = toScalarR AbsSum (vectorZipR Sub a b)
 
-a |~| b = rows a == rows b && cols a == cols b && eps > normVR (t a) (t b)
+a |~| b = rows a == rows b && cols a == cols b && epsTol > normVR (t a) (t b)
     where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat
 
 (|~~|) = aprox magnitude
 
 v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2
 
-u ~|~ v = normVR u v < eps
-
+u ~|~ v = normVR u v < epsTol
+-}
 
-eps = 1E-8::Double
+epsTol = 1E-8::Double
 
 asFortran m = (rows m >|< cols m) $ toList (fdat m)
 asC m = (rows m >< cols m) $ toList (cdat m)
@@ -81,6 +92,9 @@ asC m = (rows m >< cols m) $ toList (cdat m)
 mulC a b = multiply RowMajor a b
 mulF a b = multiply ColumnMajor a b
 
+infixl 7 <>
+a <> b = mulF a b
+
 cc = mulC ac bf
 cf = mulF af bc
 
@@ -133,14 +147,14 @@ data Sym a = Sym (Matrix a) deriving Show
 instance (Field a, Arbitrary a, Num a) => Arbitrary (Sym a) where
     arbitrary = do
         SqM m <- arbitrary
-        return $ Sym (m `addM` trans m)
+        return $ Sym (m + trans m)
     coarbitrary = undefined
 
 data Her = Her (Matrix (Complex Double)) deriving Show
 instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where
     arbitrary = do
         SqM m <- arbitrary
-        return $ Her (m `addM` (liftMatrix conj) (trans m))
+        return $ Her (m + conjTrans m)
     coarbitrary = undefined
 
 data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show
@@ -171,49 +185,40 @@ instance (Field a, Arbitrary a) => Arbitrary (PairV a) where
    coarbitrary = undefined
 
 
-addM m1 m2 = liftMatrix2 add m1 m2
-
 
 type BaseType = Double
 
-svdTestR fun prod m = u <> s <> trans v |~| m
-                  && u <> trans u |~| ident (rows m)
-                  && v <> trans v |~| ident (cols m)
+svdTestR fun 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 prod m = u <> s' <> (trans v) |~~| m
-                  && u <> (liftMatrix conj) (trans u) |~~| ident (rows m)
-                  && v <> (liftMatrix conj) (trans v) |~~| ident (cols m)
+svdTestC m = u <> s' <> (trans v) |~| m
+                  && u <> conjTrans u |~| ident (rows m)
+                  && v <> conjTrans v |~| ident (cols m)
     where (u,s,v) = svdC m
-          (<>) = prod
           s' = liftMatrix comp s
 
-eigTestC prod (SqM m) = (m <> v) |~~| (v <> diag s)
-                        && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant 1 (rows m) --normalized
+--svdg' m = (u,s',v) where 
+
+eigTestC (SqM m) = (m <> v) |~| (v <> diag s)
+                        && takeDiag (conjTrans v <> v) |~| constant 1 (rows m) --normalized
     where (s,v) = eigC m
-          (<>) = prod
 
-eigTestR prod (SqM m) = (liftMatrix comp m <> v) |~~| (v <> diag s)
-                      -- && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant 1 (rows m) --normalized ???
+eigTestR (SqM m) = (liftMatrix comp m <> v) |~| (v <> diag s)
+                      -- && takeDiag ((liftMatrix conj (trans v)) <> v) |~| constant 1 (rows m) --normalized ???
     where (s,v) = eigR m
-          (<>) = prod
 
-eigTestS prod (Sym m) = (m <> v) |~| (v <> diag s)
+eigTestS (Sym m) = (m <> v) |~| (v <> diag s)
                        && v <> trans v |~| ident (cols m)
     where (s,v) = eigS m
-          (<>) = prod
 
-eigTestH prod (Her m) = (m <> v) |~~| (v <> diag (comp s))
-                        && v <> (liftMatrix conj) (trans v) |~~| ident (cols m)
+eigTestH (Her m) = (m <> v) |~| (v <> diag (comp s))
+                        && v <> conjTrans v |~| ident (cols m)
     where (s,v) = eigH m
-          (<>) = prod
-
-linearSolveSQTest fun eqfun singu prod (PairSM a b) = singu a || (a <> fun a b) ==== b
-    where (<>) = prod
-          (====) = eqfun
 
+linearSolveSQTest fun singu (PairSM a b) = singu a || (a <> fun a b) |~| b
 
 prec = 1E-15
 
@@ -237,8 +242,7 @@ identC n = toComplex(ident n, (0::Double) <>ident n)
 
 --------------------------------------------------------------------
 
-pinvTest f feq m = (m <> f m <> m) `feq` m
-    where (<>) = mulF
+pinvTest f m = (m <> f m <> m) |~| m
 
 pinvR m = linearSolveLSR m (ident (rows m))
 pinvC m = linearSolveLSC m (ident (rows m))
@@ -252,7 +256,7 @@ pinvSVDC m = linearSolveSVDC Nothing m (ident (rows m))
 polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
 
 polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
-    where l1 |~~| l2 = eps > aproxL magnitude l1 l2
+
 
 polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4])
 
@@ -267,17 +271,17 @@ quad2 f a b g1 g2 = quad h a b
 volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) 
                         0 r (const 0) (\x->sqrt (r*r-x*x))
 
-integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < eps)
+integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol)
 
 
 ---------------------------------------------------------------------
 
 arit1 u = vectorMapValR PowVS 2 (vectorMapR Sin u)
           `add` vectorMapValR PowVS 2 (vectorMapR Cos u)
-          ~|~ constant 1 (dim u)
+          |~| constant 1 (dim u)
 
 arit2 u = (vectorMapR Cos u) `mul` (vectorMapR Tan u)
-          ~|~ vectorMapR Sin u
+          |~| vectorMapR Sin u
 
 
 -- arit3 (PairV u v) =
@@ -305,50 +309,48 @@ tests = TestList
 
 main = do
     putStrLn "--------- general -----"
-    quickCheck (\(Sym m) -> m |=| (trans m:: Matrix BaseType))    
+    quickCheck (\(Sym m) -> m == (trans m:: Matrix BaseType))    
     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 $ \m -> m == asC (m :: Matrix BaseType)
+    quickCheck $ \m -> m == asFortran (m :: Matrix BaseType)
+    quickCheck $ \m -> m == (asC . asFortran) (m :: Matrix BaseType)
     putStrLn "--------- MULTIPLY ----"
-    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 $ \(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)
     putStrLn "--------- SVD ---------"
-    quickCheck (svdTestR svdR mulC)
-    quickCheck (svdTestR svdR mulF)
-    quickCheck (svdTestR svdRdd mulC)
-    quickCheck (svdTestR svdRdd mulF)
-    quickCheck (svdTestC mulC)
-    quickCheck (svdTestC mulF)
+    quickCheck (svdTestR svdR)
+    quickCheck (svdTestR svdRdd)
+--    quickCheck (svdTestR svdg)
+    quickCheck svdTestC
     putStrLn "--------- EIG ---------"
-    quickCheck (eigTestC mulC)
-    quickCheck (eigTestC mulF)
-    quickCheck (eigTestR mulC)
-    quickCheck (eigTestR mulF)
-    quickCheck (eigTestS mulC)
-    quickCheck (eigTestS mulF)
-    quickCheck (eigTestH mulC)
-    quickCheck (eigTestH mulF)
+    quickCheck eigTestC
+    quickCheck eigTestR
+    quickCheck eigTestS
+    quickCheck eigTestH
     putStrLn "--------- SOLVE ---------"
-    quickCheck (linearSolveSQTest linearSolveR (|~|) (singular svdR') mulC)
-    quickCheck (linearSolveSQTest linearSolveC (|~~|) (singular svdC') mulF)
-    quickCheck (pinvTest pinvR (|~|))
-    quickCheck (pinvTest pinvC (|~~|))
-    quickCheck (pinvTest pinvSVDR (|~|))
-    quickCheck (pinvTest pinvSVDC (|~~|))
+    quickCheck (linearSolveSQTest linearSolveR (singular svdR'))
+    quickCheck (linearSolveSQTest linearSolveC (singular svdC'))
+    quickCheck (pinvTest pinvR)
+    quickCheck (pinvTest pinvC)
+    quickCheck (pinvTest pinvSVDR)
+    quickCheck (pinvTest pinvSVDC)
     putStrLn "--------- VEC OPER ------"
     quickCheck arit1
     quickCheck arit2
     putStrLn "--------- GSL ------"
     runTestTT tests
-    quickCheck $ \v -> ifft (fft v) ~~ v
+    quickCheck $ \v -> ifft (fft v) |~| v
 
 kk = (2><2)
  [  1.0, 0.0
  , -1.5, 1.0 ::Double]
 
-v = 11 # [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0::Double]
+v = 11 |> [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0::Double]
+
+pol =  [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5]
 
-pol =  [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5]
\ No newline at end of file
+mm = (2><2)
+ [ 0.5, 0.0
+ , 0.0, 0.0 ] :: Matrix Double