diff options
author | Alberto Ruiz <aruiz@um.es> | 2011-12-14 13:08:43 +0100 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2011-12-14 13:08:43 +0100 |
commit | 77552d080e88fc70312f55fd3303fac3464ab46e (patch) | |
tree | 1dc87dd22ce0da0f1807765568fbc04285bf3621 /lib/Numeric/LinearAlgebra | |
parent | c3bda2d38c432fb53ce456cba295b097fd4d6ad1 (diff) |
new package hmatrix-tests
Diffstat (limited to 'lib/Numeric/LinearAlgebra')
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests.hs | 723 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/Instances.hs | 249 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/Properties.hs | 272 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h | 33 |
4 files changed, 0 insertions, 1277 deletions
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs deleted file mode 100644 index e859450..0000000 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ /dev/null | |||
@@ -1,723 +0,0 @@ | |||
1 | {-# LANGUAGE CPP #-} | ||
2 | {-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-} | ||
3 | ----------------------------------------------------------------------------- | ||
4 | {- | | ||
5 | Module : Numeric.LinearAlgebra.Tests | ||
6 | Copyright : (c) Alberto Ruiz 2007-9 | ||
7 | License : GPL-style | ||
8 | |||
9 | Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
10 | Stability : provisional | ||
11 | Portability : portable | ||
12 | |||
13 | Some tests. | ||
14 | |||
15 | -} | ||
16 | |||
17 | module Numeric.LinearAlgebra.Tests( | ||
18 | -- module Numeric.LinearAlgebra.Tests.Instances, | ||
19 | -- module Numeric.LinearAlgebra.Tests.Properties, | ||
20 | qCheck, runTests, runBenchmarks, findNaN | ||
21 | --, runBigTests | ||
22 | ) where | ||
23 | |||
24 | import Data.Packed.Random | ||
25 | import Numeric.LinearAlgebra | ||
26 | import Numeric.LinearAlgebra.LAPACK | ||
27 | import Numeric.LinearAlgebra.Tests.Instances | ||
28 | import Numeric.LinearAlgebra.Tests.Properties | ||
29 | import Test.HUnit hiding ((~:),test,Testable,State) | ||
30 | import System.Info | ||
31 | import Data.List(foldl1') | ||
32 | import Numeric.GSL | ||
33 | import Prelude hiding ((^)) | ||
34 | import qualified Prelude | ||
35 | import System.CPUTime | ||
36 | import Text.Printf | ||
37 | import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr) | ||
38 | import Control.Arrow((***)) | ||
39 | import Debug.Trace | ||
40 | |||
41 | #include "Tests/quickCheckCompat.h" | ||
42 | |||
43 | debug x = trace (show x) x | ||
44 | |||
45 | a ^ b = a Prelude.^ (b :: Int) | ||
46 | |||
47 | utest str b = TestCase $ assertBool str b | ||
48 | |||
49 | a ~~ b = fromList a |~| fromList b | ||
50 | |||
51 | feye n = flipud (ident n) :: Matrix Double | ||
52 | |||
53 | ----------------------------------------------------------- | ||
54 | |||
55 | detTest1 = det m == 26 | ||
56 | && det mc == 38 :+ (-3) | ||
57 | && det (feye 2) == -1 | ||
58 | where | ||
59 | m = (3><3) | ||
60 | [ 1, 2, 3 | ||
61 | , 4, 5, 7 | ||
62 | , 2, 8, 4 :: Double | ||
63 | ] | ||
64 | mc = (3><3) | ||
65 | [ 1, 2, 3 | ||
66 | , 4, 5, 7 | ||
67 | , 2, 8, i | ||
68 | ] | ||
69 | |||
70 | detTest2 = inv1 |~| inv2 && [det1] ~~ [det2] | ||
71 | where | ||
72 | m = complex (feye 6) | ||
73 | inv1 = inv m | ||
74 | det1 = det m | ||
75 | (inv2,(lda,sa)) = invlndet m | ||
76 | det2 = sa * exp lda | ||
77 | |||
78 | -------------------------------------------------------------------- | ||
79 | |||
80 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs | ||
81 | |||
82 | polySolveProp p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) | ||
83 | |||
84 | --------------------------------------------------------------------- | ||
85 | |||
86 | quad f a b = fst $ integrateQAGS 1E-9 100 f a b | ||
87 | |||
88 | -- A multiple integral can be easily defined using partial application | ||
89 | quad2 f a b g1 g2 = quad h a b | ||
90 | where h x = quad (f x) (g1 x) (g2 x) | ||
91 | |||
92 | volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) | ||
93 | 0 r (const 0) (\x->sqrt (r*r-x*x)) | ||
94 | |||
95 | --------------------------------------------------------------------- | ||
96 | |||
97 | derivTest = abs (d (\x-> x * d (\y-> x+y) 1) 1 - 1) < 1E-10 | ||
98 | where d f x = fst $ derivCentral 0.01 f x | ||
99 | |||
100 | --------------------------------------------------------------------- | ||
101 | |||
102 | -- besselTest = utest "bessel_J0_e" ( abs (r-expected) < e ) | ||
103 | -- where (r,e) = bessel_J0_e 5.0 | ||
104 | -- expected = -0.17759677131433830434739701 | ||
105 | |||
106 | -- exponentialTest = utest "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 ) | ||
107 | -- where (v,e,_err) = exp_e10_e 30.0 | ||
108 | -- expected = exp 30.0 | ||
109 | |||
110 | --------------------------------------------------------------------- | ||
111 | |||
112 | nd1 = (3><3) [ 1/2, 1/4, 1/4 | ||
113 | , 0/1, 1/2, 1/4 | ||
114 | , 1/2, 1/4, 1/2 :: Double] | ||
115 | |||
116 | nd2 = (2><2) [1, 0, 1, 1:: Complex Double] | ||
117 | |||
118 | expmTest1 = expm nd1 :~14~: (3><3) | ||
119 | [ 1.762110887278176 | ||
120 | , 0.478085470590435 | ||
121 | , 0.478085470590435 | ||
122 | , 0.104719410945666 | ||
123 | , 1.709751181805343 | ||
124 | , 0.425725765117601 | ||
125 | , 0.851451530235203 | ||
126 | , 0.530445176063267 | ||
127 | , 1.814470592751009 ] | ||
128 | |||
129 | expmTest2 = expm nd2 :~15~: (2><2) | ||
130 | [ 2.718281828459045 | ||
131 | , 0.000000000000000 | ||
132 | , 2.718281828459045 | ||
133 | , 2.718281828459045 ] | ||
134 | |||
135 | --------------------------------------------------------------------- | ||
136 | |||
137 | minimizationTest = TestList | ||
138 | [ utest "minimization conjugatefr" (minim1 f df [5,7] ~~ [1,2]) | ||
139 | , utest "minimization nmsimplex2" (minim2 f [5,7] `elem` [24,25]) | ||
140 | ] | ||
141 | where f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30 | ||
142 | df [x,y] = [20*(x-1), 40*(y-2)] | ||
143 | minim1 g dg ini = fst $ minimizeD ConjugateFR 1E-3 30 1E-2 1E-4 g dg ini | ||
144 | minim2 g ini = rows $ snd $ minimize NMSimplex2 1E-2 30 [1,1] g ini | ||
145 | |||
146 | --------------------------------------------------------------------- | ||
147 | |||
148 | rootFindingTest = TestList [ utest "root Hybrids" (fst sol1 ~~ [1,1]) | ||
149 | , utest "root Newton" (rows (snd sol2) == 2) | ||
150 | ] | ||
151 | where sol1 = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5] | ||
152 | sol2 = rootJ Newton 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5] | ||
153 | rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ] | ||
154 | jacobian a b [x,_y] = [ [-a , 0] | ||
155 | , [-2*b*x, b] ] | ||
156 | |||
157 | --------------------------------------------------------------------- | ||
158 | |||
159 | odeTest = utest "ode" (last (toLists sol) ~~ [-1.7588880332411019, 8.364348908711941e-2]) | ||
160 | where sol = odeSolveV RK8pd 1E-6 1E-6 0 (l2v $ vanderpol 10) Nothing (fromList [1,0]) ts | ||
161 | ts = linspace 101 (0,100) | ||
162 | l2v f = \t -> fromList . f t . toList | ||
163 | vanderpol mu _t [x,y] = [y, -x + mu * y * (1-x^2) ] | ||
164 | |||
165 | --------------------------------------------------------------------- | ||
166 | |||
167 | fittingTest = utest "levmar" (ok1 && ok2) | ||
168 | where | ||
169 | xs = map return [0 .. 39] | ||
170 | sigma = 0.1 | ||
171 | ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs) | ||
172 | + scalar sigma * (randomVector 0 Gaussian 40) | ||
173 | dats = zip xs (zip ys (repeat sigma)) | ||
174 | dat = zip xs ys | ||
175 | |||
176 | expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b] | ||
177 | expModelDer [a,lambda,_b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] | ||
178 | |||
179 | sols = fst $ fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dats [1,0,0] | ||
180 | sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] | ||
181 | |||
182 | ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d | ||
183 | ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5 | ||
184 | |||
185 | ----------------------------------------------------- | ||
186 | |||
187 | mbCholTest = utest "mbCholTest" (ok1 && ok2) where | ||
188 | m1 = (2><2) [2,5,5,8 :: Double] | ||
189 | m2 = (2><2) [3,5,5,9 :: Complex Double] | ||
190 | ok1 = mbCholSH m1 == Nothing | ||
191 | ok2 = mbCholSH m2 == Just (chol m2) | ||
192 | |||
193 | --------------------------------------------------------------------- | ||
194 | |||
195 | randomTestGaussian = c :~1~: snd (meanCov dat) where | ||
196 | a = (3><3) [1,2,3, | ||
197 | 2,4,0, | ||
198 | -2,2,1] | ||
199 | m = 3 |> [1,2,3] | ||
200 | c = a <> trans a | ||
201 | dat = gaussianSample 7 (10^6) m c | ||
202 | |||
203 | randomTestUniform = c :~1~: snd (meanCov dat) where | ||
204 | c = diag $ 3 |> map ((/12).(^2)) [1,2,3] | ||
205 | dat = uniformSample 7 (10^6) [(0,1),(1,3),(3,6)] | ||
206 | |||
207 | --------------------------------------------------------------------- | ||
208 | |||
209 | rot :: Double -> Matrix Double | ||
210 | rot a = (3><3) [ c,0,s | ||
211 | , 0,1,0 | ||
212 | ,-s,0,c ] | ||
213 | where c = cos a | ||
214 | s = sin a | ||
215 | |||
216 | rotTest = fun (10^5) :~11~: rot 5E4 | ||
217 | where fun n = foldl1' (<>) (map rot angles) | ||
218 | where angles = toList $ linspace n (0,1) | ||
219 | |||
220 | --------------------------------------------------------------------- | ||
221 | -- vector <= 0.6.0.2 bug discovered by Patrick Perry | ||
222 | -- http://trac.haskell.org/vector/ticket/31 | ||
223 | |||
224 | offsetTest = y == y' where | ||
225 | x = fromList [0..3 :: Double] | ||
226 | y = subVector 1 3 x | ||
227 | (f,o,n) = unsafeToForeignPtr y | ||
228 | y' = unsafeFromForeignPtr f o n | ||
229 | |||
230 | --------------------------------------------------------------------- | ||
231 | |||
232 | normsVTest = TestList [ | ||
233 | utest "normv2CD" $ norm2PropC v | ||
234 | , utest "normv2CF" $ norm2PropC (single v) | ||
235 | #ifndef NONORMVTEST | ||
236 | , utest "normv2D" $ norm2PropR x | ||
237 | , utest "normv2F" $ norm2PropR (single x) | ||
238 | #endif | ||
239 | , utest "normv1CD" $ norm1 v == 8 | ||
240 | , utest "normv1CF" $ norm1 (single v) == 8 | ||
241 | , utest "normv1D" $ norm1 x == 6 | ||
242 | , utest "normv1F" $ norm1 (single x) == 6 | ||
243 | |||
244 | , utest "normvInfCD" $ normInf v == 5 | ||
245 | , utest "normvInfCF" $ normInf (single v) == 5 | ||
246 | , utest "normvInfD" $ normInf x == 3 | ||
247 | , utest "normvInfF" $ normInf (single x) == 3 | ||
248 | |||
249 | ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double) | ||
250 | x = fromList [1,2,-3] :: Vector Double | ||
251 | #ifndef NONORMVTEST | ||
252 | norm2PropR a = norm2 a =~= sqrt (dot a a) | ||
253 | #endif | ||
254 | norm2PropC a = norm2 a =~= realPart (sqrt (dot a (conj a))) | ||
255 | a =~= b = fromList [a] |~| fromList [b] | ||
256 | |||
257 | normsMTest = TestList [ | ||
258 | utest "norm2mCD" $ pnorm PNorm2 v =~= 8.86164970498005 | ||
259 | , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005 | ||
260 | , utest "norm2mD" $ pnorm PNorm2 x =~= 5.96667765076216 | ||
261 | , utest "norm2mF" $ pnorm PNorm2 (single x) =~= 5.96667765076216 | ||
262 | |||
263 | , utest "norm1mCD" $ pnorm PNorm1 v == 9 | ||
264 | , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9 | ||
265 | , utest "norm1mD" $ pnorm PNorm1 x == 7 | ||
266 | , utest "norm1mF" $ pnorm PNorm1 (single x) == 7 | ||
267 | |||
268 | , utest "normmInfCD" $ pnorm Infinity v == 12 | ||
269 | , utest "normmInfCF" $ pnorm Infinity (single v) == 12 | ||
270 | , utest "normmInfD" $ pnorm Infinity x == 8 | ||
271 | , utest "normmInfF" $ pnorm Infinity (single x) == 8 | ||
272 | |||
273 | , utest "normmFroCD" $ pnorm Frobenius v =~= 8.88819441731559 | ||
274 | , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559 | ||
275 | , utest "normmFroD" $ pnorm Frobenius x =~= 6.24499799839840 | ||
276 | , utest "normmFroF" $ pnorm Frobenius (single x) =~~= 6.24499799839840 | ||
277 | |||
278 | ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double) | ||
279 | x = (2><2) [1,2,-3,5] :: Matrix Double | ||
280 | a =~= b = fromList [a] :~10~: fromList [b] | ||
281 | a =~~= b = fromList [a] :~5~: fromList [b] | ||
282 | |||
283 | --------------------------------------------------------------------- | ||
284 | |||
285 | sumprodTest = TestList [ | ||
286 | utest "sumCD" $ sumElements z == 6 | ||
287 | , utest "sumCF" $ sumElements (single z) == 6 | ||
288 | , utest "sumD" $ sumElements v == 6 | ||
289 | , utest "sumF" $ sumElements (single v) == 6 | ||
290 | |||
291 | , utest "prodCD" $ prodProp z | ||
292 | , utest "prodCF" $ prodProp (single z) | ||
293 | , utest "prodD" $ prodProp v | ||
294 | , utest "prodF" $ prodProp (single v) | ||
295 | ] where v = fromList [1,2,3] :: Vector Double | ||
296 | z = fromList [1,2-i,3+i] | ||
297 | prodProp x = prodElements x == product (toList x) | ||
298 | |||
299 | --------------------------------------------------------------------- | ||
300 | |||
301 | chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where | ||
302 | ms = [ diag (fromList [1,2,3 :: Double]) | ||
303 | , konst 3 (3,5) | ||
304 | , (5><10) [1 .. ] | ||
305 | , konst 5 (10,2) | ||
306 | ] | ||
307 | |||
308 | --------------------------------------------------------------------- | ||
309 | |||
310 | conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m) | ||
311 | |||
312 | --------------------------------------------------------------------- | ||
313 | |||
314 | newtype State s a = State { runState :: s -> (a,s) } | ||
315 | |||
316 | instance Monad (State s) where | ||
317 | return a = State $ \s -> (a,s) | ||
318 | m >>= f = State $ \s -> let (a,s') = runState m s | ||
319 | in runState (f a) s' | ||
320 | |||
321 | state_get :: State s s | ||
322 | state_get = State $ \s -> (s,s) | ||
323 | |||
324 | state_put :: s -> State s () | ||
325 | state_put s = State $ \_ -> ((),s) | ||
326 | |||
327 | evalState :: State s a -> s -> a | ||
328 | evalState m s = let (a,s') = runState m s | ||
329 | in seq s' a | ||
330 | |||
331 | newtype MaybeT m a = MaybeT { runMaybeT :: m (Maybe a) } | ||
332 | |||
333 | instance Monad m => Monad (MaybeT m) where | ||
334 | return a = MaybeT $ return $ Just a | ||
335 | m >>= f = MaybeT $ do | ||
336 | res <- runMaybeT m | ||
337 | case res of | ||
338 | Nothing -> return Nothing | ||
339 | Just r -> runMaybeT (f r) | ||
340 | fail _ = MaybeT $ return Nothing | ||
341 | |||
342 | lift_maybe m = MaybeT $ do | ||
343 | res <- m | ||
344 | return $ Just res | ||
345 | |||
346 | -- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs | ||
347 | --successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool | ||
348 | successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0) | ||
349 | where stp e = do | ||
350 | ep <- lift_maybe $ state_get | ||
351 | if t e ep | ||
352 | then lift_maybe $ state_put e | ||
353 | else (fail "successive_ test failed") | ||
354 | |||
355 | -- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input | ||
356 | --successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b | ||
357 | successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0) | ||
358 | where stp e = do | ||
359 | ep <- state_get | ||
360 | state_put e | ||
361 | return $ f ep e | ||
362 | |||
363 | |||
364 | succTest = utest "successive" $ | ||
365 | successive_ (>) (fromList [1 :: Double,2,3,4]) == True | ||
366 | && successive_ (>) (fromList [1 :: Double,3,2,4]) == False | ||
367 | && successive (+) (fromList [1..10 :: Double]) == 9 |> [3,5,7,9,11,13,15,17,19] | ||
368 | |||
369 | --------------------------------------------------------------------- | ||
370 | |||
371 | findAssocTest = utest "findAssoc" ok | ||
372 | where | ||
373 | ok = m1 == m2 | ||
374 | m1 = assoc (6,6) 7 $ zip (find (>0) (ident 5 :: Matrix Float)) [10 ..] :: Matrix Double | ||
375 | m2 = diagRect 7 (fromList[10..14]) 6 6 | ||
376 | |||
377 | --------------------------------------------------------------------- | ||
378 | |||
379 | condTest = utest "cond" ok | ||
380 | where | ||
381 | ok = step v * v == cond v 0 0 0 v | ||
382 | v = fromList [-7 .. 7 ] :: Vector Float | ||
383 | |||
384 | --------------------------------------------------------------------- | ||
385 | |||
386 | conformTest = utest "conform" ok | ||
387 | where | ||
388 | ok = 1 + row [1,2,3] + col [10,20,30,40] + (4><3) [1..] | ||
389 | == (4><3) [13,15,17 | ||
390 | ,26,28,30 | ||
391 | ,39,41,43 | ||
392 | ,52,54,56] | ||
393 | row = asRow . fromList | ||
394 | col = asColumn . fromList :: [Double] -> Matrix Double | ||
395 | |||
396 | --------------------------------------------------------------------- | ||
397 | |||
398 | accumTest = utest "accum" ok | ||
399 | where | ||
400 | x = ident 3 :: Matrix Double | ||
401 | ok = accum x (+) [((1,2),7), ((2,2),3)] | ||
402 | == (3><3) [1,0,0 | ||
403 | ,0,1,7 | ||
404 | ,0,0,4] | ||
405 | && | ||
406 | toList (flatten x) == [1,0,0,0,1,0,0,0,1] | ||
407 | |||
408 | --------------------------------------------------------------------- | ||
409 | |||
410 | -- | All tests must pass with a maximum dimension of about 20 | ||
411 | -- (some tests may fail with bigger sizes due to precision loss). | ||
412 | runTests :: Int -- ^ maximum dimension | ||
413 | -> IO () | ||
414 | runTests n = do | ||
415 | setErrorHandlerOff | ||
416 | let test p = qCheck n p | ||
417 | putStrLn "------ mult Double" | ||
418 | test (multProp1 10 . rConsist) | ||
419 | test (multProp1 10 . cConsist) | ||
420 | test (multProp2 10 . rConsist) | ||
421 | test (multProp2 10 . cConsist) | ||
422 | putStrLn "------ mult Float" | ||
423 | test (multProp1 6 . (single *** single) . rConsist) | ||
424 | test (multProp1 6 . (single *** single) . cConsist) | ||
425 | test (multProp2 6 . (single *** single) . rConsist) | ||
426 | test (multProp2 6 . (single *** single) . cConsist) | ||
427 | putStrLn "------ sub-trans" | ||
428 | test (subProp . rM) | ||
429 | test (subProp . cM) | ||
430 | putStrLn "------ ctrans" | ||
431 | test (conjuTest . cM) | ||
432 | test (conjuTest . zM) | ||
433 | putStrLn "------ lu" | ||
434 | test (luProp . rM) | ||
435 | test (luProp . cM) | ||
436 | putStrLn "------ inv (linearSolve)" | ||
437 | test (invProp . rSqWC) | ||
438 | test (invProp . cSqWC) | ||
439 | putStrLn "------ luSolve" | ||
440 | test (linearSolveProp (luSolve.luPacked) . rSqWC) | ||
441 | test (linearSolveProp (luSolve.luPacked) . cSqWC) | ||
442 | putStrLn "------ cholSolve" | ||
443 | test (linearSolveProp (cholSolve.chol) . rPosDef) | ||
444 | test (linearSolveProp (cholSolve.chol) . cPosDef) | ||
445 | putStrLn "------ luSolveLS" | ||
446 | test (linearSolveProp linearSolveLS . rSqWC) | ||
447 | test (linearSolveProp linearSolveLS . cSqWC) | ||
448 | test (linearSolveProp2 linearSolveLS . rConsist) | ||
449 | test (linearSolveProp2 linearSolveLS . cConsist) | ||
450 | putStrLn "------ pinv (linearSolveSVD)" | ||
451 | test (pinvProp . rM) | ||
452 | test (pinvProp . cM) | ||
453 | putStrLn "------ det" | ||
454 | test (detProp . rSqWC) | ||
455 | test (detProp . cSqWC) | ||
456 | putStrLn "------ svd" | ||
457 | test (svdProp1 . rM) | ||
458 | test (svdProp1 . cM) | ||
459 | test (svdProp1a svdR) | ||
460 | test (svdProp1a svdC) | ||
461 | test (svdProp1a svdRd) | ||
462 | test (svdProp1b svdR) | ||
463 | test (svdProp1b svdC) | ||
464 | test (svdProp1b svdRd) | ||
465 | test (svdProp2 thinSVDR) | ||
466 | test (svdProp2 thinSVDC) | ||
467 | test (svdProp2 thinSVDRd) | ||
468 | test (svdProp2 thinSVDCd) | ||
469 | test (svdProp3 . rM) | ||
470 | test (svdProp3 . cM) | ||
471 | test (svdProp4 . rM) | ||
472 | test (svdProp4 . cM) | ||
473 | test (svdProp5a) | ||
474 | test (svdProp5b) | ||
475 | test (svdProp6a) | ||
476 | test (svdProp6b) | ||
477 | test (svdProp7 . rM) | ||
478 | test (svdProp7 . cM) | ||
479 | putStrLn "------ svdCd" | ||
480 | #ifdef NOZGESDD | ||
481 | putStrLn "Omitted" | ||
482 | #else | ||
483 | test (svdProp1a svdCd) | ||
484 | test (svdProp1b svdCd) | ||
485 | #endif | ||
486 | putStrLn "------ eig" | ||
487 | test (eigSHProp . rHer) | ||
488 | test (eigSHProp . cHer) | ||
489 | test (eigProp . rSq) | ||
490 | test (eigProp . cSq) | ||
491 | test (eigSHProp2 . rHer) | ||
492 | test (eigSHProp2 . cHer) | ||
493 | test (eigProp2 . rSq) | ||
494 | test (eigProp2 . cSq) | ||
495 | putStrLn "------ nullSpace" | ||
496 | test (nullspaceProp . rM) | ||
497 | test (nullspaceProp . cM) | ||
498 | putStrLn "------ qr" | ||
499 | test (qrProp . rM) | ||
500 | test (qrProp . cM) | ||
501 | test (rqProp . rM) | ||
502 | test (rqProp . cM) | ||
503 | test (rqProp1 . cM) | ||
504 | test (rqProp2 . cM) | ||
505 | test (rqProp3 . cM) | ||
506 | putStrLn "------ hess" | ||
507 | test (hessProp . rSq) | ||
508 | test (hessProp . cSq) | ||
509 | putStrLn "------ schur" | ||
510 | test (schurProp2 . rSq) | ||
511 | test (schurProp1 . cSq) | ||
512 | putStrLn "------ chol" | ||
513 | test (cholProp . rPosDef) | ||
514 | test (cholProp . cPosDef) | ||
515 | test (exactProp . rPosDef) | ||
516 | test (exactProp . cPosDef) | ||
517 | putStrLn "------ expm" | ||
518 | test (expmDiagProp . complex. rSqWC) | ||
519 | test (expmDiagProp . cSqWC) | ||
520 | putStrLn "------ fft" | ||
521 | test (\v -> ifft (fft v) |~| v) | ||
522 | putStrLn "------ vector operations - Double" | ||
523 | test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM)) | ||
524 | test $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::CM)) . liftMatrix makeUnitary | ||
525 | test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) | ||
526 | test (\u -> cos u * tan u |~| sin (u::RM)) | ||
527 | test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary | ||
528 | putStrLn "------ vector operations - Float" | ||
529 | test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM)) | ||
530 | test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary | ||
531 | test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM)) | ||
532 | test (\u -> cos u * tan u |~~| sin (u::FM)) | ||
533 | test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary | ||
534 | putStrLn "------ read . show" | ||
535 | test (\m -> (m::RM) == read (show m)) | ||
536 | test (\m -> (m::CM) == read (show m)) | ||
537 | test (\m -> toRows (m::RM) == read (show (toRows m))) | ||
538 | test (\m -> toRows (m::CM) == read (show (toRows m))) | ||
539 | test (\m -> (m::FM) == read (show m)) | ||
540 | test (\m -> (m::ZM) == read (show m)) | ||
541 | test (\m -> toRows (m::FM) == read (show (toRows m))) | ||
542 | test (\m -> toRows (m::ZM) == read (show (toRows m))) | ||
543 | putStrLn "------ some unit tests" | ||
544 | _ <- runTestTT $ TestList | ||
545 | [ utest "1E5 rots" rotTest | ||
546 | , utest "det1" detTest1 | ||
547 | , utest "invlndet" detTest2 | ||
548 | , utest "expm1" (expmTest1) | ||
549 | , utest "expm2" (expmTest2) | ||
550 | , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) | ||
551 | , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM) | ||
552 | , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0 | ||
553 | , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] | ||
554 | -- , utest "gamma" (gamma 5 == 24.0) | ||
555 | -- , besselTest | ||
556 | -- , exponentialTest | ||
557 | , utest "deriv" derivTest | ||
558 | , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8) | ||
559 | , utest "polySolve" (polySolveProp [1,2,3,4]) | ||
560 | , minimizationTest | ||
561 | , rootFindingTest | ||
562 | , utest "randomGaussian" randomTestGaussian | ||
563 | , utest "randomUniform" randomTestUniform | ||
564 | , utest "buildVector/Matrix" $ | ||
565 | complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral | ||
566 | && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0) | ||
567 | , utest "rank" $ rank ((2><3)[1,0,0,1,6*eps,0]) == 1 | ||
568 | && rank ((2><3)[1,0,0,1,7*eps,0]) == 2 | ||
569 | , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM) | ||
570 | , odeTest | ||
571 | , fittingTest | ||
572 | , mbCholTest | ||
573 | , utest "offset" offsetTest | ||
574 | , normsVTest | ||
575 | , normsMTest | ||
576 | , sumprodTest | ||
577 | , chainTest | ||
578 | , succTest | ||
579 | , findAssocTest | ||
580 | , condTest | ||
581 | , conformTest | ||
582 | , accumTest | ||
583 | ] | ||
584 | return () | ||
585 | |||
586 | |||
587 | -- single precision approximate equality | ||
588 | infixl 4 |~~| | ||
589 | a |~~| b = a :~6~: b | ||
590 | |||
591 | makeUnitary v | realPart n > 1 = v / scalar n | ||
592 | | otherwise = v | ||
593 | where n = sqrt (conj v <.> v) | ||
594 | |||
595 | -- -- | Some additional tests on big matrices. They take a few minutes. | ||
596 | -- runBigTests :: IO () | ||
597 | -- runBigTests = undefined | ||
598 | |||
599 | -- testcase for nonempty fpu stack | ||
600 | findNaN :: Int -> Bool | ||
601 | findNaN n = all (bugProp . eye) (take n $ cycle [1..20]) | ||
602 | where eye m = ident m :: Matrix ( Double) | ||
603 | |||
604 | -------------------------------------------------------------------------------- | ||
605 | |||
606 | -- | Performance measurements. | ||
607 | runBenchmarks :: IO () | ||
608 | runBenchmarks = do | ||
609 | --cholBench | ||
610 | solveBench | ||
611 | subBench | ||
612 | multBench | ||
613 | svdBench | ||
614 | eigBench | ||
615 | putStrLn "" | ||
616 | |||
617 | -------------------------------- | ||
618 | |||
619 | time msg act = do | ||
620 | putStr (msg++" ") | ||
621 | t0 <- getCPUTime | ||
622 | act `seq` putStr " " | ||
623 | t1 <- getCPUTime | ||
624 | printf "%6.2f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO () | ||
625 | return () | ||
626 | |||
627 | -------------------------------- | ||
628 | |||
629 | manymult n = foldl1' (<>) (map rot2 angles) where | ||
630 | angles = toList $ linspace n (0,1) | ||
631 | rot2 :: Double -> Matrix Double | ||
632 | rot2 a = (3><3) [ c,0,s | ||
633 | , 0,1,0 | ||
634 | ,-s,0,c ] | ||
635 | where c = cos a | ||
636 | s = sin a | ||
637 | |||
638 | multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double)) | ||
639 | |||
640 | -------------------------------- | ||
641 | |||
642 | subBench = do | ||
643 | putStrLn "" | ||
644 | let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v)) | ||
645 | time "0.1M subVector " (g (constant 1 (1+10^5) :: Vector Double) @> 0) | ||
646 | let f = foldl1' (.) (replicate (10^5) (fromRows.toRows)) | ||
647 | time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0)) | ||
648 | time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0)) | ||
649 | |||
650 | -------------------------------- | ||
651 | |||
652 | multBench = do | ||
653 | let a = ident 1000 :: Matrix Double | ||
654 | let b = ident 2000 :: Matrix Double | ||
655 | a `seq` b `seq` putStrLn "" | ||
656 | time "product of 1M different 3x3 matrices" (manymult (10^6)) | ||
657 | putStrLn "" | ||
658 | time "product of 1M constant 1x1 matrices" (multb 1) | ||
659 | time "product of 1M constant 3x3 matrices" (multb 3) | ||
660 | --time "product of 1M constant 5x5 matrices" (multb 5) | ||
661 | time "product of 1M const. 10x10 matrices" (multb 10) | ||
662 | --time "product of 1M const. 15x15 matrices" (multb 15) | ||
663 | time "product of 1M const. 20x20 matrices" (multb 20) | ||
664 | --time "product of 1M const. 25x25 matrices" (multb 25) | ||
665 | putStrLn "" | ||
666 | time "product (1000 x 1000)<>(1000 x 1000)" (a<>a) | ||
667 | time "product (2000 x 2000)<>(2000 x 2000)" (b<>b) | ||
668 | |||
669 | -------------------------------- | ||
670 | |||
671 | eigBench = do | ||
672 | let m = reshape 1000 (randomVector 777 Uniform (1000*1000)) | ||
673 | s = m + trans m | ||
674 | m `seq` s `seq` putStrLn "" | ||
675 | time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m) | ||
676 | time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m) | ||
677 | time "eigenvalues general 1000x1000" (eigenvalues m) | ||
678 | time "eigenvectors general 1000x1000" (snd $ eig m) | ||
679 | |||
680 | -------------------------------- | ||
681 | |||
682 | svdBench = do | ||
683 | let a = reshape 500 (randomVector 777 Uniform (3000*500)) | ||
684 | b = reshape 1000 (randomVector 777 Uniform (1000*1000)) | ||
685 | fv (_,_,v) = v@@>(0,0) | ||
686 | a `seq` b `seq` putStrLn "" | ||
687 | time "singular values 3000x500" (singularValues a) | ||
688 | time "thin svd 3000x500" (fv $ thinSVD a) | ||
689 | time "full svd 3000x500" (fv $ svd a) | ||
690 | time "singular values 1000x1000" (singularValues b) | ||
691 | time "full svd 1000x1000" (fv $ svd b) | ||
692 | |||
693 | -------------------------------- | ||
694 | |||
695 | solveBenchN n = do | ||
696 | let x = uniformSample 777 (2*n) (replicate n (-1,1)) | ||
697 | a = trans x <> x | ||
698 | b = asColumn $ randomVector 666 Uniform n | ||
699 | a `seq` b `seq` putStrLn "" | ||
700 | time ("svd solve " ++ show n) (linearSolveSVD a b) | ||
701 | time (" ls solve " ++ show n) (linearSolveLS a b) | ||
702 | time (" solve " ++ show n) (linearSolve a b) | ||
703 | time ("cholSolve " ++ show n) (cholSolve (chol a) b) | ||
704 | |||
705 | solveBench = do | ||
706 | solveBenchN 500 | ||
707 | solveBenchN 1000 | ||
708 | -- solveBenchN 1500 | ||
709 | |||
710 | -------------------------------- | ||
711 | |||
712 | cholBenchN n = do | ||
713 | let x = uniformSample 777 (2*n) (replicate n (-1,1)) | ||
714 | a = trans x <> x | ||
715 | a `seq` putStrLn "" | ||
716 | time ("chol " ++ show n) (chol a) | ||
717 | |||
718 | cholBench = do | ||
719 | cholBenchN 1200 | ||
720 | cholBenchN 600 | ||
721 | cholBenchN 300 | ||
722 | -- cholBenchN 150 | ||
723 | -- cholBenchN 50 | ||
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs deleted file mode 100644 index 6dd9cfe..0000000 --- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs +++ /dev/null | |||
@@ -1,249 +0,0 @@ | |||
1 | {-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-} | ||
2 | {-# OPTIONS_GHC -fno-warn-unused-imports #-} | ||
3 | ----------------------------------------------------------------------------- | ||
4 | {- | | ||
5 | Module : Numeric.LinearAlgebra.Tests.Instances | ||
6 | Copyright : (c) Alberto Ruiz 2008 | ||
7 | License : GPL-style | ||
8 | |||
9 | Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
10 | Stability : provisional | ||
11 | Portability : portable | ||
12 | |||
13 | Arbitrary instances for vectors, matrices. | ||
14 | |||
15 | -} | ||
16 | |||
17 | module Numeric.LinearAlgebra.Tests.Instances( | ||
18 | Sq(..), rSq,cSq, | ||
19 | Rot(..), rRot,cRot, | ||
20 | Her(..), rHer,cHer, | ||
21 | WC(..), rWC,cWC, | ||
22 | SqWC(..), rSqWC, cSqWC, | ||
23 | PosDef(..), rPosDef, cPosDef, | ||
24 | Consistent(..), rConsist, cConsist, | ||
25 | RM,CM, rM,cM, | ||
26 | FM,ZM, fM,zM | ||
27 | ) where | ||
28 | |||
29 | import System.Random | ||
30 | |||
31 | import Numeric.LinearAlgebra | ||
32 | import Control.Monad(replicateM) | ||
33 | #include "quickCheckCompat.h" | ||
34 | |||
35 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
36 | shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] | ||
37 | shrinkListElementwise [] = [] | ||
38 | shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ] | ||
39 | ++ [ x:ys | ys <- shrinkListElementwise xs ] | ||
40 | |||
41 | shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)] | ||
42 | shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ] | ||
43 | #endif | ||
44 | |||
45 | #if MIN_VERSION_QuickCheck(2,1,1) | ||
46 | #else | ||
47 | instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where | ||
48 | arbitrary = do | ||
49 | re <- arbitrary | ||
50 | im <- arbitrary | ||
51 | return (re :+ im) | ||
52 | |||
53 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
54 | shrink (re :+ im) = | ||
55 | [ u :+ v | (u,v) <- shrinkPair (re,im) ] | ||
56 | #else | ||
57 | -- this has been moved to the 'Coarbitrary' class in QuickCheck 2 | ||
58 | coarbitrary = undefined | ||
59 | #endif | ||
60 | |||
61 | #endif | ||
62 | |||
63 | chooseDim = sized $ \m -> choose (1,max 1 m) | ||
64 | |||
65 | instance (Field a, Arbitrary a) => Arbitrary (Vector a) where | ||
66 | arbitrary = do m <- chooseDim | ||
67 | l <- vector m | ||
68 | return $ fromList l | ||
69 | |||
70 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
71 | -- shrink any one of the components | ||
72 | shrink = map fromList . shrinkListElementwise . toList | ||
73 | |||
74 | #else | ||
75 | coarbitrary = undefined | ||
76 | #endif | ||
77 | |||
78 | instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where | ||
79 | arbitrary = do | ||
80 | m <- chooseDim | ||
81 | n <- chooseDim | ||
82 | l <- vector (m*n) | ||
83 | return $ (m><n) l | ||
84 | |||
85 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
86 | -- shrink any one of the components | ||
87 | shrink a = map (rows a >< cols a) | ||
88 | . shrinkListElementwise | ||
89 | . concat . toLists | ||
90 | $ a | ||
91 | #else | ||
92 | coarbitrary = undefined | ||
93 | #endif | ||
94 | |||
95 | |||
96 | -- a square matrix | ||
97 | newtype (Sq a) = Sq (Matrix a) deriving Show | ||
98 | instance (Element a, Arbitrary a) => Arbitrary (Sq a) where | ||
99 | arbitrary = do | ||
100 | n <- chooseDim | ||
101 | l <- vector (n*n) | ||
102 | return $ Sq $ (n><n) l | ||
103 | |||
104 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
105 | shrink (Sq a) = [ Sq b | b <- shrink a ] | ||
106 | #else | ||
107 | coarbitrary = undefined | ||
108 | #endif | ||
109 | |||
110 | |||
111 | -- a unitary matrix | ||
112 | newtype (Rot a) = Rot (Matrix a) deriving Show | ||
113 | instance (Field a, Arbitrary a) => Arbitrary (Rot a) where | ||
114 | arbitrary = do | ||
115 | Sq m <- arbitrary | ||
116 | let (q,_) = qr m | ||
117 | return (Rot q) | ||
118 | |||
119 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
120 | #else | ||
121 | coarbitrary = undefined | ||
122 | #endif | ||
123 | |||
124 | |||
125 | -- a complex hermitian or real symmetric matrix | ||
126 | newtype (Her a) = Her (Matrix a) deriving Show | ||
127 | instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where | ||
128 | arbitrary = do | ||
129 | Sq m <- arbitrary | ||
130 | let m' = m/2 | ||
131 | return $ Her (m' + ctrans m') | ||
132 | |||
133 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
134 | #else | ||
135 | coarbitrary = undefined | ||
136 | #endif | ||
137 | |||
138 | class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a | ||
139 | instance ArbitraryField Double | ||
140 | instance ArbitraryField (Complex Double) | ||
141 | |||
142 | |||
143 | -- a well-conditioned general matrix (the singular values are between 1 and 100) | ||
144 | newtype (WC a) = WC (Matrix a) deriving Show | ||
145 | instance (ArbitraryField a) => Arbitrary (WC a) where | ||
146 | arbitrary = do | ||
147 | m <- arbitrary | ||
148 | let (u,_,v) = svd m | ||
149 | r = rows m | ||
150 | c = cols m | ||
151 | n = min r c | ||
152 | sv' <- replicateM n (choose (1,100)) | ||
153 | let s = diagRect 0 (fromList sv') r c | ||
154 | return $ WC (u <> real s <> trans v) | ||
155 | |||
156 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
157 | #else | ||
158 | coarbitrary = undefined | ||
159 | #endif | ||
160 | |||
161 | |||
162 | -- a well-conditioned square matrix (the singular values are between 1 and 100) | ||
163 | newtype (SqWC a) = SqWC (Matrix a) deriving Show | ||
164 | instance (ArbitraryField a) => Arbitrary (SqWC a) where | ||
165 | arbitrary = do | ||
166 | Sq m <- arbitrary | ||
167 | let (u,_,v) = svd m | ||
168 | n = rows m | ||
169 | sv' <- replicateM n (choose (1,100)) | ||
170 | let s = diag (fromList sv') | ||
171 | return $ SqWC (u <> real s <> trans v) | ||
172 | |||
173 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
174 | #else | ||
175 | coarbitrary = undefined | ||
176 | #endif | ||
177 | |||
178 | |||
179 | -- a positive definite square matrix (the eigenvalues are between 0 and 100) | ||
180 | newtype (PosDef a) = PosDef (Matrix a) deriving Show | ||
181 | instance (ArbitraryField a, Num (Vector a)) | ||
182 | => Arbitrary (PosDef a) where | ||
183 | arbitrary = do | ||
184 | Her m <- arbitrary | ||
185 | let (_,v) = eigSH m | ||
186 | n = rows m | ||
187 | l <- replicateM n (choose (0,100)) | ||
188 | let s = diag (fromList l) | ||
189 | p = v <> real s <> ctrans v | ||
190 | return $ PosDef (0.5 * p + 0.5 * ctrans p) | ||
191 | |||
192 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
193 | #else | ||
194 | coarbitrary = undefined | ||
195 | #endif | ||
196 | |||
197 | |||
198 | -- a pair of matrices that can be multiplied | ||
199 | newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show | ||
200 | instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where | ||
201 | arbitrary = do | ||
202 | n <- chooseDim | ||
203 | k <- chooseDim | ||
204 | m <- chooseDim | ||
205 | la <- vector (n*k) | ||
206 | lb <- vector (k*m) | ||
207 | return $ Consistent ((n><k) la, (k><m) lb) | ||
208 | |||
209 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
210 | shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ] | ||
211 | #else | ||
212 | coarbitrary = undefined | ||
213 | #endif | ||
214 | |||
215 | |||
216 | |||
217 | type RM = Matrix Double | ||
218 | type CM = Matrix (Complex Double) | ||
219 | type FM = Matrix Float | ||
220 | type ZM = Matrix (Complex Float) | ||
221 | |||
222 | |||
223 | rM m = m :: RM | ||
224 | cM m = m :: CM | ||
225 | fM m = m :: FM | ||
226 | zM m = m :: ZM | ||
227 | |||
228 | |||
229 | rHer (Her m) = m :: RM | ||
230 | cHer (Her m) = m :: CM | ||
231 | |||
232 | rRot (Rot m) = m :: RM | ||
233 | cRot (Rot m) = m :: CM | ||
234 | |||
235 | rSq (Sq m) = m :: RM | ||
236 | cSq (Sq m) = m :: CM | ||
237 | |||
238 | rWC (WC m) = m :: RM | ||
239 | cWC (WC m) = m :: CM | ||
240 | |||
241 | rSqWC (SqWC m) = m :: RM | ||
242 | cSqWC (SqWC m) = m :: CM | ||
243 | |||
244 | rPosDef (PosDef m) = m :: RM | ||
245 | cPosDef (PosDef m) = m :: CM | ||
246 | |||
247 | rConsist (Consistent (a,b)) = (a,b::RM) | ||
248 | cConsist (Consistent (a,b)) = (a,b::CM) | ||
249 | |||
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs deleted file mode 100644 index fe13544..0000000 --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ /dev/null | |||
@@ -1,272 +0,0 @@ | |||
1 | {-# LANGUAGE CPP, FlexibleContexts #-} | ||
2 | {-# OPTIONS_GHC -fno-warn-unused-imports #-} | ||
3 | ----------------------------------------------------------------------------- | ||
4 | {- | | ||
5 | Module : Numeric.LinearAlgebra.Tests.Properties | ||
6 | Copyright : (c) Alberto Ruiz 2008 | ||
7 | License : GPL-style | ||
8 | |||
9 | Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
10 | Stability : provisional | ||
11 | Portability : portable | ||
12 | |||
13 | Testing properties. | ||
14 | |||
15 | -} | ||
16 | |||
17 | module Numeric.LinearAlgebra.Tests.Properties ( | ||
18 | dist, (|~|), (~:), Aprox((:~)), | ||
19 | zeros, ones, | ||
20 | square, | ||
21 | unitary, | ||
22 | hermitian, | ||
23 | wellCond, | ||
24 | positiveDefinite, | ||
25 | upperTriang, | ||
26 | upperHessenberg, | ||
27 | luProp, | ||
28 | invProp, | ||
29 | pinvProp, | ||
30 | detProp, | ||
31 | nullspaceProp, | ||
32 | bugProp, | ||
33 | svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, | ||
34 | svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, | ||
35 | eigProp, eigSHProp, eigProp2, eigSHProp2, | ||
36 | qrProp, rqProp, rqProp1, rqProp2, rqProp3, | ||
37 | hessProp, | ||
38 | schurProp1, schurProp2, | ||
39 | cholProp, exactProp, | ||
40 | expmDiagProp, | ||
41 | multProp1, multProp2, | ||
42 | subProp, | ||
43 | linearSolveProp, linearSolveProp2 | ||
44 | ) where | ||
45 | |||
46 | import Numeric.LinearAlgebra --hiding (real,complex) | ||
47 | import Numeric.LinearAlgebra.LAPACK | ||
48 | import Debug.Trace | ||
49 | #include "quickCheckCompat.h" | ||
50 | |||
51 | |||
52 | --real x = real'' x | ||
53 | --complex x = complex'' x | ||
54 | |||
55 | debug x = trace (show x) x | ||
56 | |||
57 | -- relative error | ||
58 | dist :: (Normed c t, Num (c t)) => c t -> c t -> Double | ||
59 | dist a b = realToFrac r | ||
60 | where norm = pnorm Infinity | ||
61 | na = norm a | ||
62 | nb = norm b | ||
63 | nab = norm (a-b) | ||
64 | mx = max na nb | ||
65 | mn = min na nb | ||
66 | r = if mn < peps | ||
67 | then mx | ||
68 | else nab/mx | ||
69 | |||
70 | infixl 4 |~| | ||
71 | a |~| b = a :~10~: b | ||
72 | --a |~| b = dist a b < 10^^(-10) | ||
73 | |||
74 | data Aprox a = (:~) a Int | ||
75 | -- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool | ||
76 | a :~n~: b = dist a b < 10^^(-n) | ||
77 | |||
78 | ------------------------------------------------------ | ||
79 | |||
80 | square m = rows m == cols m | ||
81 | |||
82 | -- orthonormal columns | ||
83 | orthonormal m = ctrans m <> m |~| ident (cols m) | ||
84 | |||
85 | unitary m = square m && orthonormal m | ||
86 | |||
87 | hermitian m = square m && m |~| ctrans m | ||
88 | |||
89 | wellCond m = rcond m > 1/100 | ||
90 | |||
91 | positiveDefinite m = minimum (toList e) > 0 | ||
92 | where (e,_v) = eigSH m | ||
93 | |||
94 | upperTriang m = rows m == 1 || down == z | ||
95 | where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) | ||
96 | z = constant 0 (dim down) | ||
97 | |||
98 | upperHessenberg m = rows m < 3 || down == z | ||
99 | where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) | ||
100 | z = constant 0 (dim down) | ||
101 | |||
102 | zeros (r,c) = reshape c (constant 0 (r*c)) | ||
103 | |||
104 | ones (r,c) = zeros (r,c) + 1 | ||
105 | |||
106 | ----------------------------------------------------- | ||
107 | |||
108 | luProp m = m |~| p <> l <> u && f (det p) |~| f s | ||
109 | where (l,u,p,s) = lu m | ||
110 | f x = fromList [x] | ||
111 | |||
112 | invProp m = m <> inv m |~| ident (rows m) | ||
113 | |||
114 | pinvProp m = m <> p <> m |~| m | ||
115 | && p <> m <> p |~| p | ||
116 | && hermitian (m<>p) | ||
117 | && hermitian (p<>m) | ||
118 | where p = pinv m | ||
119 | |||
120 | detProp m = s d1 |~| s d2 | ||
121 | where d1 = det m | ||
122 | d2 = det' * det q | ||
123 | det' = product $ toList $ takeDiag r | ||
124 | (q,r) = qr m | ||
125 | s x = fromList [x] | ||
126 | |||
127 | nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) | ||
128 | && orthonormal (fromColumns nl)) | ||
129 | where nl = nullspacePrec 1 m | ||
130 | n = fromColumns nl | ||
131 | r = rows m | ||
132 | c = cols m - rank m | ||
133 | |||
134 | ------------------------------------------------------------------ | ||
135 | |||
136 | -- testcase for nonempty fpu stack | ||
137 | -- uncommenting unitary' signature eliminates the problem | ||
138 | bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v | ||
139 | where (u,d,v) = fullSVD m | ||
140 | -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool | ||
141 | unitary' a = unitary a | ||
142 | |||
143 | ------------------------------------------------------------------ | ||
144 | |||
145 | -- fullSVD | ||
146 | svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v | ||
147 | where (u,d,v) = fullSVD m | ||
148 | |||
149 | svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where | ||
150 | (u,s,v) = svdfun m | ||
151 | d = diagRect 0 s (rows m) (cols m) | ||
152 | |||
153 | svdProp1b svdfun m = unitary u && unitary v where | ||
154 | (u,_,v) = svdfun m | ||
155 | |||
156 | -- thinSVD | ||
157 | svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) | ||
158 | where (u,s,v) = thinSVDfun m | ||
159 | |||
160 | -- compactSVD | ||
161 | svdProp3 m = (m |~| u <> real (diag s) <> trans v | ||
162 | && orthonormal u && orthonormal v) | ||
163 | where (u,s,v) = compactSVD m | ||
164 | |||
165 | svdProp4 m' = m |~| u <> real (diag s) <> trans v | ||
166 | && orthonormal u && orthonormal v | ||
167 | && (dim s == r || r == 0 && dim s == 1) | ||
168 | where (u,s,v) = compactSVD m | ||
169 | m = fromBlocks [[m'],[m']] | ||
170 | r = rank m' | ||
171 | |||
172 | svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where | ||
173 | s1 = svR m | ||
174 | s2 = svRd m | ||
175 | (_,s3,_) = svdR m | ||
176 | (_,s4,_) = svdRd m | ||
177 | (_,s5,_) = thinSVDR m | ||
178 | (_,s6,_) = thinSVDRd m | ||
179 | |||
180 | svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where | ||
181 | s1 = svC m | ||
182 | s2 = svCd m | ||
183 | (_,s3,_) = svdC m | ||
184 | (_,s4,_) = svdCd m | ||
185 | (_,s5,_) = thinSVDC m | ||
186 | (_,s6,_) = thinSVDCd m | ||
187 | |||
188 | svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' | ||
189 | where (u,s,v) = svdR m | ||
190 | (s',v') = rightSVR m | ||
191 | (u',s'') = leftSVR m | ||
192 | |||
193 | svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' | ||
194 | where (u,s,v) = svdC m | ||
195 | (s',v') = rightSVC m | ||
196 | (u',s'') = leftSVC m | ||
197 | |||
198 | svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' | ||
199 | where (u,s,v) = svd m | ||
200 | (s',v') = rightSV m | ||
201 | (u',_s'') = leftSV m | ||
202 | s''' = singularValues m | ||
203 | |||
204 | ------------------------------------------------------------------ | ||
205 | |||
206 | eigProp m = complex m <> v |~| v <> diag s | ||
207 | where (s, v) = eig m | ||
208 | |||
209 | eigSHProp m = m <> v |~| v <> real (diag s) | ||
210 | && unitary v | ||
211 | && m |~| v <> real (diag s) <> ctrans v | ||
212 | where (s, v) = eigSH m | ||
213 | |||
214 | eigProp2 m = fst (eig m) |~| eigenvalues m | ||
215 | |||
216 | eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m | ||
217 | |||
218 | ------------------------------------------------------------------ | ||
219 | |||
220 | qrProp m = q <> r |~| m && unitary q && upperTriang r | ||
221 | where (q,r) = qr m | ||
222 | |||
223 | rqProp m = r <> q |~| m && unitary q && upperTriang' r | ||
224 | where (r,q) = rq m | ||
225 | |||
226 | rqProp1 m = r <> q |~| m | ||
227 | where (r,q) = rq m | ||
228 | |||
229 | rqProp2 m = unitary q | ||
230 | where (r,q) = rq m | ||
231 | |||
232 | rqProp3 m = upperTriang' r | ||
233 | where (r,q) = rq m | ||
234 | |||
235 | upperTriang' r = upptr (rows r) (cols r) * r |~| r | ||
236 | where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1 | ||
237 | where t = f-c | ||
238 | |||
239 | hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h | ||
240 | where (p,h) = hess m | ||
241 | |||
242 | schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s | ||
243 | where (u,s) = schur m | ||
244 | |||
245 | schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme | ||
246 | where (u,s) = schur m | ||
247 | |||
248 | cholProp m = m |~| ctrans c <> c && upperTriang c | ||
249 | where c = chol m | ||
250 | |||
251 | exactProp m = chol m == chol (m+0) | ||
252 | |||
253 | expmDiagProp m = expm (logm m) :~ 7 ~: complex m | ||
254 | where logm = matFunc log | ||
255 | |||
256 | -- reference multiply | ||
257 | mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ] | ||
258 | where doth u v = sum $ zipWith (*) (toList u) (toList v) | ||
259 | |||
260 | multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) | ||
261 | |||
262 | multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) | ||
263 | |||
264 | linearSolveProp f m = f m m |~| ident (rows m) | ||
265 | |||
266 | linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) | ||
267 | where q = min (rows a) (cols a) | ||
268 | b = a <> x | ||
269 | wc = rank a == q | ||
270 | |||
271 | subProp m = m == (trans . fromColumns . toRows) m | ||
272 | |||
diff --git a/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h b/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h deleted file mode 100644 index 714587b..0000000 --- a/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h +++ /dev/null | |||
@@ -1,33 +0,0 @@ | |||
1 | #ifndef MIN_VERSION_QuickCheck | ||
2 | #define MIN_VERSION_QuickCheck(A,B,C) 1 | ||
3 | #endif | ||
4 | |||
5 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
6 | import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector | ||
7 | ,sized,classify,Testable,Property | ||
8 | |||
9 | ,quickCheckWith,maxSize,stdArgs,shrink) | ||
10 | |||
11 | #else | ||
12 | import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector | ||
13 | ,sized,classify,Testable,Property | ||
14 | |||
15 | ,check,configSize,defaultConfig,trivial) | ||
16 | #endif | ||
17 | |||
18 | |||
19 | |||
20 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
21 | trivial :: Testable a => Bool -> a -> Property | ||
22 | trivial = (`classify` "trivial") | ||
23 | #else | ||
24 | #endif | ||
25 | |||
26 | |||
27 | -- define qCheck, which used to be in Tests.hs | ||
28 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
29 | qCheck n = quickCheckWith stdArgs {maxSize = n} | ||
30 | #else | ||
31 | qCheck n = check defaultConfig {configSize = const n} | ||
32 | #endif | ||
33 | |||