diff --git a/hmatrix-tests.cabal b/hmatrix-tests.cabal --- a/hmatrix-tests.cabal +++ b/hmatrix-tests.cabal @@ -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 diff --git a/src/Numeric/GSL/Tests.hs b/src/Numeric/GSL/Tests.hs --- a/src/Numeric/GSL/Tests.hs +++ b/src/Numeric/GSL/Tests.hs @@ -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]) diff --git a/src/Numeric/LinearAlgebra/Tests.hs b/src/Numeric/LinearAlgebra/Tests.hs --- a/src/Numeric/LinearAlgebra/Tests.hs +++ b/src/Numeric/LinearAlgebra/Tests.hs @@ -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> 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" + + diff --git a/src/Numeric/LinearAlgebra/Tests/Instances.hs b/src/Numeric/LinearAlgebra/Tests/Instances.hs --- a/src/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/src/Numeric/LinearAlgebra/Tests/Instances.hs @@ -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>< 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> 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> 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