hmatrix-tests 0.4.1.0 → 0.5.0.0

raw patch · 5 files changed

+473/−293 lines, 5 filesdep +deepseqdep ~hmatrixdep ~hmatrix-gslPVP ok

version bump matches the API change (PVP)

Dependencies added: deepseq

Dependency ranges changed: hmatrix, hmatrix-gsl

API changes (from Hackage documentation)

Files

• hmatrix-tests.cabal +4/−4
• src/Numeric/GSL/Tests.hs +59/−3
• src/Numeric/LinearAlgebra/Tests.hs +297/−122
• src/Numeric/LinearAlgebra/Tests/Instances.hs +22/−87
• src/Numeric/LinearAlgebra/Tests/Properties.hs +91/−77

@@ -1,5 +1,5 @@ Name:               hmatrix-tests-Version:            0.4.1.0+Version:            0.5.0.0 License:            BSD3 License-file:       LICENSE Author:             Alberto Ruiz@@ -26,11 +26,11 @@  library -    Build-Depends:      base >= 4 && < 5,+    Build-Depends:      base >= 4 && < 5, deepseq,                         QuickCheck >= 2, HUnit, random,-                        hmatrix >= 0.16+                        hmatrix >= 0.17     if flag(gsl)-      Build-Depends:    hmatrix-gsl >= 0.16+      Build-Depends:    hmatrix-gsl >= 0.17      hs-source-dirs:     src 

@@ -19,10 +19,11 @@  import Test.HUnit (runTestTT, failures, Test(..), errors) -import Numeric.LinearAlgebra+import Numeric.LinearAlgebra.HMatrix import Numeric.GSL+import Numeric.GSL.SimulatedAnnealing import Numeric.LinearAlgebra.Tests (qCheck, utest)-import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~))+import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~), (~=))  --------------------------------------------------------------------- @@ -42,7 +43,7 @@     sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]      ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d-    ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5+    ok2 = norm_2 (fromList (map fst sols) - fromList sol) < 1E-5  --------------------------------------------------------------------- @@ -66,8 +67,61 @@           jacobian a b [x,_y] = [ [-a    , 0]                                 , [-2*b*x, b] ] +--------------------------------------------------------------------++interpolationTest = TestList [+    utest "interpolation evaluateV" (esol ~= ev)+  , utest "interpolation evaluate" (esol ~= eval)+  , utest "interpolation evaluateDerivativeV" (desol ~= dev)+  , utest "interpolation evaluateDerivative" (desol ~= de)+  , utest "interpolation evaluateDerivative2V" (d2esol ~= d2ev)+  , utest "interpolation evaluateDerivative2" (d2esol ~= d2e)+  , utest "interpolation evaluateIntegralV" (intesol ~= intev)+  , utest "interpolation evaluateIntegral" (intesol ~= inte)+  ]+  where+    xtest = 2.2+    applyVec f = f Akima xs ys xtest+    applyList f = f Akima (zip xs' ys') xtest++    esol = xtest**2+    ev = applyVec evaluateV+    eval = applyList evaluate++    desol = 2*xtest+    dev = applyVec evaluateDerivativeV+    de = applyList evaluateDerivative++    d2esol = 2+    d2ev = applyVec evaluateDerivative2V+    d2e = applyList evaluateDerivative2++    intesol = 1/3 * xtest**3+    intev = evaluateIntegralV Akima xs ys 0 xtest+    inte = evaluateIntegral Akima (zip xs' ys') (0, xtest)++    xs' = [-1..10]+    ys' = map (**2) xs'+    xs = vector xs'+    ys = vector ys'+ --------------------------------------------------------------------- +simanTest = TestList [+  -- We use a slightly more relaxed tolerance here because the+  -- simulated annealer is randomized+  utest "simulated annealing manual example" $ abs (result - 1.3631300) < 1e-6+  ]+  where+    -- This is the example from the GSL manual.+    result = simanSolve 0 1 exampleParams 15.5 exampleE exampleM exampleS Nothing+    exampleParams = SimulatedAnnealingParams 200 10000 1.0 1.0 0.008 1.003 2.0e-6+    exampleE x = exp (-(x - 1)**2) * sin (8 * x)+    exampleM x y = abs $ x - y+    exampleS rands stepSize current = (rands ! 0) * 2 * stepSize - stepSize + current++---------------------------------------------------------------------+ minimizationTest = TestList     [ utest "minimization conjugatefr" (minim1 f df [5,7] ~~ [1,2])     , utest "minimization nmsimplex2"  (minim2 f [5,7] `elem` [24,25])@@ -123,6 +177,8 @@         , odeTest         , rootFindingTest         , minimizationTest+        , interpolationTest+        , simanTest         , utest "deriv" derivTest         , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5**3) < 1E-8)         , utest "polySolve" (polySolveProp [1,2,3,4])

