diff options
author | Alberto Ruiz <aruiz@um.es> | 2007-06-25 07:32:56 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2007-06-25 07:32:56 +0000 |
commit | 1871acb835b4fc164bcff3f6e7467884b87fbd0f (patch) | |
tree | ac1028d40778bbae532c3915276b5af21ba5f5cb /examples/tests.hs | |
parent | 3d5d6f06598aac00906c93ac5358e68697c47fc7 (diff) |
l.a. algorithms, etc.
Diffstat (limited to 'examples/tests.hs')
-rw-r--r-- | examples/tests.hs | 152 |
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 | |||
16 | import GSL.Polynomials | 16 | import GSL.Polynomials |
17 | import LAPACK | 17 | import LAPACK |
18 | import Test.QuickCheck | 18 | import Test.QuickCheck |
19 | import Test.HUnit | 19 | import Test.HUnit hiding ((~:)) |
20 | import Complex | 20 | import Complex |
21 | import LinearAlgebra.Algorithms | ||
22 | import GSL.Matrix | ||
23 | import Data.Packed.Instances hiding ((<>)) | ||
24 | |||
25 | dist :: (Normed t, Num t) => t -> t -> Double | ||
26 | dist a b = norm (a-b) | ||
27 | |||
28 | infixl 4 |~| | ||
29 | a |~| b = a :~8~: b | ||
30 | |||
31 | data Aprox a = (:~) a Int | ||
32 | |||
33 | (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool | ||
34 | a :~n~: b = dist a b < 10^^(-n) | ||
21 | 35 | ||
22 | 36 | ||
23 | {- | 37 | {- |
24 | -- Bravo por quickCheck! | 38 | -- Bravo por quickCheck! |
25 | 39 | ||
26 | pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m ~~ ident (cols m) | 40 | pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m ~~ ident (cols m) |
27 | where infix 2 ~~ | 41 | where infix 2 ~~ |
@@ -32,7 +46,7 @@ pinvProp2 tol m = 0 < r && r <= c ==> (r==c) `trivial` (m <> pinv m <> m ~~ m) | |||
32 | c = cols m | 46 | c = cols m |
33 | infix 2 ~~ | 47 | infix 2 ~~ |
34 | (~~) = approxEqual tol | 48 | (~~) = approxEqual tol |
35 | 49 | ||
36 | nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros | 50 | nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros |
37 | where nt = trans (nullspace m) | 51 | where nt = trans (nullspace m) |
38 | cr = corank m | 52 | cr = corank m |
@@ -49,31 +63,28 @@ mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)] | |||
49 | af = (2>|<3) [1,4,2,5,3,6::Double] | 63 | af = (2>|<3) [1,4,2,5,3,6::Double] |
50 | bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double] | 64 | bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double] |
51 | 65 | ||
52 | a |=| b = rows a == rows b && | ||
53 | cols a == cols b && | ||
54 | toList (cdat a) == toList (cdat b) && | ||
55 | toList (fdat a) == toList (fdat b) | ||
56 | 66 | ||
67 | {- | ||
57 | aprox fun a b = rows a == rows b && | 68 | aprox fun a b = rows a == rows b && |
58 | cols a == cols b && | 69 | cols a == cols b && |
59 | eps > aproxL fun (toList (t a)) (toList (t b)) | 70 | epsTol > aproxL fun (toList (t a)) (toList (t b)) |
60 | where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat | 71 | where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat |
61 | 72 | ||
62 | aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length v1) | 73 | aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length v1) |
63 | 74 | ||
64 | normVR a b = toScalarR AbsSum (vectorZipR Sub a b) | 75 | normVR a b = toScalarR AbsSum (vectorZipR Sub a b) |
65 | 76 | ||
66 | a |~| b = rows a == rows b && cols a == cols b && eps > normVR (t a) (t b) | 77 | a |~| b = rows a == rows b && cols a == cols b && epsTol > normVR (t a) (t b) |
67 | where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat | 78 | where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat |
68 | 79 | ||
69 | (|~~|) = aprox magnitude | 80 | (|~~|) = aprox magnitude |
70 | 81 | ||
71 | v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2 | 82 | v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2 |
72 | 83 | ||
73 | u ~|~ v = normVR u v < eps | 84 | u ~|~ v = normVR u v < epsTol |
74 | 85 | -} | |
75 | 86 | ||
76 | eps = 1E-8::Double | 87 | epsTol = 1E-8::Double |
77 | 88 | ||
78 | asFortran m = (rows m >|< cols m) $ toList (fdat m) | 89 | asFortran m = (rows m >|< cols m) $ toList (fdat m) |
79 | asC m = (rows m >< cols m) $ toList (cdat m) | 90 | asC m = (rows m >< cols m) $ toList (cdat m) |
@@ -81,6 +92,9 @@ asC m = (rows m >< cols m) $ toList (cdat m) | |||
81 | mulC a b = multiply RowMajor a b | 92 | mulC a b = multiply RowMajor a b |
82 | mulF a b = multiply ColumnMajor a b | 93 | mulF a b = multiply ColumnMajor a b |
83 | 94 | ||
95 | infixl 7 <> | ||
96 | a <> b = mulF a b | ||
97 | |||
84 | cc = mulC ac bf | 98 | cc = mulC ac bf |
85 | cf = mulF af bc | 99 | cf = mulF af bc |
86 | 100 | ||
@@ -133,14 +147,14 @@ data Sym a = Sym (Matrix a) deriving Show | |||
133 | instance (Field a, Arbitrary a, Num a) => Arbitrary (Sym a) where | 147 | instance (Field a, Arbitrary a, Num a) => Arbitrary (Sym a) where |
134 | arbitrary = do | 148 | arbitrary = do |
135 | SqM m <- arbitrary | 149 | SqM m <- arbitrary |
136 | return $ Sym (m `addM` trans m) | 150 | return $ Sym (m + trans m) |
137 | coarbitrary = undefined | 151 | coarbitrary = undefined |
138 | 152 | ||
139 | data Her = Her (Matrix (Complex Double)) deriving Show | 153 | data Her = Her (Matrix (Complex Double)) deriving Show |
140 | instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where | 154 | instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where |
141 | arbitrary = do | 155 | arbitrary = do |
142 | SqM m <- arbitrary | 156 | SqM m <- arbitrary |
143 | return $ Her (m `addM` (liftMatrix conj) (trans m)) | 157 | return $ Her (m + conjTrans m) |
144 | coarbitrary = undefined | 158 | coarbitrary = undefined |
145 | 159 | ||
146 | data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show | 160 | data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show |
@@ -171,49 +185,40 @@ instance (Field a, Arbitrary a) => Arbitrary (PairV a) where | |||
171 | coarbitrary = undefined | 185 | coarbitrary = undefined |
172 | 186 | ||
173 | 187 | ||
174 | addM m1 m2 = liftMatrix2 add m1 m2 | ||
175 | |||
176 | 188 | ||
177 | type BaseType = Double | 189 | type BaseType = Double |
178 | 190 | ||
179 | svdTestR fun prod m = u <> s <> trans v |~| m | 191 | svdTestR fun m = u <> s <> trans v |~| m |
180 | && u <> trans u |~| ident (rows m) | 192 | && u <> trans u |~| ident (rows m) |
181 | && v <> trans v |~| ident (cols m) | 193 | && v <> trans v |~| ident (cols m) |
182 | where (u,s,v) = fun m | 194 | where (u,s,v) = fun m |
183 | (<>) = prod | ||
184 | 195 | ||
185 | 196 | ||
186 | svdTestC prod m = u <> s' <> (trans v) |~~| m | 197 | svdTestC m = u <> s' <> (trans v) |~| m |
187 | && u <> (liftMatrix conj) (trans u) |~~| ident (rows m) | 198 | && u <> conjTrans u |~| ident (rows m) |
188 | && v <> (liftMatrix conj) (trans v) |~~| ident (cols m) | 199 | && v <> conjTrans v |~| ident (cols m) |
189 | where (u,s,v) = svdC m | 200 | where (u,s,v) = svdC m |
190 | (<>) = prod | ||
191 | s' = liftMatrix comp s | 201 | s' = liftMatrix comp s |
192 | 202 | ||
193 | eigTestC prod (SqM m) = (m <> v) |~~| (v <> diag s) | 203 | --svdg' m = (u,s',v) where |
194 | && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant 1 (rows m) --normalized | 204 | |
205 | eigTestC (SqM m) = (m <> v) |~| (v <> diag s) | ||
206 | && takeDiag (conjTrans v <> v) |~| constant 1 (rows m) --normalized | ||
195 | where (s,v) = eigC m | 207 | where (s,v) = eigC m |
196 | (<>) = prod | ||
197 | 208 | ||
198 | eigTestR prod (SqM m) = (liftMatrix comp m <> v) |~~| (v <> diag s) | 209 | eigTestR (SqM m) = (liftMatrix comp m <> v) |~| (v <> diag s) |
199 | -- && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant 1 (rows m) --normalized ??? | 210 | -- && takeDiag ((liftMatrix conj (trans v)) <> v) |~| constant 1 (rows m) --normalized ??? |
200 | where (s,v) = eigR m | 211 | where (s,v) = eigR m |
201 | (<>) = prod | ||
202 | 212 | ||
203 | eigTestS prod (Sym m) = (m <> v) |~| (v <> diag s) | 213 | eigTestS (Sym m) = (m <> v) |~| (v <> diag s) |
204 | && v <> trans v |~| ident (cols m) | 214 | && v <> trans v |~| ident (cols m) |
205 | where (s,v) = eigS m | 215 | where (s,v) = eigS m |
206 | (<>) = prod | ||
207 | 216 | ||
208 | eigTestH prod (Her m) = (m <> v) |~~| (v <> diag (comp s)) | 217 | eigTestH (Her m) = (m <> v) |~| (v <> diag (comp s)) |
209 | && v <> (liftMatrix conj) (trans v) |~~| ident (cols m) | 218 | && v <> conjTrans v |~| ident (cols m) |
210 | where (s,v) = eigH m | 219 | where (s,v) = eigH m |
211 | (<>) = prod | ||
212 | |||
213 | linearSolveSQTest fun eqfun singu prod (PairSM a b) = singu a || (a <> fun a b) ==== b | ||
214 | where (<>) = prod | ||
215 | (====) = eqfun | ||
216 | 220 | ||
221 | linearSolveSQTest fun singu (PairSM a b) = singu a || (a <> fun a b) |~| b | ||
217 | 222 | ||
218 | prec = 1E-15 | 223 | prec = 1E-15 |
219 | 224 | ||
@@ -237,8 +242,7 @@ identC n = toComplex(ident n, (0::Double) <>ident n) | |||
237 | 242 | ||
238 | -------------------------------------------------------------------- | 243 | -------------------------------------------------------------------- |
239 | 244 | ||
240 | pinvTest f feq m = (m <> f m <> m) `feq` m | 245 | pinvTest f m = (m <> f m <> m) |~| m |
241 | where (<>) = mulF | ||
242 | 246 | ||
243 | pinvR m = linearSolveLSR m (ident (rows m)) | 247 | pinvR m = linearSolveLSR m (ident (rows m)) |
244 | pinvC m = linearSolveLSC m (ident (rows m)) | 248 | pinvC m = linearSolveLSC m (ident (rows m)) |
@@ -252,7 +256,7 @@ pinvSVDC m = linearSolveSVDC Nothing m (ident (rows m)) | |||
252 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs | 256 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs |
253 | 257 | ||
254 | polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) | 258 | polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) |
255 | where l1 |~~| l2 = eps > aproxL magnitude l1 l2 | 259 | |
256 | 260 | ||
257 | polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4]) | 261 | polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4]) |
258 | 262 | ||
@@ -267,17 +271,17 @@ quad2 f a b g1 g2 = quad h a b | |||
267 | volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) | 271 | volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) |
268 | 0 r (const 0) (\x->sqrt (r*r-x*x)) | 272 | 0 r (const 0) (\x->sqrt (r*r-x*x)) |
269 | 273 | ||
270 | integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < eps) | 274 | integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol) |
271 | 275 | ||
272 | 276 | ||
273 | --------------------------------------------------------------------- | 277 | --------------------------------------------------------------------- |
274 | 278 | ||
275 | arit1 u = vectorMapValR PowVS 2 (vectorMapR Sin u) | 279 | arit1 u = vectorMapValR PowVS 2 (vectorMapR Sin u) |
276 | `add` vectorMapValR PowVS 2 (vectorMapR Cos u) | 280 | `add` vectorMapValR PowVS 2 (vectorMapR Cos u) |
277 | ~|~ constant 1 (dim u) | 281 | |~| constant 1 (dim u) |
278 | 282 | ||
279 | arit2 u = (vectorMapR Cos u) `mul` (vectorMapR Tan u) | 283 | arit2 u = (vectorMapR Cos u) `mul` (vectorMapR Tan u) |
280 | ~|~ vectorMapR Sin u | 284 | |~| vectorMapR Sin u |
281 | 285 | ||
282 | 286 | ||
283 | -- arit3 (PairV u v) = | 287 | -- arit3 (PairV u v) = |
@@ -305,50 +309,48 @@ tests = TestList | |||
305 | 309 | ||
306 | main = do | 310 | main = do |
307 | putStrLn "--------- general -----" | 311 | putStrLn "--------- general -----" |
308 | quickCheck (\(Sym m) -> m |=| (trans m:: Matrix BaseType)) | 312 | quickCheck (\(Sym m) -> m == (trans m:: Matrix BaseType)) |
309 | quickCheck $ \l -> null l || (toList . fromList) l == (l :: [BaseType]) | 313 | quickCheck $ \l -> null l || (toList . fromList) l == (l :: [BaseType]) |
310 | 314 | ||
311 | quickCheck $ \m -> m |=| asC (m :: Matrix BaseType) | 315 | quickCheck $ \m -> m == asC (m :: Matrix BaseType) |
312 | quickCheck $ \m -> m |=| asFortran (m :: Matrix BaseType) | 316 | quickCheck $ \m -> m == asFortran (m :: Matrix BaseType) |
313 | quickCheck $ \m -> m |=| (asC . asFortran) (m :: Matrix BaseType) | 317 | quickCheck $ \m -> m == (asC . asFortran) (m :: Matrix BaseType) |
314 | putStrLn "--------- MULTIPLY ----" | 318 | putStrLn "--------- MULTIPLY ----" |
315 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| mulF m1 (m2 :: Matrix BaseType) | 319 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: Matrix BaseType) |
316 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| trans (mulF (trans m2) (trans m1 :: Matrix BaseType)) | 320 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: Matrix BaseType)) |
317 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| multiplyG m1 (m2 :: Matrix BaseType) | 321 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: Matrix BaseType) |
318 | putStrLn "--------- SVD ---------" | 322 | putStrLn "--------- SVD ---------" |
319 | quickCheck (svdTestR svdR mulC) | 323 | quickCheck (svdTestR svdR) |
320 | quickCheck (svdTestR svdR mulF) | 324 | quickCheck (svdTestR svdRdd) |
321 | quickCheck (svdTestR svdRdd mulC) | 325 | -- quickCheck (svdTestR svdg) |
322 | quickCheck (svdTestR svdRdd mulF) | 326 | quickCheck svdTestC |
323 | quickCheck (svdTestC mulC) | ||
324 | quickCheck (svdTestC mulF) | ||
325 | putStrLn "--------- EIG ---------" | 327 | putStrLn "--------- EIG ---------" |
326 | quickCheck (eigTestC mulC) | 328 | quickCheck eigTestC |
327 | quickCheck (eigTestC mulF) | 329 | quickCheck eigTestR |
328 | quickCheck (eigTestR mulC) | 330 | quickCheck eigTestS |
329 | quickCheck (eigTestR mulF) | 331 | quickCheck eigTestH |
330 | quickCheck (eigTestS mulC) | ||
331 | quickCheck (eigTestS mulF) | ||
332 | quickCheck (eigTestH mulC) | ||
333 | quickCheck (eigTestH mulF) | ||
334 | putStrLn "--------- SOLVE ---------" | 332 | putStrLn "--------- SOLVE ---------" |
335 | quickCheck (linearSolveSQTest linearSolveR (|~|) (singular svdR') mulC) | 333 | quickCheck (linearSolveSQTest linearSolveR (singular svdR')) |
336 | quickCheck (linearSolveSQTest linearSolveC (|~~|) (singular svdC') mulF) | 334 | quickCheck (linearSolveSQTest linearSolveC (singular svdC')) |
337 | quickCheck (pinvTest pinvR (|~|)) | 335 | quickCheck (pinvTest pinvR) |
338 | quickCheck (pinvTest pinvC (|~~|)) | 336 | quickCheck (pinvTest pinvC) |
339 | quickCheck (pinvTest pinvSVDR (|~|)) | 337 | quickCheck (pinvTest pinvSVDR) |
340 | quickCheck (pinvTest pinvSVDC (|~~|)) | 338 | quickCheck (pinvTest pinvSVDC) |
341 | putStrLn "--------- VEC OPER ------" | 339 | putStrLn "--------- VEC OPER ------" |
342 | quickCheck arit1 | 340 | quickCheck arit1 |
343 | quickCheck arit2 | 341 | quickCheck arit2 |
344 | putStrLn "--------- GSL ------" | 342 | putStrLn "--------- GSL ------" |
345 | runTestTT tests | 343 | runTestTT tests |
346 | quickCheck $ \v -> ifft (fft v) ~~ v | 344 | quickCheck $ \v -> ifft (fft v) |~| v |
347 | 345 | ||
348 | kk = (2><2) | 346 | kk = (2><2) |
349 | [ 1.0, 0.0 | 347 | [ 1.0, 0.0 |
350 | , -1.5, 1.0 ::Double] | 348 | , -1.5, 1.0 ::Double] |
351 | 349 | ||
352 | 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] | 350 | 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] |
351 | |||
352 | pol = [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5] | ||
353 | 353 | ||
354 | 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 | 354 | mm = (2><2) |
355 | [ 0.5, 0.0 | ||
356 | , 0.0, 0.0 ] :: Matrix Double | ||