summaryrefslogtreecommitdiff
path: root/examples
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2007-09-21 18:10:59 +0000
committerAlberto Ruiz <aruiz@um.es>2007-09-21 18:10:59 +0000
commitedfaf9e0d1dcfccc9015476510a23e8cf64288be (patch)
tree2b11f0a933a7cce2362aed26ac160312e6b9431a /examples
parent6cafd2f26a89008cc0db02e70e39f92a50ec4b4d (diff)
algorithms refactoring
Diffstat (limited to 'examples')
-rw-r--r--examples/oldtests.hs120
-rw-r--r--examples/tests.hs357
2 files changed, 0 insertions, 477 deletions
diff --git a/examples/oldtests.hs b/examples/oldtests.hs
deleted file mode 100644
index 7d4701c..0000000
--- a/examples/oldtests.hs
+++ /dev/null
@@ -1,120 +0,0 @@
1import Test.HUnit
2import LinearAlgebra
3import GSL hiding (exp)
4import System.Random(randomRs,mkStdGen)
5
6realMatrix = fromLists :: [[Double]] -> Matrix Double
7realVector = fromList :: [Double] -> Vector Double
8
9
10
11infixl 2 =~=
12a =~= b = pnorm PNorm1 (flatten (a - b)) < 1E-6
13
14randomMatrix seed (n,m) = reshape m $ realVector $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed
15
16randomMatrixC seed (n,m) = toComplex (randomMatrix seed (n,m), randomMatrix (seed+1) (n,m))
17
18besselTest = do
19 let (r,e) = bessel_J0_e 5.0
20 let expected = -0.17759677131433830434739701
21 assertBool "bessel_J0_e" ( abs (r-expected) < e )
22
23exponentialTest = do
24 let (v,e,err) = exp_e10_e 30.0
25 let expected = exp 30.0
26 assertBool "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )
27
28disp m = putStrLn (format " " show m)
29
30ms = realMatrix [[1,2,3]
31 ,[-4,1,7]]
32
33ms' = randomMatrix 27 (50,100)
34
35ms'' = toComplex (randomMatrix 100 (50,100),randomMatrix 101 (50,100))
36
37fullsvdTest method mat msg = do
38 let (u,s,vt) = method mat
39 assertBool msg (u <> s <> trans vt =~= mat)
40
41svdg' m = (u, diag s, v) where (u,s,v) = svdg m
42
43full_svd_Rd = svdRdd
44
45--------------------------------------------------------------------
46
47mcu = toComplex (randomMatrix 33 (20,20),randomMatrix 34 (20,20))
48
49mcur = randomMatrix 35 (40,40)
50
51-- eigenvectors are columns
52eigTest method m msg = do
53 let (s,v) = method m
54 assertBool msg $ m <> v =~= v <> diag s
55
56bigmat = m + trans m where m = randomMatrix 18 (1000,1000)
57bigmatc = mc + conjTrans mc where mc = toComplex(m,m)
58 m = randomMatrix 19 (1000,1000)
59
60--------------------------------------------------------------------
61
62invTest msg m = do
63 assertBool msg $ m <> inv m =~= ident (rows m)
64
65invComplexTest msg m = do
66 assertBool msg $ m <> invC m =~= identC (rows m)
67
68invC m = linearSolveC m (identC (rows m))
69
70identC = comp . ident
71
72--------------------------------------------------------------------
73
74pinvTest f msg m = do
75 assertBool msg $ m <> f m <> m =~= m
76
77pinvC m = linearSolveLSC m (identC (rows m))
78
79pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m))
80
81pinvSVDC m = linearSolveSVDC Nothing m (identC (rows m))
82
83--------------------------------------------------------------------
84
85
86tests = TestList [
87 TestCase $ besselTest
88 , TestCase $ exponentialTest
89 , TestCase $ invTest "inv 100x100" (randomMatrix 18 (100,100))
90 , TestCase $ invComplexTest "complex inv 100x100" (randomMatrixC 18 (100,100))
91 , TestCase $ pinvTest (pinvTolg 1) "pinvg 100x50" (randomMatrix 18 (100,50))
92 , TestCase $ pinvTest pinv "pinv 100x50" (randomMatrix 18 (100,50))
93 , TestCase $ pinvTest pinv "pinv 50x100" (randomMatrix 18 (50,100))
94 , TestCase $ pinvTest pinvSVDR "pinvSVDR 100x50" (randomMatrix 18 (100,50))
95 , TestCase $ pinvTest pinvSVDR "pinvSVDR 50x100" (randomMatrix 18 (50,100))
96 , TestCase $ pinvTest pinvC "pinvC 100x50" (randomMatrixC 18 (100,50))
97 , TestCase $ pinvTest pinvC "pinvC 50x100" (randomMatrixC 18 (50,100))
98 , TestCase $ pinvTest pinvSVDC "pinvSVDC 100x50" (randomMatrixC 18 (100,50))
99 , TestCase $ pinvTest pinvSVDC "pinvSVDC 50x100" (randomMatrixC 18 (50,100))
100 , TestCase $ eigTest eigC mcu "eigC"
101 , TestCase $ eigTest eigR mcur "eigR"
102 , TestCase $ eigTest eigS (mcur+trans mcur) "eigS"
103 , TestCase $ eigTest eigSg (mcur+trans mcur) "eigSg"
104 , TestCase $ eigTest eigH (mcu+ (conjTrans) mcu) "eigH"
105 , TestCase $ eigTest eigHg (mcu+ (conjTrans) mcu) "eigHg"
106 , TestCase $ fullsvdTest svdg' ms "GSL svd small"
107 , TestCase $ fullsvdTest svdR ms "fullsvdR small"
108 , TestCase $ fullsvdTest svdR (trans ms) "fullsvdR small"
109 , TestCase $ fullsvdTest svdR ms' "fullsvdR"
110 , TestCase $ fullsvdTest svdR (trans ms') "fullsvdR"
111 , TestCase $ fullsvdTest full_svd_Rd ms' "fullsvdRd"
112 , TestCase $ fullsvdTest full_svd_Rd (trans ms') "fullsvdRd"
113 , TestCase $ fullsvdTest svdC ms'' "fullsvdC"
114 , TestCase $ fullsvdTest svdC (trans ms'') "fullsvdC"
115 , TestCase $ eigTest eigS bigmat "big eigS"
116 , TestCase $ eigTest eigH bigmatc "big eigH"
117 , TestCase $ eigTest eigR bigmat "big eigR"
118 ]
119
120main = runTestTT tests
diff --git a/examples/tests.hs b/examples/tests.hs
deleted file mode 100644
index dcc3cbf..0000000
--- a/examples/tests.hs
+++ /dev/null
@@ -1,357 +0,0 @@
1{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
2
3--
4-- QuickCheck tests
5--
6
7-----------------------------------------------------------------------------
8
9import Data.Packed.Internal((>|<), fdat, cdat, multiply', multiplyG, MatrixOrder(..))
10import GSL hiding (sin,cos,exp,choose)
11import LinearAlgebra hiding ((<>))
12import Test.QuickCheck
13import Test.HUnit hiding ((~:))
14
15
16dist :: (Normed t, Num t) => t -> t -> Double
17dist a b = pnorm Infinity (a-b)
18
19infixl 4 |~|
20a |~| b = a :~8~: b
21
22data Aprox a = (:~) a Int
23
24(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
25a :~n~: b = dist a b < 10^^(-n)
26
27
28{-
29-- Bravo por quickCheck!
30
31pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m ~~ ident (cols m)
32 where infix 2 ~~
33 (~~) = approxEqual tol
34
35pinvProp2 tol m = 0 < r && r <= c ==> (r==c) `trivial` (m <> pinv m <> m ~~ m)
36 where r = rank m
37 c = cols m
38 infix 2 ~~
39 (~~) = approxEqual tol
40
41nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros
42 where nt = trans (nullspace m)
43 cr = corank m
44 r = rows m
45 zeros = create [r,cr] $ replicate (r*cr) 0
46
47-}
48
49ac = (2><3) [1 .. 6::Double]
50bc = (3><4) [7 .. 18::Double]
51
52mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)]
53
54af = (2>|<3) [1,4,2,5,3,6::Double]
55bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]
56
57
58{-
59aprox fun a b = rows a == rows b &&
60 cols a == cols b &&
61 epsTol > aproxL fun (toList (t a)) (toList (t b))
62 where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat
63
64aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length v1)
65
66normVR a b = toScalarR AbsSum (vectorZipR Sub a b)
67
68a |~| b = rows a == rows b && cols a == cols b && epsTol > normVR (t a) (t b)
69 where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat
70
71(|~~|) = aprox magnitude
72
73v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2
74
75u ~|~ v = normVR u v < epsTol
76-}
77
78epsTol = 1E-8::Double
79
80asFortran m = (rows m >|< cols m) $ toList (fdat m)
81asC m = (rows m >< cols m) $ toList (cdat m)
82
83mulC a b = multiply' RowMajor a b
84mulF a b = multiply' ColumnMajor a b
85
86identC = comp . ident
87
88infixl 7 <>
89a <> b = mulF a b
90
91cc = mulC ac bf
92cf = mulF af bc
93
94r = mulC cc (trans cf)
95
96rd = (2><2)
97 [ 27736.0, 65356.0
98 , 65356.0, 154006.0 ::Double]
99
100instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
101 arbitrary = do
102 r <- arbitrary
103 i <- arbitrary
104 return (r:+i)
105 coarbitrary = undefined
106
107instance (Field a, Arbitrary a) => Arbitrary (Matrix a) where
108 arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
109 m <- choose (1,10)
110 n <- choose (1,10)
111 l <- vector (m*n)
112 ctype <- arbitrary
113 let h = if ctype then (m><n) else (m>|<n)
114 trMode <- arbitrary
115 let tr = if trMode then trans else id
116 return $ tr (h l)
117 coarbitrary = undefined
118
119data PairM a = PairM (Matrix a) (Matrix a) deriving Show
120instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where
121 arbitrary = do
122 a <- choose (1,10)
123 b <- choose (1,10)
124 c <- choose (1,10)
125 l1 <- vector (a*b)
126 l2 <- vector (b*c)
127 return $ PairM ((a><b) (map fromIntegral (l1::[Int]))) ((b><c) (map fromIntegral (l2::[Int])))
128 --return $ PairM ((a><b) l1) ((b><c) l2)
129 coarbitrary = undefined
130
131data SqM a = SqM (Matrix a) deriving Show
132instance (Field a, Arbitrary a) => Arbitrary (SqM a) where
133 arbitrary = do
134 n <- choose (1,10)
135 l <- vector (n*n)
136 return $ SqM $ (n><n) l
137 coarbitrary = undefined
138
139data Sym a = Sym (Matrix a) deriving Show
140instance (Linear Vector a, Arbitrary a) => Arbitrary (Sym a) where
141 arbitrary = do
142 SqM m <- arbitrary
143 return $ Sym (m + trans m)
144 coarbitrary = undefined
145
146data Her = Her (Matrix (Complex Double)) deriving Show
147instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where
148 arbitrary = do
149 SqM m <- arbitrary
150 return $ Her (m + conjTrans m)
151 coarbitrary = undefined
152
153data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show
154instance (Num a, Field a, Arbitrary a) => Arbitrary (PairSM a) where
155 arbitrary = do
156 a <- choose (1,10)
157 c <- choose (1,10)
158 l1 <- vector (a*a)
159 l2 <- vector (a*c)
160 return $ PairSM ((a><a) (map fromIntegral (l1::[Int]))) ((a><c) (map fromIntegral (l2::[Int])))
161 --return $ PairSM ((a><a) l1) ((a><c) l2)
162 coarbitrary = undefined
163
164instance (Field a, Arbitrary a) => Arbitrary (Vector a) where
165 arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
166 m <- choose (1,100)
167 l <- vector m
168 return $ fromList l
169 coarbitrary = undefined
170
171data PairV a = PairV (Vector a) (Vector a)
172instance (Field a, Arbitrary a) => Arbitrary (PairV a) where
173 arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
174 m <- choose (1,100)
175 l1 <- vector m
176 l2 <- vector m
177 return $ PairV (fromList l1) (fromList l2)
178 coarbitrary = undefined
179
180
181
182type BaseType = Complex Double
183
184svdTestR fun m = u <> s <> trans v |~| m
185 && u <> trans u |~| ident (rows m)
186 && v <> trans v |~| ident (cols m)
187 where (u,s,v) = fun m
188
189
190svdTestC m = u <> s' <> (trans v) |~| m
191 && u <> conjTrans u |~| identC (rows m)
192 && v <> conjTrans v |~| identC (cols m)
193 where (u,s,v) = svdC m
194 s' = liftMatrix comp s
195
196--svdg' m = (u,s',v) where
197
198eigTestC (SqM m) = (m <> v) |~| (v <> diag s)
199 && takeDiag (conjTrans v <> v) |~| comp (constant 1 (rows m)) --normalized
200 where (s,v) = eigC m
201
202eigTestR (SqM m) = (liftMatrix comp m <> v) |~| (v <> diag s)
203 -- && takeDiag ((liftMatrix conj (trans v)) <> v) |~| constant 1 (rows m) --normalized ???
204 where (s,v) = eigR m
205
206eigTestS (Sym m) = (m <> v) |~| (v <> diag s)
207 && v <> trans v |~| ident (cols m)
208 where (s,v) = eigS m
209
210eigTestH (Her m) = (m <> v) |~| (v <> diag (comp s))
211 && v <> conjTrans v |~| identC (cols m)
212 where (s,v) = eigH m
213
214linearSolveSQTest fun singu (PairSM a b) = singu a || (a <> fun a b) |~| b
215
216prec = 1E-15
217
218singular fun m = s1 < prec || s2/s1 < prec
219 where (_,ss,v) = fun m
220 s = toList ss
221 s1 = maximum s
222 s2 = minimum s
223
224{-
225invTest msg m = do
226 assertBool msg $ m <> inv m =~= ident (rows m)
227
228invComplexTest msg m = do
229 assertBool msg $ m <> invC m =~= identC (rows m)
230
231invC m = linearSolveC m (identC (rows m))
232
233identC n = toComplex(ident n, (0::Double) <>ident n)
234-}
235
236--------------------------------------------------------------------
237
238pinvTest f m = (m <> f m <> m) |~| m
239
240pinvR m = linearSolveLSR m (ident (rows m))
241pinvC m = linearSolveLSC m (identC (rows m))
242
243pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m))
244
245pinvSVDC m = linearSolveSVDC Nothing m (identC (rows m))
246
247--------------------------------------------------------------------
248
249polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
250
251polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
252
253
254polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4])
255
256---------------------------------------------------------------------
257
258quad f a b = fst $ integrateQAGS 1E-9 100 f a b
259
260-- A multiple integral can be easily defined using partial application
261quad2 f a b g1 g2 = quad h a b
262 where h x = quad (f x) (g1 x) (g2 x)
263
264volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y))
265 0 r (const 0) (\x->sqrt (r*r-x*x))
266
267integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol)
268
269
270---------------------------------------------------------------------
271
272arit1 u = sin u ^ 2 + cos u ^ 2 |~| 1
273 where _ = u :: Vector Double
274
275arit2 u = sin u ** 2 + cos u ** 2 |~| 1
276 where _ = u :: Vector Double
277
278arit3 u = cos u * tan u |~| sin u
279 where _ = u :: Vector Double
280
281arit4 u = (cos u * tan u) :~6~: sin u
282 where _ = u :: Vector (Complex Double)
283
284---------------------------------------------------------------------
285
286besselTest = do
287 let (r,e) = bessel_J0_e 5.0
288 let expected = -0.17759677131433830434739701
289 assertBool "bessel_J0_e" ( abs (r-expected) < e )
290
291exponentialTest = do
292 let (v,e,err) = exp_e10_e 30.0
293 let expected = exp 30.0
294 assertBool "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )
295
296gammaTest = do
297 assertBool "gamma" (gamma 5 == 24.0)
298
299tests = TestList
300 [ TestCase $ besselTest
301 , TestCase $ exponentialTest
302 , TestCase $ gammaTest
303 , TestCase $ polySolveTest
304 , TestCase $ integrateTest
305 ]
306
307----------------------------------------------------------------------
308
309main = do
310 putStrLn "--------- general -----"
311 quickCheck (\(Sym m) -> m == (trans m:: Matrix BaseType))
312 quickCheck $ \l -> null l || (toList . fromList) l == (l :: [BaseType])
313
314 quickCheck $ \m -> m == asC (m :: Matrix BaseType)
315 quickCheck $ \m -> m == asFortran (m :: Matrix BaseType)
316 quickCheck $ \m -> m == (asC . asFortran) (m :: Matrix BaseType)
317 putStrLn "--------- MULTIPLY ----"
318 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: Matrix BaseType)
319 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: Matrix BaseType))
320 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: Matrix BaseType)
321 putStrLn "--------- SVD ---------"
322 quickCheck (svdTestR svdR)
323 quickCheck (svdTestR svdRdd)
324-- quickCheck (svdTestR svdg)
325 quickCheck svdTestC
326 putStrLn "--------- EIG ---------"
327 quickCheck eigTestC
328 quickCheck eigTestR
329 quickCheck eigTestS
330 quickCheck eigTestH
331 putStrLn "--------- SOLVE ---------"
332 quickCheck (linearSolveSQTest linearSolveR (singular svdR'))
333 quickCheck (linearSolveSQTest linearSolveC (singular svdC'))
334 quickCheck (pinvTest pinvR)
335 quickCheck (pinvTest pinvC)
336 quickCheck (pinvTest pinvSVDR)
337 quickCheck (pinvTest pinvSVDC)
338 putStrLn "--------- VEC OPER ------"
339 quickCheck arit1
340 quickCheck arit2
341 quickCheck arit3
342 quickCheck arit4
343 putStrLn "--------- GSL ------"
344 runTestTT tests
345 quickCheck $ \v -> ifft (fft v) |~| v
346
347kk = (2><2)
348 [ 1.0, 0.0
349 , -1.5, 1.0 ::Double]
350
351v = 11 |> [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0::Double]
352
353pol = [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5]
354
355mm = (2><2)
356 [ 0.5, 0.0
357 , 0.0, 0.0 ] :: Matrix Double
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-29 09:53:35 +0000