summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra/Tests.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/LinearAlgebra/Tests.hs')
-rw-r--r--lib/Numeric/LinearAlgebra/Tests.hs723
1 files changed, 0 insertions, 723 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{- |
5Module : Numeric.LinearAlgebra.Tests
6Copyright : (c) Alberto Ruiz 2007-9
7License : GPL-style
8
9Maintainer : Alberto Ruiz (aruiz at um dot es)
10Stability : provisional
11Portability : portable
12
13Some tests.
14
15-}
16
17module 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
24import Data.Packed.Random
25import Numeric.LinearAlgebra
26import Numeric.LinearAlgebra.LAPACK
27import Numeric.LinearAlgebra.Tests.Instances
28import Numeric.LinearAlgebra.Tests.Properties
29import Test.HUnit hiding ((~:),test,Testable,State)
30import System.Info
31import Data.List(foldl1')
32import Numeric.GSL
33import Prelude hiding ((^))
34import qualified Prelude
35import System.CPUTime
36import Text.Printf
37import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr)
38import Control.Arrow((***))
39import Debug.Trace
40
41#include "Tests/quickCheckCompat.h"
42
43debug x = trace (show x) x
44
45a ^ b = a Prelude.^ (b :: Int)
46
47utest str b = TestCase $ assertBool str b
48
49a ~~ b = fromList a |~| fromList b
50
51feye n = flipud (ident n) :: Matrix Double
52
53-----------------------------------------------------------
54
55detTest1 = 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
70detTest2 = 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
80polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
81
82polySolveProp p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
83
84---------------------------------------------------------------------
85
86quad f a b = fst $ integrateQAGS 1E-9 100 f a b
87
88-- A multiple integral can be easily defined using partial application
89quad2 f a b g1 g2 = quad h a b
90 where h x = quad (f x) (g1 x) (g2 x)
91
92volSphere 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
97derivTest = 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
112nd1 = (3><3) [ 1/2, 1/4, 1/4
113 , 0/1, 1/2, 1/4
114 , 1/2, 1/4, 1/2 :: Double]
115
116nd2 = (2><2) [1, 0, 1, 1:: Complex Double]
117
118expmTest1 = 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
129expmTest2 = expm nd2 :~15~: (2><2)
130 [ 2.718281828459045
131 , 0.000000000000000
132 , 2.718281828459045
133 , 2.718281828459045 ]
134
135---------------------------------------------------------------------
136
137minimizationTest = 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
148rootFindingTest = 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
159odeTest = 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
167fittingTest = 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
187mbCholTest = 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
195randomTestGaussian = 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
203randomTestUniform = 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
209rot :: Double -> Matrix Double
210rot a = (3><3) [ c,0,s
211 , 0,1,0
212 ,-s,0,c ]
213 where c = cos a
214 s = sin a
215
216rotTest = 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
224offsetTest = 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
232normsVTest = 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
257normsMTest = 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
285sumprodTest = 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
301chainTest = 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
310conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m)
311
312---------------------------------------------------------------------
313
314newtype State s a = State { runState :: s -> (a,s) }
315
316instance 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
321state_get :: State s s
322state_get = State $ \s -> (s,s)
323
324state_put :: s -> State s ()
325state_put s = State $ \_ -> ((),s)
326
327evalState :: State s a -> s -> a
328evalState m s = let (a,s') = runState m s
329 in seq s' a
330
331newtype MaybeT m a = MaybeT { runMaybeT :: m (Maybe a) }
332
333instance 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
342lift_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
348successive_ 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
357successive 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
364succTest = 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
371findAssocTest = 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
379condTest = 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
386conformTest = 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
398accumTest = 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).
412runTests :: Int -- ^ maximum dimension
413 -> IO ()
414runTests 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
588infixl 4 |~~|
589a |~~| b = a :~6~: b
590
591makeUnitary 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
600findNaN :: Int -> Bool
601findNaN n = all (bugProp . eye) (take n $ cycle [1..20])
602 where eye m = ident m :: Matrix ( Double)
603
604--------------------------------------------------------------------------------
605
606-- | Performance measurements.
607runBenchmarks :: IO ()
608runBenchmarks = do
609 --cholBench
610 solveBench
611 subBench
612 multBench
613 svdBench
614 eigBench
615 putStrLn ""
616
617--------------------------------
618
619time 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
629manymult 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
638multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double))
639
640--------------------------------
641
642subBench = 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
652multBench = 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
671eigBench = 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
682svdBench = 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
695solveBenchN 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
705solveBench = do
706 solveBenchN 500
707 solveBenchN 1000
708 -- solveBenchN 1500
709
710--------------------------------
711
712cholBenchN 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
718cholBench = do
719 cholBenchN 1200
720 cholBenchN 600
721 cholBenchN 300
722-- cholBenchN 150
723-- cholBenchN 50
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-07 07:12:19 +0000