summaryrefslogtreecommitdiff
path: root/lib/Numeric
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2011-12-14 13:08:43 +0100
committerAlberto Ruiz <aruiz@um.es>2011-12-14 13:08:43 +0100
commit77552d080e88fc70312f55fd3303fac3464ab46e (patch)
tree1dc87dd22ce0da0f1807765568fbc04285bf3621 /lib/Numeric
parentc3bda2d38c432fb53ce456cba295b097fd4d6ad1 (diff)
new package hmatrix-tests
Diffstat (limited to 'lib/Numeric')
-rw-r--r--lib/Numeric/LinearAlgebra/Tests.hs723
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Instances.hs249
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Properties.hs272
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h33
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{- |
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
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{- |
5Module : Numeric.LinearAlgebra.Tests.Instances
6Copyright : (c) Alberto Ruiz 2008
7License : GPL-style
8
9Maintainer : Alberto Ruiz (aruiz at um dot es)
10Stability : provisional
11Portability : portable
12
13Arbitrary instances for vectors, matrices.
14
15-}
16
17module 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
29import System.Random
30
31import Numeric.LinearAlgebra
32import Control.Monad(replicateM)
33#include "quickCheckCompat.h"
34
35#if MIN_VERSION_QuickCheck(2,0,0)
36shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]]
37shrinkListElementwise [] = []
38shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ]
39 ++ [ x:ys | ys <- shrinkListElementwise xs ]
40
41shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)]
42shrinkPair (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
47instance (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
63chooseDim = sized $ \m -> choose (1,max 1 m)
64
65instance (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
78instance (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
97newtype (Sq a) = Sq (Matrix a) deriving Show
98instance (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
112newtype (Rot a) = Rot (Matrix a) deriving Show
113instance (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
126newtype (Her a) = Her (Matrix a) deriving Show
127instance (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
138class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
139instance ArbitraryField Double
140instance ArbitraryField (Complex Double)
141
142
143-- a well-conditioned general matrix (the singular values are between 1 and 100)
144newtype (WC a) = WC (Matrix a) deriving Show
145instance (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)
163newtype (SqWC a) = SqWC (Matrix a) deriving Show
164instance (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)
180newtype (PosDef a) = PosDef (Matrix a) deriving Show
181instance (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
199newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show
200instance (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
217type RM = Matrix Double
218type CM = Matrix (Complex Double)
219type FM = Matrix Float
220type ZM = Matrix (Complex Float)
221
222
223rM m = m :: RM
224cM m = m :: CM
225fM m = m :: FM
226zM m = m :: ZM
227
228
229rHer (Her m) = m :: RM
230cHer (Her m) = m :: CM
231
232rRot (Rot m) = m :: RM
233cRot (Rot m) = m :: CM
234
235rSq (Sq m) = m :: RM
236cSq (Sq m) = m :: CM
237
238rWC (WC m) = m :: RM
239cWC (WC m) = m :: CM
240
241rSqWC (SqWC m) = m :: RM
242cSqWC (SqWC m) = m :: CM
243
244rPosDef (PosDef m) = m :: RM
245cPosDef (PosDef m) = m :: CM
246
247rConsist (Consistent (a,b)) = (a,b::RM)
248cConsist (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{- |
5Module : Numeric.LinearAlgebra.Tests.Properties
6Copyright : (c) Alberto Ruiz 2008
7License : GPL-style
8
9Maintainer : Alberto Ruiz (aruiz at um dot es)
10Stability : provisional
11Portability : portable
12
13Testing properties.
14
15-}
16
17module 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
46import Numeric.LinearAlgebra --hiding (real,complex)
47import Numeric.LinearAlgebra.LAPACK
48import Debug.Trace
49#include "quickCheckCompat.h"
50
51
52--real x = real'' x
53--complex x = complex'' x
54
55debug x = trace (show x) x
56
57-- relative error
58dist :: (Normed c t, Num (c t)) => c t -> c t -> Double
59dist 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
70infixl 4 |~|
71a |~| b = a :~10~: b
72--a |~| b = dist a b < 10^^(-10)
73
74data Aprox a = (:~) a Int
75-- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
76a :~n~: b = dist a b < 10^^(-n)
77
78------------------------------------------------------
79
80square m = rows m == cols m
81
82-- orthonormal columns
83orthonormal m = ctrans m <> m |~| ident (cols m)
84
85unitary m = square m && orthonormal m
86
87hermitian m = square m && m |~| ctrans m
88
89wellCond m = rcond m > 1/100
90
91positiveDefinite m = minimum (toList e) > 0
92 where (e,_v) = eigSH m
93
94upperTriang m = rows m == 1 || down == z
95 where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
96 z = constant 0 (dim down)
97
98upperHessenberg m = rows m < 3 || down == z
99 where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))
100 z = constant 0 (dim down)
101
102zeros (r,c) = reshape c (constant 0 (r*c))
103
104ones (r,c) = zeros (r,c) + 1
105
106-----------------------------------------------------
107
108luProp m = m |~| p <> l <> u && f (det p) |~| f s
109 where (l,u,p,s) = lu m
110 f x = fromList [x]
111
112invProp m = m <> inv m |~| ident (rows m)
113
114pinvProp m = m <> p <> m |~| m
115 && p <> m <> p |~| p
116 && hermitian (m<>p)
117 && hermitian (p<>m)
118 where p = pinv m
119
120detProp 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
127nullspaceProp 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
138bugProp 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
146svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v
147 where (u,d,v) = fullSVD m
148
149svdProp1a 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
153svdProp1b svdfun m = unitary u && unitary v where
154 (u,_,v) = svdfun m
155
156-- thinSVD
157svdProp2 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
161svdProp3 m = (m |~| u <> real (diag s) <> trans v
162 && orthonormal u && orthonormal v)
163 where (u,s,v) = compactSVD m
164
165svdProp4 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
172svdProp5a 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
180svdProp5b 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
188svdProp6a 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
193svdProp6b 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
198svdProp7 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
206eigProp m = complex m <> v |~| v <> diag s
207 where (s, v) = eig m
208
209eigSHProp 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
214eigProp2 m = fst (eig m) |~| eigenvalues m
215
216eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
217
218------------------------------------------------------------------
219
220qrProp m = q <> r |~| m && unitary q && upperTriang r
221 where (q,r) = qr m
222
223rqProp m = r <> q |~| m && unitary q && upperTriang' r
224 where (r,q) = rq m
225
226rqProp1 m = r <> q |~| m
227 where (r,q) = rq m
228
229rqProp2 m = unitary q
230 where (r,q) = rq m
231
232rqProp3 m = upperTriang' r
233 where (r,q) = rq m
234
235upperTriang' 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
239hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
240 where (p,h) = hess m
241
242schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s
243 where (u,s) = schur m
244
245schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme
246 where (u,s) = schur m
247
248cholProp m = m |~| ctrans c <> c && upperTriang c
249 where c = chol m
250
251exactProp m = chol m == chol (m+0)
252
253expmDiagProp m = expm (logm m) :~ 7 ~: complex m
254 where logm = matFunc log
255
256-- reference multiply
257mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
258 where doth u v = sum $ zipWith (*) (toList u) (toList v)
259
260multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)
261
262multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a)
263
264linearSolveProp f m = f m m |~| ident (rows m)
265
266linearSolveProp2 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
271subProp 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)
6import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
7 ,sized,classify,Testable,Property
8
9 ,quickCheckWith,maxSize,stdArgs,shrink)
10
11#else
12import 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)
21trivial :: Testable a => Bool -> a -> Property
22trivial = (`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)
29qCheck n = quickCheckWith stdArgs {maxSize = n}
30#else
31qCheck n = check defaultConfig {configSize = const n}
32#endif
33