diff options
author | Alberto Ruiz <aruiz@um.es> | 2007-09-21 18:10:59 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2007-09-21 18:10:59 +0000 |
commit | edfaf9e0d1dcfccc9015476510a23e8cf64288be (patch) | |
tree | 2b11f0a933a7cce2362aed26ac160312e6b9431a /examples/tests.hs | |
parent | 6cafd2f26a89008cc0db02e70e39f92a50ec4b4d (diff) |
algorithms refactoring
Diffstat (limited to 'examples/tests.hs')
-rw-r--r-- | examples/tests.hs | 357 |
1 files changed, 0 insertions, 357 deletions
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 | |||
9 | import Data.Packed.Internal((>|<), fdat, cdat, multiply', multiplyG, MatrixOrder(..)) | ||
10 | import GSL hiding (sin,cos,exp,choose) | ||
11 | import LinearAlgebra hiding ((<>)) | ||
12 | import Test.QuickCheck | ||
13 | import Test.HUnit hiding ((~:)) | ||
14 | |||
15 | |||
16 | dist :: (Normed t, Num t) => t -> t -> Double | ||
17 | dist a b = pnorm Infinity (a-b) | ||
18 | |||
19 | infixl 4 |~| | ||
20 | a |~| b = a :~8~: b | ||
21 | |||
22 | data Aprox a = (:~) a Int | ||
23 | |||
24 | (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool | ||
25 | a :~n~: b = dist a b < 10^^(-n) | ||
26 | |||
27 | |||
28 | {- | ||
29 | -- Bravo por quickCheck! | ||
30 | |||
31 | pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m ~~ ident (cols m) | ||
32 | where infix 2 ~~ | ||
33 | (~~) = approxEqual tol | ||
34 | |||
35 | pinvProp2 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 | |||
41 | nullspaceProp 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 | |||
49 | ac = (2><3) [1 .. 6::Double] | ||
50 | bc = (3><4) [7 .. 18::Double] | ||
51 | |||
52 | mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)] | ||
53 | |||
54 | af = (2>|<3) [1,4,2,5,3,6::Double] | ||
55 | bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double] | ||
56 | |||
57 | |||
58 | {- | ||
59 | aprox 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 | |||
64 | aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length v1) | ||
65 | |||
66 | normVR a b = toScalarR AbsSum (vectorZipR Sub a b) | ||
67 | |||
68 | a |~| 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 | |||
73 | v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2 | ||
74 | |||
75 | u ~|~ v = normVR u v < epsTol | ||
76 | -} | ||
77 | |||
78 | epsTol = 1E-8::Double | ||
79 | |||
80 | asFortran m = (rows m >|< cols m) $ toList (fdat m) | ||
81 | asC m = (rows m >< cols m) $ toList (cdat m) | ||
82 | |||
83 | mulC a b = multiply' RowMajor a b | ||
84 | mulF a b = multiply' ColumnMajor a b | ||
85 | |||
86 | identC = comp . ident | ||
87 | |||
88 | infixl 7 <> | ||
89 | a <> b = mulF a b | ||
90 | |||
91 | cc = mulC ac bf | ||
92 | cf = mulF af bc | ||
93 | |||
94 | r = mulC cc (trans cf) | ||
95 | |||
96 | rd = (2><2) | ||
97 | [ 27736.0, 65356.0 | ||
98 | , 65356.0, 154006.0 ::Double] | ||
99 | |||
100 | instance (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 | |||
107 | instance (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 | |||
119 | data PairM a = PairM (Matrix a) (Matrix a) deriving Show | ||
120 | instance (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 | |||
131 | data SqM a = SqM (Matrix a) deriving Show | ||
132 | instance (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 | |||
139 | data Sym a = Sym (Matrix a) deriving Show | ||
140 | instance (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 | |||
146 | data Her = Her (Matrix (Complex Double)) deriving Show | ||
147 | instance {-(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 | |||
153 | data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show | ||
154 | instance (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 | |||
164 | instance (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 | |||
171 | data PairV a = PairV (Vector a) (Vector a) | ||
172 | instance (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 | |||
182 | type BaseType = Complex Double | ||
183 | |||
184 | svdTestR 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 | |||
190 | svdTestC 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 | |||
198 | eigTestC (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 | |||
202 | eigTestR (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 | |||
206 | eigTestS (Sym m) = (m <> v) |~| (v <> diag s) | ||
207 | && v <> trans v |~| ident (cols m) | ||
208 | where (s,v) = eigS m | ||
209 | |||
210 | eigTestH (Her m) = (m <> v) |~| (v <> diag (comp s)) | ||
211 | && v <> conjTrans v |~| identC (cols m) | ||
212 | where (s,v) = eigH m | ||
213 | |||
214 | linearSolveSQTest fun singu (PairSM a b) = singu a || (a <> fun a b) |~| b | ||
215 | |||
216 | prec = 1E-15 | ||
217 | |||
218 | singular 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 | {- | ||
225 | invTest msg m = do | ||
226 | assertBool msg $ m <> inv m =~= ident (rows m) | ||
227 | |||
228 | invComplexTest msg m = do | ||
229 | assertBool msg $ m <> invC m =~= identC (rows m) | ||
230 | |||
231 | invC m = linearSolveC m (identC (rows m)) | ||
232 | |||
233 | identC n = toComplex(ident n, (0::Double) <>ident n) | ||
234 | -} | ||
235 | |||
236 | -------------------------------------------------------------------- | ||
237 | |||
238 | pinvTest f m = (m <> f m <> m) |~| m | ||
239 | |||
240 | pinvR m = linearSolveLSR m (ident (rows m)) | ||
241 | pinvC m = linearSolveLSC m (identC (rows m)) | ||
242 | |||
243 | pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m)) | ||
244 | |||
245 | pinvSVDC m = linearSolveSVDC Nothing m (identC (rows m)) | ||
246 | |||
247 | -------------------------------------------------------------------- | ||
248 | |||
249 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs | ||
250 | |||
251 | polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) | ||
252 | |||
253 | |||
254 | polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4]) | ||
255 | |||
256 | --------------------------------------------------------------------- | ||
257 | |||
258 | quad f a b = fst $ integrateQAGS 1E-9 100 f a b | ||
259 | |||
260 | -- A multiple integral can be easily defined using partial application | ||
261 | quad2 f a b g1 g2 = quad h a b | ||
262 | where h x = quad (f x) (g1 x) (g2 x) | ||
263 | |||
264 | volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) | ||
265 | 0 r (const 0) (\x->sqrt (r*r-x*x)) | ||
266 | |||
267 | integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol) | ||
268 | |||
269 | |||
270 | --------------------------------------------------------------------- | ||
271 | |||
272 | arit1 u = sin u ^ 2 + cos u ^ 2 |~| 1 | ||
273 | where _ = u :: Vector Double | ||
274 | |||
275 | arit2 u = sin u ** 2 + cos u ** 2 |~| 1 | ||
276 | where _ = u :: Vector Double | ||
277 | |||
278 | arit3 u = cos u * tan u |~| sin u | ||
279 | where _ = u :: Vector Double | ||
280 | |||
281 | arit4 u = (cos u * tan u) :~6~: sin u | ||
282 | where _ = u :: Vector (Complex Double) | ||
283 | |||
284 | --------------------------------------------------------------------- | ||
285 | |||
286 | besselTest = 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 | |||
291 | exponentialTest = 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 | |||
296 | gammaTest = do | ||
297 | assertBool "gamma" (gamma 5 == 24.0) | ||
298 | |||
299 | tests = TestList | ||
300 | [ TestCase $ besselTest | ||
301 | , TestCase $ exponentialTest | ||
302 | , TestCase $ gammaTest | ||
303 | , TestCase $ polySolveTest | ||
304 | , TestCase $ integrateTest | ||
305 | ] | ||
306 | |||
307 | ---------------------------------------------------------------------- | ||
308 | |||
309 | main = 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 | |||
347 | kk = (2><2) | ||
348 | [ 1.0, 0.0 | ||
349 | , -1.5, 1.0 ::Double] | ||
350 | |||
351 | 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] | ||
352 | |||
353 | pol = [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5] | ||
354 | |||
355 | mm = (2><2) | ||
356 | [ 0.5, 0.0 | ||
357 | , 0.0, 0.0 ] :: Matrix Double | ||