@@ -1,6 +1,11 @@ {-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-} {-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}  ----------------------------------------------------------------------------- {- |@@ -25,12 +30,9 @@ --, runBigTests ) where -import Numeric.LinearAlgebra-import Numeric.LinearAlgebra.HMatrix hiding ((<>),linearSolve)+import Numeric.LinearAlgebra hiding (unitary)+import Numeric.LinearAlgebra.Devel import Numeric.LinearAlgebra.Static(L)-import Numeric.LinearAlgebra.Util(col,row)-import Data.Packed-import Numeric.LinearAlgebra.LAPACK import Numeric.LinearAlgebra.Tests.Instances import Numeric.LinearAlgebra.Tests.Properties import Test.HUnit hiding ((~:),test,Testable,State)@@ -41,19 +43,18 @@ import System.CPUTime import System.Exit import Text.Printf-import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr)+import Numeric.LinearAlgebra.Devel(unsafeFromForeignPtr,unsafeToForeignPtr) import Control.Arrow((***)) import Debug.Trace import Control.Monad(when)-import Numeric.LinearAlgebra.Util hiding (ones,row,col) import Control.Applicative import Control.Monad(ap)--import Data.Packed.ST+import Control.DeepSeq ( NFData(..) )  import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector                       ,sized,classify,Testable,Property                       ,quickCheckWithResult,maxSize,stdArgs,shrink)+import qualified Test.QuickCheck as T  import Test.QuickCheck.Test(isSuccess) @@ -77,7 +78,7 @@         && det mc == 38 :+ (-3)         && det (feye 2) == -1     where-        m = (3><3) +        m = (3><3)             [ 1, 2, 3             , 4, 5, 7             , 2, 8, 4 :: Double@@ -85,7 +86,7 @@         mc = (3><3)             [ 1, 2, 3             , 4, 5, 7-            , 2, 8, i+            , 2, 8, iC             ]  detTest2 = inv1 |~| inv2 && [det1] ~~ [det2]@@ -126,8 +127,8 @@ mbCholTest = utest "mbCholTest" (ok1 && ok2) where     m1 = (2><2) [2,5,5,8 :: Double]     m2 = (2><2) [3,5,5,9 :: Complex Double]-    ok1 = mbCholSH m1 == Nothing-    ok2 = mbCholSH m2 == Just (chol m2)+    ok1 = mbChol (trustSym m1) == Nothing+    ok2 = mbChol (trustSym m2) == Just (chol $ trustSym m2)  --------------------------------------------------------------------- @@ -136,7 +137,7 @@                 2,4,0,                -2,2,1]     m = 3 |> [1,2,3]-    c = a <> trans a+    c = a <> tr a     dat = gaussianSample 7 (10^6) m c  randomTestUniform = c :~1~: snd (meanCov dat) where@@ -170,54 +171,54 @@  normsVTest = TestList [     utest "normv2CD" $ norm2PropC v-  , utest "normv2CF" $ norm2PropC (single v)+--  , utest "normv2CF" $ norm2PropC (single v) #ifndef NONORMVTEST   , utest "normv2D"  $ norm2PropR x-  , utest "normv2F"  $ norm2PropR (single x)+--  , utest "normv2F"  $ norm2PropR (single x) #endif-  , utest "normv1CD" $ norm1 v          == 8-  , utest "normv1CF" $ norm1 (single v) == 8-  , utest "normv1D"  $ norm1 x          == 6-  , utest "normv1F"  $ norm1 (single x) == 6+  , utest "normv1CD" $ norm_1 v          == 8+--  , utest "normv1CF" $ norm_1 (single v) == 8+  , utest "normv1D"  $ norm_1 x          == 6+--  , utest "normv1F"  $ norm_1 (single x) == 6 -  , utest "normvInfCD" $ normInf v          == 5-  , utest "normvInfCF" $ normInf (single v) == 5-  , utest "normvInfD"  $ normInf x          == 3-  , utest "normvInfF"  $ normInf (single x) == 3+  , utest "normvInfCD" $ norm_Inf v          == 5+--  , utest "normvInfCF" $ norm_Inf (single v) == 5+  , utest "normvInfD"  $ norm_Inf x          == 3+--  , utest "normvInfF"  $ norm_Inf (single x) == 3   ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double)          x = fromList [1,2,-3] :: Vector Double #ifndef NONORMVTEST-         norm2PropR a = norm2 a =~= sqrt (udot a a)+         norm2PropR a = norm_2 a =~= sqrt (udot a a) #endif-         norm2PropC a = norm2 a =~= realPart (sqrt (a <.> a))+         norm2PropC a = norm_2 a =~= realPart (sqrt (a `dot` a))          a =~= b = fromList [a] |~| fromList [b]  normsMTest = TestList [-    utest "norm2mCD" $ pnorm PNorm2 v          =~= 8.86164970498005-  , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005-  , utest "norm2mD"  $ pnorm PNorm2 x          =~= 5.96667765076216-  , utest "norm2mF"  $ pnorm PNorm2 (single x) =~= 5.96667765076216+    utest "norm2mCD" $ norm_2 v          =~= 8.86164970498005+--  , utest "norm2mCF" $ norm_2 (single v) =~= 8.86164970498005+  , utest "norm2mD"  $ norm_2 x          =~= 5.96667765076216+--  , utest "norm2mF"  $ norm_2 (single x) =~= 5.96667765076216 -  , utest "norm1mCD" $ pnorm PNorm1 v          == 9-  , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9-  , utest "norm1mD"  $ pnorm PNorm1 x          == 7-  , utest "norm1mF"  $ pnorm PNorm1 (single x) == 7+  , utest "norm1mCD" $ norm_1 v          == 9+--  , utest "norm1mCF" $ norm_1 (single v) == 9+  , utest "norm1mD"  $ norm_1 x          == 7+--  , utest "norm1mF"  $ norm_1 (single x) == 7 -  , utest "normmInfCD" $ pnorm Infinity v          == 12-  , utest "normmInfCF" $ pnorm Infinity (single v) == 12-  , utest "normmInfD"  $ pnorm Infinity x          == 8-  , utest "normmInfF"  $ pnorm Infinity (single x) == 8+  , utest "normmInfCD" $ norm_Inf v          == 12+--  , utest "normmInfCF" $ norm_Inf (single v) == 12+  , utest "normmInfD"  $ norm_Inf x          == 8+--  , utest "normmInfF"  $ norm_Inf (single x) == 8 -  , utest "normmFroCD" $ pnorm Frobenius v          =~= 8.88819441731559-  , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559-  , utest "normmFroD"  $ pnorm Frobenius x          =~= 6.24499799839840-  , utest "normmFroF"  $ pnorm Frobenius (single x) =~~= 6.24499799839840+  , utest "normmFroCD" $ norm_Frob v          =~= 8.88819441731559+--  , utest "normmFroCF" $ norm_Frob (single v) =~~= 8.88819441731559+  , utest "normmFroD"  $ norm_Frob x          =~= 6.24499799839840+--  , utest "normmFroF"  $ norm_Frob (single x) =~~= 6.24499799839840 - ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double)+ ] where v = (2><2) [1,-2*iC,3:+4,7] :: Matrix (Complex Double)          x = (2><2) [1,2,-3,5] :: Matrix Double          a =~= b = fromList [a] :~10~: fromList [b]-         a =~~= b = fromList [a] :~5~: fromList [b]+--       a =~~= b = fromList [a] :~5~: fromList [b]  --------------------------------------------------------------------- @@ -232,7 +233,7 @@   , utest "prodD"  $ prodProp v   , utest "prodF"  $ prodProp (single v)  ] where v = fromList [1,2,3] :: Vector Double-         z = fromList [1,2-i,3+i]+         z = fromList [1,2-iC,3+iC]          prodProp x = prodElements x == product (toList x)  ---------------------------------------------------------------------@@ -246,7 +247,7 @@  --------------------------------------------------------------------- -conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m)+conjuTest m = cmap conjugate (flatten (conj (tr m))) == flatten (tr m)  --------------------------------------------------------------------- @@ -302,7 +303,7 @@  -- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs --successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool-successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0)+successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (size v - 1) v))) (v ! 0)    where stp e  = do                   ep <- lift_maybe $ state_get                   if t e ep@@ -311,7 +312,7 @@  -- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input --successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b-successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0)+successive f v = evalState (mapVectorM stp (subVector 1 (size v - 1) v)) (v ! 0)    where stp  e = do                   ep <- state_get                   state_put e@@ -358,7 +359,7 @@                    ,0,1,7                    ,0,0,4]          &&-         toList (flatten x) == [1,0,0,0,1,0,0,0,1] +         toList (flatten x) == [1,0,0,0,1,0,0,0,1]  -------------------------------------------------------------------------------- @@ -373,31 +374,22 @@  -------------------------------------------------------------------------------- -kroneckerTest = utest "kronecker" ok-  where-    a,x,b :: Matrix Double-    a = (3><4) [1..]-    x = (4><2) [3,5..]-    b = (2><5) [0,5..]-    v1 = vec (a <> x <> b)-    v2 = (trans b `kronecker` a) <> vec x-    s = trans b <> b-    v3 = vec s-    v4 = (dup 5 :: Matrix Double) <> vech s-    ok = v1 == v2 && v3 == v4-      && vtrans 1 a == trans a-      && vtrans (rows a) a == asColumn (vec a)+sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix))  -------------------------------------------------------------------------------- -sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix))+staticTest = utest "static" (fst $ checkT (undefined :: L 3 5))  -------------------------------------------------------------------------------- -staticTest = utest "static" (fst $ checkT (undefined :: L 3 5))+intTest = utest "int ops" (fst $ checkT (undefined :: Matrix I))  -------------------------------------------------------------------------------- +modularTest = utest "modular ops" (fst $ checkT (undefined :: Matrix (Mod 13 I)))++--------------------------------------------------------------------------------+ indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3   where     l = map g (toList (f x))@@ -410,12 +402,157 @@  -------------------------------------------------------------------------------- +sliceTest = utest "slice test" $ and+    [ testSlice (chol . trustSym)  (gen 5 :: Matrix R)+    , testSlice (chol . trustSym)  (gen 5 :: Matrix C)+    , testSlice qr    (rec :: Matrix R)+    , testSlice qr    (rec :: Matrix C)+    , testSlice hess  (agen 5 :: Matrix R)+    , testSlice hess  (agen 5 :: Matrix C)+    , testSlice schur (agen 5 :: Matrix R)+    , testSlice schur (agen 5 :: Matrix C)+    , testSlice lu    (agen 5 :: Matrix R)+    , testSlice lu    (agen 5 :: Matrix C)+    , testSlice (luSolve (luPacked (agen 5 :: Matrix R))) (agen 5)+    , testSlice (luSolve (luPacked (agen 5 :: Matrix C))) (agen 5)+    , test_lus (agen 5 :: Matrix R)+    , test_lus (agen 5 :: Matrix C)++    , testSlice eig   (agen 5 :: Matrix R)+    , testSlice eig   (agen 5 :: Matrix C)+    , testSlice (eigSH . trustSym) (gen 5 :: Matrix R)+    , testSlice (eigSH . trustSym) (gen 5 :: Matrix C)+    , testSlice eigenvalues   (agen 5 :: Matrix R)+    , testSlice eigenvalues   (agen 5 :: Matrix C)+    , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix R)+    , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix C)++    , testSlice svd           (rec :: Matrix R)+    , testSlice thinSVD       (rec :: Matrix R)+    , testSlice compactSVD     (rec :: Matrix R)+    , testSlice leftSV        (rec :: Matrix R)+    , testSlice rightSV       (rec :: Matrix R)+    , testSlice singularValues (rec :: Matrix R)++    , testSlice svd           (rec :: Matrix C)+    , testSlice thinSVD       (rec :: Matrix C)+    , testSlice compactSVD     (rec :: Matrix C)+    , testSlice leftSV        (rec :: Matrix C)+    , testSlice rightSV       (rec :: Matrix C)+    , testSlice singularValues (rec :: Matrix C)++    , testSlice (linearSolve (agen 5:: Matrix R)) (agen 5)+    , testSlice (flip linearSolve (agen 5:: Matrix R)) (agen 5)++    , testSlice (linearSolve (agen 5:: Matrix C)) (agen 5)+    , testSlice (flip linearSolve (agen 5:: Matrix C)) (agen 5)++    , testSlice (linearSolveLS (ogen 5:: Matrix R)) (ogen 5)+    , testSlice (flip linearSolveLS (ogen 5:: Matrix R)) (ogen 5)++    , testSlice (linearSolveLS (ogen 5:: Matrix C)) (ogen 5)+    , testSlice (flip linearSolveLS (ogen 5:: Matrix C)) (ogen 5)++    , testSlice (linearSolveSVD (ogen 5:: Matrix R)) (ogen 5)+    , testSlice (flip linearSolveSVD (ogen 5:: Matrix R)) (ogen 5)++    , testSlice (linearSolveSVD (ogen 5:: Matrix C)) (ogen 5)+    , testSlice (flip linearSolveSVD (ogen 5:: Matrix C)) (ogen 5)++    , testSlice (linearSolveLS (ugen 5:: Matrix R)) (ugen 5)+    , testSlice (flip linearSolveLS (ugen 5:: Matrix R)) (ugen 5)++    , testSlice (linearSolveLS (ugen 5:: Matrix C)) (ugen 5)+    , testSlice (flip linearSolveLS (ugen 5:: Matrix C)) (ugen 5)++    , testSlice (linearSolveSVD (ugen 5:: Matrix R)) (ugen 5)+    , testSlice (flip linearSolveSVD (ugen 5:: Matrix R)) (ugen 5)++    , testSlice (linearSolveSVD (ugen 5:: Matrix C)) (ugen 5)+    , testSlice (flip linearSolveSVD (ugen 5:: Matrix C)) (ugen 5)++    , testSlice ((<>) (ogen 5:: Matrix R)) (gen 5)+    , testSlice (flip (<>) (gen 5:: Matrix R)) (ogen 5)+    , testSlice ((<>) (ogen 5:: Matrix C)) (gen 5)+    , testSlice (flip (<>) (gen 5:: Matrix C)) (ogen 5)+    , testSlice ((<>) (ogen 5:: Matrix Float)) (gen 5)+    , testSlice (flip (<>) (gen 5:: Matrix Float)) (ogen 5)+    , testSlice ((<>) (ogen 5:: Matrix (Complex Float))) (gen 5)+    , testSlice (flip (<>) (gen 5:: Matrix (Complex Float))) (ogen 5)+    , testSlice ((<>) (ogen 5:: Matrix I)) (gen 5)+    , testSlice (flip (<>) (gen 5:: Matrix I)) (ogen 5)+    , testSlice ((<>) (ogen 5:: Matrix Z)) (gen 5)+    , testSlice (flip (<>) (gen 5:: Matrix Z)) (ogen 5)++    , testSlice ((<>) (ogen 5:: Matrix (I ./. 7))) (gen 5)+    , testSlice (flip (<>) (gen 5:: Matrix (I ./. 7))) (ogen 5)+    , testSlice ((<>) (ogen 5:: Matrix (Z ./. 7))) (gen 5)+    , testSlice (flip (<>) (gen 5:: Matrix (Z ./. 7))) (ogen 5)++    , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ trustSym $ gen 5)+    , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ trustSym $ gen 5)+    , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix R)) (agen 5)+    , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix C)) (agen 5)++    , ok_qrgr        (rec :: Matrix R)+    , ok_qrgr        (rec :: Matrix C)+    , testSlice (test_qrgr 4 tau1) qrr1+    , testSlice (test_qrgr 4 tau2) qrr2+    ]+  where+    QR qrr1 tau1 = qrRaw (rec :: Matrix R)+    QR qrr2 tau2 = qrRaw (rec :: Matrix C)++    test_qrgr n t x = qrgr n (QR x t)++    ok_qrgr x = simeq 1E-15 q q'+      where+        (q,_) = qr x+        atau = qrRaw x+        q' = qrgr (rows q) atau++    simeq eps a b =  not $ magnit eps (norm_1 $ flatten (a-b))++    test_lus m = testSlice f lup+      where+        f x = luSolve (LU x p) m+        (LU lup p) = luPacked m++    gen :: Numeric t => Int -> Matrix t+    gen n = diagRect 1 (konst 5 n) n n++    agen :: (Numeric t, Num (Vector t))=> Int -> Matrix t+    agen n = gen n + fromInt ((n><n)[0..])++    ogen :: (Numeric t, Num (Vector t))=> Int -> Matrix t+    ogen n = gen n === gen n++    ugen :: (Numeric t, Num (Vector t))=> Int -> Matrix t+    ugen n = takeRows 3 (gen n)+++    rec :: Numeric t => Matrix t+    rec = subMatrix (0,0) (4,5) (gen 5)++    testSlice f x@(size->sz@(r,c)) = all (==f x) (map f (g y1 ++ g y2))+      where+        subm = subMatrix+        g y = [ subm (a*r,b*c) sz y | a <-[0..2], b <- [0..2]]+        h z = fromBlocks (replicate 3 (replicate 3 z))+        y1  = h x+        y2  = (tr . h . tr) x++++--------------------------------------------------------------------------------+ -- | All tests must pass with a maximum dimension of about 20 --  (some tests may fail with bigger sizes due to precision loss). runTests :: Int  -- ^ maximum dimension          -> IO () runTests n = do-    let test p = qCheck n p+    let test :: forall t . T.Testable t => t -> IO ()+        test p = qCheck n p     putStrLn "------ index"     test( \m -> indexProp id flatten (single (m :: RM)) )     test( \v -> indexProp id id (single (v :: Vector Double)) )@@ -430,11 +567,11 @@     test (multProp1 10 . cConsist)     test (multProp2 10 . rConsist)     test (multProp2 10 . cConsist)-    putStrLn "------ mult Float"-    test (multProp1  6 . (single *** single) . rConsist)-    test (multProp1  6 . (single *** single) . cConsist)-    test (multProp2  6 . (single *** single) . rConsist)-    test (multProp2  6 . (single *** single) . cConsist)+--    putStrLn "------ mult Float"+--    test (multProp1  6 . (single *** single) . rConsist)+--    test (multProp1  6 . (single *** single) . cConsist)+--    test (multProp2  6 . (single *** single) . rConsist)+--    test (multProp2  6 . (single *** single) . cConsist)     putStrLn "------ sub-trans"     test (subProp . rM)     test (subProp . cM)@@ -450,9 +587,12 @@     putStrLn "------ luSolve"     test (linearSolveProp (luSolve.luPacked) . rSqWC)     test (linearSolveProp (luSolve.luPacked) . cSqWC)+    putStrLn "------ ldlSolve"+    test (linearSolvePropH (ldlSolve.ldlPacked) . rSymWC)+    test (linearSolvePropH (ldlSolve.ldlPacked) . cSymWC)     putStrLn "------ cholSolve"-    test (linearSolveProp (cholSolve.chol) . rPosDef)-    test (linearSolveProp (cholSolve.chol) . cPosDef)+    test (linearSolveProp (cholSolve.chol.trustSym) . rPosDef)+    test (linearSolveProp (cholSolve.chol.trustSym) . cPosDef)     putStrLn "------ luSolveLS"     test (linearSolveProp linearSolveLS . rSqWC)     test (linearSolveProp linearSolveLS . cSqWC)@@ -467,16 +607,16 @@     putStrLn "------ svd"     test (svdProp1  . rM)     test (svdProp1  . cM)-    test (svdProp1a svdR)-    test (svdProp1a svdC)-    test (svdProp1a svdRd)-    test (svdProp1b svdR)-    test (svdProp1b svdC)-    test (svdProp1b svdRd)-    test (svdProp2 thinSVDR)-    test (svdProp2 thinSVDC)-    test (svdProp2 thinSVDRd)-    test (svdProp2 thinSVDCd)+    test (svdProp1a svd . rM)+    test (svdProp1a svd . cM)+--    test (svdProp1a svdRd)+    test (svdProp1b svd . rM)+    test (svdProp1b svd . cM)+--    test (svdProp1b svdRd)+    test (svdProp2 thinSVD . rM)+    test (svdProp2 thinSVD . cM)+--    test (svdProp2 thinSVDRd)+--    test (svdProp2 thinSVDCd)     test (svdProp3  . rM)     test (svdProp3  . cM)     test (svdProp4  . rM)@@ -487,12 +627,12 @@     test (svdProp6b)     test (svdProp7  . rM)     test (svdProp7  . cM)-    putStrLn "------ svdCd"+--    putStrLn "------ svdCd" #ifdef NOZGESDD-    putStrLn "Omitted"+--    putStrLn "Omitted" #else-    test (svdProp1a svdCd)-    test (svdProp1b svdCd)+--    test (svdProp1a svdCd)+--    test (svdProp1b svdCd) #endif     putStrLn "------ eig"     test (eigSHProp . rHer)@@ -510,10 +650,10 @@     test (qrProp     . rM)     test (qrProp     . cM)     test (rqProp     . rM)-    test (rqProp     . cM)+--    test (rqProp     . cM)     test (rqProp1     . cM)     test (rqProp2     . cM)-    test (rqProp3     . cM)+--    test (rqProp3     . cM)     putStrLn "------ hess"     test (hessProp   . rSq)     test (hessProp   . cSq)@@ -523,8 +663,8 @@     putStrLn "------ chol"     test (cholProp   . rPosDef)     test (cholProp   . cPosDef)-    test (exactProp  . rPosDef)-    test (exactProp  . cPosDef)+--    test (exactProp  . rPosDef)+--    test (exactProp  . cPosDef)     putStrLn "------ expm"     test (expmDiagProp . complex. rSqWC)     test (expmDiagProp . cSqWC)@@ -534,12 +674,12 @@     test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))     test (\u -> cos u * tan u |~| sin (u::RM))     test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary-    putStrLn "------ vector operations - Float"-    test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))-    test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary-    test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))-    test (\u -> cos u * tan u |~~| sin (u::FM))-    test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary+--    putStrLn "------ vector operations - Float"+--    test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))+--    test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary+--    test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))+--    test (\u -> cos u * tan u |~~| sin (u::FM))+--    test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary     putStrLn "------ read . show"     test (\m -> (m::RM) == read (show m))     test (\m -> (m::CM) == read (show m))@@ -557,8 +697,8 @@         , utest "expm1" (expmTest1)         , utest "expm2" (expmTest2)         , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)-        , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM)-        , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0+        , utest "arith2" $ ((scalar (1+iC) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*iC-51) :: CM)+        , utest "arith3" $ exp (scalar iC * ones(10,10)*pi) + 1 |~| 0         , utest "<\\>"   $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] --        , utest "gamma" (gamma 5 == 24.0) --        , besselTest@@ -566,10 +706,10 @@         , utest "randomGaussian" randomTestGaussian         , utest "randomUniform" randomTestUniform         , utest "buildVector/Matrix" $-                        complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral-                     && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0)-        , utest "rank" $  rank ((2><3)[1,0,0,1,5*eps,0]) == 1-                       && rank ((2><3)[1,0,0,1,7*eps,0]) == 2+                        complex (10 |> [0::Double ..]) == build 10 id+                     && ident 5 == build (5,5) (\r c -> if r==c then 1::Double else 0)+        , utest "rank" $  rank ((2><3)[1,0,0,1,5*peps,0::Double]) == 1+                       && rank ((2><3)[1,0,0,1,7*peps,0::Double]) == 2         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)         , mbCholTest         , utest "offset" offsetTest@@ -583,21 +723,23 @@         , conformTest         , accumTest         , convolutionTest-        , kroneckerTest         , sparseTest         , staticTest+        , intTest+        , modularTest+        , sliceTest         ]     when (errors c + failures c > 0) exitFailure     return ()   -- single precision approximate equality-infixl 4 |~~|-a |~~| b = a :~6~: b+-- infixl 4 |~~|+-- a |~~| b = a :~6~: b  makeUnitary v | realPart n > 1    = v / scalar n               | otherwise = v-    where n = sqrt (v <.> v)+    where n = sqrt (v `dot` v)  -- -- | Some additional tests on big matrices. They take a few minutes. -- runBigTests :: IO ()@@ -620,6 +762,8 @@     mkVecBench     multBench     cholBench+    luBench+    luBench_2     svdBench     eigBench     putStrLn ""@@ -663,9 +807,9 @@   manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs-manyvec3 xs = sum $ map (pnorm PNorm2 . (\x -> fromList [x,x**2,x**3])) xs+manyvec3 xs = sum $ map (norm_2 . (\x -> fromList [x,x**2,x**3])) xs -manyvec4 xs = sum $ map (pnorm PNorm2 . (\x -> vec3 x (x**2) (x**3))) xs+manyvec4 xs = sum $ map (norm_2 . (\x -> vec3 x (x**2) (x**3))) xs  vec3 :: Double -> Double -> Double -> Vector Double vec3 a b c = runSTVector $ do@@ -690,11 +834,11 @@  subBench = do     putStrLn ""-    let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v))-    time "0.1M subVector   " (g (konst 1 (1+10^5) :: Vector Double) @> 0)+    let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (size v -1) v))+    time "0.1M subVector   " (g (konst 1 (1+10^5) :: Vector Double) ! 0)     let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))-    time "subVector-join  3" (f (ident  3 :: Matrix Double) @@>(0,0))-    time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0))+    time "subVector-join  3" (f (ident  3 :: Matrix Double) `atIndex` (0,0))+    time "subVector-join 10" (f (ident 10 :: Matrix Double) `atIndex` (0,0))  -------------------------------- @@ -719,10 +863,10 @@  eigBench = do     let m = reshape 1000 (randomVector 777 Uniform (1000*1000))-        s = m + trans m+        s = m + tr m     m `seq` s `seq` putStrLn ""-    time "eigenvalues  symmetric 1000x1000" (eigenvaluesSH' m)-    time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m)+    time "eigenvalues  symmetric 1000x1000" (eigenvaluesSH (trustSym m))+    time "eigenvectors symmetric 1000x1000" (snd $ eigSH (trustSym m))     time "eigenvalues  general   1000x1000" (eigenvalues m)     time "eigenvectors general   1000x1000" (snd $ eig m) @@ -731,7 +875,7 @@ svdBench = do     let a = reshape 500  (randomVector 777 Uniform (3000*500))         b = reshape 1000 (randomVector 777 Uniform (1000*1000))-        fv (_,_,v) = v@@>(0,0)+        fv (_,_,v) = v `atIndex` (0,0)     a `seq` b `seq` putStrLn ""     time "singular values  3000x500" (singularValues a)     time "thin svd         3000x500" (fv $ thinSVD a)@@ -743,26 +887,28 @@  solveBenchN n = do     let x = uniformSample 777 (2*n) (replicate n (-1,1))-        a = trans x <> x+        a = tr x <> x         b = asColumn $ randomVector 666 Uniform n     a `seq` b `seq` putStrLn ""     time ("svd solve " ++ show n) (linearSolveSVD a b)     time (" ls solve " ++ show n) (linearSolveLS a b)     time ("    solve " ++ show n) (linearSolve a b)-    time ("cholSolve " ++ show n) (cholSolve (chol a) b)+--    time (" LU solve " ++ show n) (luSolve (luPacked a) b)+    time ("LDL solve " ++ show n) (ldlSolve (ldlPacked (trustSym a)) b)+    time ("cholSolve " ++ show n) (cholSolve (chol $ trustSym a) b)  solveBench = do     solveBenchN 500     solveBenchN 1000-    -- solveBenchN 1500+    solveBenchN 1500  --------------------------------  cholBenchN n = do     let x = uniformSample 777 (2*n) (replicate n (-1,1))-        a = trans x <> x+        a = tr x <> x     a `seq` putStr ""-    time ("chol " ++ show n) (chol a)+    time ("chol " ++ show n) (chol $ trustSym a)  cholBench = do     putStrLn ""@@ -771,3 +917,32 @@     cholBenchN 300 --    cholBenchN 150 --    cholBenchN 50++--------------------------------------------------------------------------------++luBenchN f n x msg = do+    let m = diagRect 1 (fromList (replicate n x)) n n+    m `seq` putStr ""+    time (msg ++ " "++ show n) (rnf $ f m)++luBench = do+    putStrLn ""+    luBenchN luPacked  1000 (5::R)          "luPacked  Double    "+    luBenchN luPacked' 1000 (5::R)          "luPacked' Double    "+    luBenchN luPacked' 1000 (5::Mod 9973 I) "luPacked' I mod 9973"+    luBenchN luPacked' 1000 (5::Mod 9973 Z) "luPacked' Z mod 9973"++luBenchN_2 f g n x msg = do+    let m = diagRect 1 (fromList (replicate n x)) n n+        b = flipud m+    m `seq` b `seq` putStr ""+    time (msg ++ " "++ show n) (f (g m) b)++luBench_2 = do+    putStrLn ""+    luBenchN_2 luSolve  luPacked  500 (5::R)          "luSolve .luPacked  Double    "+    luBenchN_2 luSolve' luPacked' 500 (5::R)          "luSolve'.luPacked' Double    "+    luBenchN_2 luSolve' luPacked' 500 (5::Mod 9973 I) "luSolve'.luPacked' I mod 9973"+    luBenchN_2 luSolve' luPacked' 500 (5::Mod 9973 Z) "luSolve'.luPacked' Z mod 9973"++

@@ -1,5 +1,4 @@-{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-}-{-# OPTIONS_GHC -fno-warn-unused-imports #-}+{-# LANGUAGE FlexibleContexts, UndecidableInstances, FlexibleInstances #-} ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra.Tests.Instances@@ -15,9 +14,9 @@ module Numeric.LinearAlgebra.Tests.Instances(     Sq(..),     rSq,cSq,     Rot(..),    rRot,cRot,-    Her(..),    rHer,cHer,+                rHer,cHer,     WC(..),     rWC,cWC,-    SqWC(..),   rSqWC, cSqWC,+    SqWC(..),   rSqWC, cSqWC, rSymWC, cSymWC,     PosDef(..), rPosDef, cPosDef,     Consistent(..), rConsist, cConsist,     RM,CM, rM,cM,@@ -26,15 +25,11 @@  import System.Random -import Numeric.LinearAlgebra-import Numeric.LinearAlgebra.Devel-import Numeric.Container+import Numeric.LinearAlgebra.HMatrix hiding (vector) import Control.Monad(replicateM)-import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector-                      ,sized,classify,Testable,Property-                      ,quickCheckWith,maxSize,stdArgs,shrink)+import Test.QuickCheck(Arbitrary,arbitrary,choose,vector,sized,shrink) -#if MIN_VERSION_QuickCheck(2,0,0)+ shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] shrinkListElementwise []     = [] shrinkListElementwise (x:xs) = [ y:xs | y  <- shrink x                 ]@@ -42,41 +37,16 @@  shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)] shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ]-#endif -#if MIN_VERSION_QuickCheck(2,1,1)-#else-instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where-    arbitrary = do-        re <- arbitrary-        im <- arbitrary-        return (re :+ im)--#if MIN_VERSION_QuickCheck(2,0,0)-    shrink (re :+ im) = -        [ u :+ v | (u,v) <- shrinkPair (re,im) ]-#else-    -- this has been moved to the 'Coarbitrary' class in QuickCheck 2-    coarbitrary = undefined -#endif--#endif- chooseDim = sized $ \m -> choose (1,max 1 m)  instance (Field a, Arbitrary a) => Arbitrary (Vector a) where      arbitrary = do m <- chooseDim                    l <- vector m                    return $ fromList l--#if MIN_VERSION_QuickCheck(2,0,0)     -- shrink any one of the components     shrink = map fromList . shrinkListElementwise . toList -#else-    coarbitrary = undefined-#endif- instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where      arbitrary = do         m <- chooseDim@@ -84,17 +54,12 @@         l <- vector (m*n)         return $ (m><n) l -#if MIN_VERSION_QuickCheck(2,0,0)     -- shrink any one of the components     shrink a = map (rows a >< cols a)                . shrinkListElementwise                . concat . toLists                       $ a-#else-    coarbitrary = undefined-#endif - -- a square matrix newtype (Sq a) = Sq (Matrix a) deriving Show instance (Element a, Arbitrary a) => Arbitrary (Sq a) where@@ -103,11 +68,7 @@         l <- vector (n*n)         return $ Sq $ (n><n) l -#if MIN_VERSION_QuickCheck(2,0,0)     shrink (Sq a) = [ Sq b | b <- shrink a ]-#else-    coarbitrary = undefined-#endif   -- a unitary matrix@@ -118,24 +79,14 @@         let (q,_) = qr m         return (Rot q) -#if MIN_VERSION_QuickCheck(2,0,0)-#else-    coarbitrary = undefined-#endif - -- a complex hermitian or real symmetric matrix-newtype (Her a) = Her (Matrix a) deriving Show-instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where+instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Herm a) where     arbitrary = do         Sq m <- arbitrary         let m' = m/2-        return $ Her (m' + ctrans m')+        return $ sym m' -#if MIN_VERSION_QuickCheck(2,0,0)-#else-    coarbitrary = undefined-#endif  class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a instance ArbitraryField Double@@ -144,7 +95,7 @@  -- a well-conditioned general matrix (the singular values are between 1 and 100) newtype (WC a) = WC (Matrix a) deriving Show-instance (ArbitraryField a) => Arbitrary (WC a) where+instance (Numeric a, ArbitraryField a) => Arbitrary (WC a) where     arbitrary = do         m <- arbitrary         let (u,_,v) = svd m@@ -153,48 +104,33 @@             n = min r c         sv' <- replicateM n (choose (1,100))         let s = diagRect 0 (fromList sv') r c-        return $ WC (u `mXm` real s `mXm` trans v)--#if MIN_VERSION_QuickCheck(2,0,0)-#else-    coarbitrary = undefined-#endif+        return $ WC (u <> real s <> tr v)   -- a well-conditioned square matrix (the singular values are between 1 and 100) newtype (SqWC a) = SqWC (Matrix a) deriving Show-instance (ArbitraryField a) => Arbitrary (SqWC a) where+instance (ArbitraryField a, Numeric a) => Arbitrary (SqWC a) where     arbitrary = do         Sq m <- arbitrary         let (u,_,v) = svd m             n = rows m         sv' <- replicateM n (choose (1,100))         let s = diag (fromList sv')-        return $ SqWC (u `mXm` real s `mXm` trans v)--#if MIN_VERSION_QuickCheck(2,0,0)-#else-    coarbitrary = undefined-#endif+        return $ SqWC (u <> real s <> tr v)   -- a positive definite square matrix (the eigenvalues are between 0 and 100) newtype (PosDef a) = PosDef (Matrix a) deriving Show-instance (ArbitraryField a, Num (Vector a)) +instance (Numeric a, ArbitraryField a, Num (Vector a))      => Arbitrary (PosDef a) where     arbitrary = do-        Her m <- arbitrary+        m <- arbitrary         let (_,v) = eigSH m-            n = rows m+            n = rows (unSym m)         l <- replicateM n (choose (0,100))         let s = diag (fromList l)-            p = v `mXm` real s `mXm` ctrans v-        return $ PosDef (0.5 * p + 0.5 * ctrans p)--#if MIN_VERSION_QuickCheck(2,0,0)-#else-    coarbitrary = undefined-#endif+            p = v <> real s <> tr v+        return $ PosDef (0.5 * p + 0.5 * tr p)   -- a pair of matrices that can be multiplied@@ -208,11 +144,7 @@         lb <- vector (k*m)         return $ Consistent ((n><k) la, (k><m) lb) -#if MIN_VERSION_QuickCheck(2,0,0)     shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ]-#else-    coarbitrary = undefined-#endif   @@ -228,8 +160,8 @@ zM m = m :: ZM  -rHer (Her m) = m :: RM-cHer (Her m) = m :: CM+rHer m = unSym m :: RM+cHer m = unSym m :: CM  rRot (Rot m) = m :: RM cRot (Rot m) = m :: CM@@ -242,6 +174,9 @@  rSqWC (SqWC m) = m :: RM cSqWC (SqWC m) = m :: CM++rSymWC (SqWC m) = sym m :: Herm R+cSymWC (SqWC m) = sym m :: Herm C  rPosDef (PosDef m) = m :: RM cPosDef (PosDef m) = m :: CM

@@ -1,5 +1,6 @@-{-# LANGUAGE CPP, FlexibleContexts #-}-{-# OPTIONS_GHC -fno-warn-unused-imports #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE TypeFamilies #-}+ ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra.Tests.Properties@@ -13,7 +14,7 @@ -}  module Numeric.LinearAlgebra.Tests.Properties (-    dist, (|~|), (~~), (~:), Aprox((:~)),+    dist, (|~|), (~~), (~:), Aprox((:~)), (~=),     zeros, ones,     square,     unitary,@@ -27,7 +28,7 @@     pinvProp,     detProp,     nullspaceProp,-    bugProp,+--    bugProp,     svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,     svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,     eigProp, eigSHProp, eigProp2, eigSHProp2,@@ -38,23 +39,21 @@     expmDiagProp,     multProp1, multProp2,     subProp,-    linearSolveProp, linearSolveProp2+    linearSolveProp, linearSolvePropH, linearSolveProp2 ) where -import Numeric.Container-import Numeric.LinearAlgebra --hiding (real,complex)-import Numeric.LinearAlgebra.LAPACK-import Debug.Trace-import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector-                      ,sized,classify,Testable,Property-                      ,quickCheckWith,maxSize,stdArgs,shrink)+import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary)+import Test.QuickCheck +(~=) :: Double -> Double -> Bool+a ~= b = abs (a - b) < 1e-10+ trivial :: Testable a => Bool -> a -> Property trivial = (`classify` "trivial")  -- relative error-dist :: (Normed c t, Num (c t)) => c t -> c t -> Double-dist = relativeError Infinity+dist :: (Num a, Normed a) => a -> a -> Double+dist = relativeError norm_Inf  infixl 4 |~| a |~| b = a :~10~: b@@ -71,11 +70,11 @@ square m = rows m == cols m  -- orthonormal columns-orthonormal m = ctrans m <> m |~| ident (cols m)+orthonormal m = tr m <> m |~| ident (cols m)  unitary m = square m && orthonormal m -hermitian m = square m && m |~| ctrans m+hermitian m = square m && m |~| tr m  wellCond m = rcond m > 1/100 @@ -83,12 +82,12 @@     where (e,_v) = eigSH m  upperTriang m = rows m == 1 || down == z-    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))-          z = konst 0 (dim down)+    where down = fromList $ concat $ zipWith drop [1..] (toLists (tr m))+          z = konst 0 (size down)  upperHessenberg m = rows m < 3 || down == z-    where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))-          z = konst 0 (dim down)+    where down = fromList $ concat $ zipWith drop [2..] (toLists (tr m))+          z = konst 0 (size down)  zeros (r,c) = reshape c (konst 0 (r*c)) @@ -116,81 +115,94 @@           s x = fromList [x]  nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)-                                     && orthonormal (fromColumns nl))-    where nl = nullspacePrec 1 m-          n = fromColumns nl+                                     && orthonormal n)+    where n = nullspaceSVD (Left (1*peps)) m (rightSV m)+          nl = toColumns n           r = rows m           c = cols m - rank m  -------------------------------------------------------------------+{- -- testcase for nonempty fpu stack -- uncommenting unitary' signature eliminates the problem-bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v-    where (u,d,v) = fullSVD m+bugProp m = m |~| u <> real d <> tr v && unitary' u && unitary' v+    where (u,d,v) = svd m           -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool           unitary' a = unitary a-+-} ------------------------------------------------------------------  -- fullSVD-svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v-    where (u,d,v) = fullSVD m+svdProp1 m = m |~| u <> real d <> tr v && unitary u && unitary v+  where+    (u,s,v) = svd m+    d = diagRect 0 s (rows m) (cols m) -svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where+svdProp1a svdfun m = m |~| u <> real d <> tr v && unitary u && unitary v+  where     (u,s,v) = svdfun m     d = diagRect 0 s (rows m) (cols m) -svdProp1b svdfun m = unitary u && unitary v where+svdProp1b svdfun m = unitary u && unitary v+  where     (u,_,v) = svdfun m  -- thinSVD-svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m)-    where (u,s,v) = thinSVDfun m+svdProp2 thinSVDfun m+    =  m |~| u <> diag (real s) <> tr v+    && orthonormal u && orthonormal v+    && size s == min (rows m) (cols m)+  where+    (u,s,v) = thinSVDfun m  -- compactSVD-svdProp3 m = (m |~| u <> real (diag s) <> trans v+svdProp3 m = (m |~| u <> real (diag s) <> tr v              && orthonormal u && orthonormal v)-    where (u,s,v) = compactSVD m+  where+    (u,s,v) = compactSVD m -svdProp4 m' = m |~| u <> real (diag s) <> trans v+svdProp4 m' = m |~| u <> real (diag s) <> tr v            && orthonormal u && orthonormal v-           && (dim s == r || r == 0 && dim s == 1)-    where (u,s,v) = compactSVD m-          m = fromBlocks [[m'],[m']]-          r = rank m'+           && (size s == r || r == 0 && size s == 1)+  where+    (u,s,v) = compactSVD m+    m = fromBlocks [[m'],[m']]+    r = rank m' -svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where-    s1       = svR  m-    s2       = svRd m-    (_,s3,_) = svdR m-    (_,s4,_) = svdRd m-    (_,s5,_) = thinSVDR m-    (_,s6,_) = thinSVDRd m+svdProp5a m = all (s1|~|) [s3,s5] where+    s1       = singularValues (m :: Matrix Double)+--  s2       = svRd m+    (_,s3,_) = svd m+--  (_,s4,_) = svdRd m+    (_,s5,_) = thinSVD m+--  (_,s6,_) = thinSVDRd m -svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where-    s1       = svC  m-    s2       = svCd m-    (_,s3,_) = svdC m-    (_,s4,_) = svdCd m-    (_,s5,_) = thinSVDC m-    (_,s6,_) = thinSVDCd m+svdProp5b m = all (s1|~|) [s3,s5] where+    s1       = singularValues (m :: Matrix (Complex Double))+--  s2       = svCd m+    (_,s3,_) = svd m+--  (_,s4,_) = svdCd m+    (_,s5,_) = thinSVD m+--  (_,s6,_) = thinSVDCd m  svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'-    where (u,s,v) = svdR m-          (s',v') = rightSVR m-          (u',s'') = leftSVR m+  where+    (u,s,v) = svd (m :: Matrix Double)+    (s',v') = rightSV m+    (u',s'') = leftSV m  svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'-    where (u,s,v) = svdC m-          (s',v') = rightSVC m-          (u',s'') = leftSVC m+  where+    (u,s,v) = svd (m :: Matrix (Complex Double))+    (s',v') = rightSV m+    (u',s'') = leftSV m  svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''-    where (u,s,v) = svd m-          (s',v') = rightSV m-          (u',_s'') = leftSV m-          s''' = singularValues m+  where+    (u,s,v) = svd m+    (s',v') = rightSV m+    (u',_s'') = leftSV m+    s''' = singularValues m  ------------------------------------------------------------------ @@ -199,12 +211,12 @@  eigSHProp m = m <> v |~| v <> real (diag s)               && unitary v-              && m |~| v <> real (diag s) <> ctrans v-    where (s, v) = eigSH m+              && m |~| v <> real (diag s) <> tr v+    where (s, v) = eigSH' m  eigProp2 m = fst (eig m) |~| eigenvalues m -eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m+eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m  ------------------------------------------------------------------ @@ -224,22 +236,22 @@     where (r,_q) = rq m  upperTriang' r = upptr (rows r) (cols r) * r |~| r-    where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1-              where t = f-c+    where upptr f c = build (f,c) $ \r' c' -> if r'-t > c' then 0 else 1+              where t = fromIntegral (f-c) -hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h+hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h     where (p,h) = hess m -schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s+schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s     where (u,s) = schur m -schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme+schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme     where (u,s) = schur m -cholProp m = m |~| ctrans c <> c && upperTriang c-    where c = chol m+cholProp m = m |~| tr c <> c && upperTriang c+    where c = chol (trustSym m) -exactProp m = chol m == chol (m+0)+exactProp m = chol (trustSym m) == chol (trustSym (m+0))  expmDiagProp m = expm (logm m) :~ 7 ~: complex m     where logm = matFunc log@@ -250,14 +262,16 @@  multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) -multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a)+multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)  linearSolveProp f m = f m m |~| ident (rows m) +linearSolvePropH f m = f m (unSym m) |~| ident (rows (unSym m))+ linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)     where q = min (rows a) (cols a)           b = a <> x           wc = rank a == q -subProp m = m == (trans . fromColumns . toRows) m+subProp m = m == (conj . tr . fromColumns . toRows) m