{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE GADTs #-} module Test.Hermitian (testsVar, genHermitian) where import qualified Test.Mosaic as Mosaic import qualified Test.Divide as Divide import qualified Test.Generic as Generic import qualified Test.Indexed as Indexed import qualified Test.Generator as Gen import qualified Test.Logic as Logic import qualified Test.Utility as Util import Test.Mosaic (repack) import Test.Generator ((<-*#>), (<#*|>), (<.*#>), (<#*#>), (<#\#>), (<#=#>)) import Test.Utility (approxReal, approxArray, approxArrayTol, approxMatrix, equalArray, Tagged, genOrder, (!===)) import qualified Numeric.LAPACK.Orthogonal.Householder as HH import qualified Numeric.LAPACK.Matrix.HermitianPositiveDefinite as HermitianPD import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian import qualified Numeric.LAPACK.Matrix.Triangular as Triangular import qualified Numeric.LAPACK.Matrix.Square as Square import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape import qualified Numeric.LAPACK.Matrix.Layout as Layout import qualified Numeric.LAPACK.Matrix as Matrix import qualified Numeric.LAPACK.Vector as Vector import Numeric.LAPACK.Matrix.Square (Square) import Numeric.LAPACK.Matrix.Layout (Order) import Numeric.LAPACK.Matrix (General, ShapeInt, (#+#), (|||)) import Numeric.LAPACK.Vector (Vector, (.*|)) import Numeric.LAPACK.Scalar (RealOf, selectReal) import qualified Numeric.Netlib.Class as Class import qualified Type.Data.Bool as TBool import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Shape ((::+)) import Control.Applicative (liftA2, (<$>)) import qualified Data.NonEmpty.Class as NonEmptyC import qualified Data.NonEmpty as NonEmpty import qualified Data.List as List import Data.Semigroup ((<>)) import Data.Tuple.HT (uncurry3, mapFst) import qualified Test.QuickCheck as QC type FlexHermitianP pack neg zero pos sh = ArrMatrix.FullQuadratic pack (Omni.Hermitian neg zero pos) sh type HermitianP pack sh = ArrMatrix.FullQuadratic pack Omni.HermitianUnknownDefiniteness sh type HermitianPosDefP pack sh = ArrMatrix.FullQuadratic pack Omni.HermitianPositiveDefinite sh type HermitianNegDefP pack sh = ArrMatrix.FullQuadratic pack Omni.HermitianNegativeDefinite sh genHermitian :: (Logic.Dim sh, Shape.Indexed sh, Shape.Index sh ~ ix, Eq ix, Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Gen.Square sh a (HermitianP pack sh a) genHermitian p = repack p <$> Gen.hermitian generalFromHermitian :: (Layout.Packing pack) => (TBool.C neg, TBool.C zero, TBool.C pos, Shape.C sh, Class.Floating a) => FlexHermitianP pack neg zero pos sh a -> General sh sh a generalFromHermitian = Matrix.fromFull . Hermitian.toSquare stack :: (Layout.Packing pack, Class.Floating a) => (HermitianP pack ShapeInt a, General ShapeInt ShapeInt a, HermitianP pack ShapeInt a) -> Bool stack (a,b,c) = let abc = generalFromHermitian $ Hermitian.stack a b c in equalArray abc $ (Matrix.fromFull (Hermitian.toSquare a) ||| b !=== Matrix.adjoint b ||| Matrix.fromFull (Hermitian.toSquare c)) split :: (Layout.Packing pack, Class.Floating a) => HermitianP pack (ShapeInt::+ShapeInt) a -> Bool split abc = equalArray abc $ uncurry3 Hermitian.stack $ Hermitian.split abc gramian :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> General ShapeInt ShapeInt a -> Bool gramian pack x = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.gramian x) (Matrix.adjoint x <> x) gramianAdjoint :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> General ShapeInt ShapeInt a -> Bool gramianAdjoint pack x = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.gramianAdjoint x) (Matrix.adaptOrder x $ x <> Matrix.adjoint x) gramianNonAdjoint :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> General ShapeInt ShapeInt a -> Bool gramianNonAdjoint pack x = approxArray (Matrix.forceOrder (ArrMatrix.order x) $ Hermitian.gramian $ Matrix.adjoint x) (ArrMatrix.requirePacking pack $ Hermitian.gramianAdjoint x) congruenceDiagonal :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> (Vector ShapeInt ar, General ShapeInt ShapeInt a) -> Bool congruenceDiagonal pack (d,a) = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.congruenceDiagonal d a) (Matrix.adjoint a <> Matrix.scaleRowsReal d a) congruenceDiagonalAdjoint :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> (General ShapeInt ShapeInt a, Vector ShapeInt ar) -> Bool congruenceDiagonalAdjoint pack (a,d) = approxMatrix 1e-5 (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.congruenceDiagonalAdjoint a d) (Matrix.scaleColumnsReal d a <> Matrix.adjoint a) congruenceDiagonalGramian :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> General ShapeInt ShapeInt a -> Bool congruenceDiagonalGramian pack a = approxArray (Hermitian.congruenceDiagonal (Vector.one $ Matrix.height a) a) (Hermitian.relaxIndefinite $ ArrMatrix.requirePacking pack $ Hermitian.gramian a) congruence :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => (HermitianP pack ShapeInt a, General ShapeInt ShapeInt a) -> Bool congruence (b,a) = approxArray (Hermitian.toSquare $ Hermitian.congruence b a) (Square.congruence (Hermitian.toSquare b) a) congruenceAdjoint :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => (General ShapeInt ShapeInt a, HermitianP pack ShapeInt a) -> Bool congruenceAdjoint (a,b) = approxMatrix 1e-5 (Hermitian.toSquare $ Hermitian.congruenceAdjoint a b) (Square.congruenceAdjoint a $ Hermitian.toSquare b) congruenceCongruenceDiagonal :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Order -> (Vector ShapeInt ar, General ShapeInt ShapeInt a) -> Bool congruenceCongruenceDiagonal pack order (d,a) = approxArray (Hermitian.congruenceDiagonal d a) (Hermitian.congruence (repack pack $ Hermitian.diagonal order d) a) anticommutator :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool anticommutator pack (a,b) = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.anticommutator a b) ((Matrix.adjoint b <> a) #+# (Matrix.adjoint a <> b)) anticommutatorCommutative :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool anticommutatorCommutative pack (a,b) = approxMatrix 1e-5 (ArrMatrix.requirePacking pack $ Hermitian.anticommutator a b) (Hermitian.anticommutator b a) anticommutatorAdjoint :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool anticommutatorAdjoint pack (a,b) = approxArray (Matrix.forceOrder (ArrMatrix.order b) $ Hermitian.anticommutator (Matrix.adjoint a) (Matrix.adjoint b)) (ArrMatrix.requirePacking pack $ Hermitian.anticommutatorAdjoint a b) outer :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Order -> Vector ShapeInt a -> Bool outer pack order x = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.outer order x) (Matrix.outer order x x) genScaledVectors :: (NonEmptyC.Gen f, Class.Floating a, RealOf a ~ ar, Class.Real ar) => Gen.VectorInt a (ShapeInt, f (ar, Vector ShapeInt a)) genScaledVectors = Gen.listOfVector ((,) <$> Gen.scalarReal <.*#> Gen.vector) sumRank1 :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Order -> (ShapeInt, [(ar, Vector ShapeInt a)]) -> Bool sumRank1 pack order (sh,xs) = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.sumRank1 order sh xs) (Util.addMatrices (MatrixShape.general order sh sh) $ fmap (rank1 order) xs) sumRank1NonEmpty :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Order -> NonEmpty.T [] (ar, Vector ShapeInt a) -> Bool sumRank1NonEmpty pack order xs = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.sumRank1NonEmpty order xs) (NonEmpty.foldl1 (ArrMatrix.lift2 Vector.add) $ fmap (rank1 order) xs) rank1 :: (Eq size, Shape.C size, Class.Floating a) => Order -> (RealOf a, Vector size a) -> Matrix.General size size a rank1 order (r,x) = Matrix.scaleReal r $ Matrix.outer order x x genScaledVectorPairs :: (NonEmptyC.Gen f, Class.Floating a) => Gen.VectorInt a (ShapeInt, f (a, (Vector ShapeInt a, Vector ShapeInt a))) genScaledVectorPairs = flip Gen.mapGen Gen.vectorDim $ \maxElem size -> fmap ((,) size) $ NonEmptyC.genOf $ liftA2 (,) (Util.genElement maxElem) $ liftA2 (,) (Util.genVector maxElem size) (Util.genVector maxElem size) sumRank2 :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Order -> (ShapeInt, [(a, (Vector ShapeInt a, Vector ShapeInt a))]) -> Bool sumRank2 pack order (sh,xys) = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.sumRank2 order sh xys) (Util.addMatrices (MatrixShape.general order sh sh) $ fmap (rank2 order) xys) sumRank2NonEmpty :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Order -> NonEmpty.T [] (a, (Vector ShapeInt a, Vector ShapeInt a)) -> Bool sumRank2NonEmpty pack order xys = approxArray (generalFromHermitian $ ArrMatrix.requirePacking pack $ Hermitian.sumRank2NonEmpty order xys) (NonEmpty.foldl1 (ArrMatrix.lift2 Vector.add) $ fmap (rank2 order) xys) rank2 :: (Eq size, Shape.C size, Class.Floating a) => Order -> (a, (Vector size a, Vector size a)) -> Matrix.General size size a rank2 order (a,(x,y)) = let ax = a.*|x in ArrMatrix.lift2 Vector.add (Matrix.outer order ax y) (Matrix.outer order y ax) addAdjoint :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Square ShapeInt a -> Bool addAdjoint pack x = approxArray (Hermitian.toSquare $ ArrMatrix.requirePacking pack $ Hermitian.addAdjoint x) (Matrix.adjoint x #+# x) choleskyQR :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Matrix.Tall ShapeInt ShapeInt a -> QC.Property choleskyQR pack a = let qr = HH.fromMatrix a r = HH.tallExtractR qr in HH.determinantAbsolute qr > 0.1 QC.==> approxArrayTol 1e-1 (Matrix.scaleRows (Array.map signum $ Triangular.takeDiagonal r) $ Triangular.toSquare r) (Triangular.toSquare $ HermitianPD.decompose $ ArrMatrix.requirePacking pack $ gramianPosDef a) gramianPosDef :: (Layout.Packing pack, Class.Floating a) => Matrix.Tall ShapeInt ShapeInt a -> HermitianPosDefP pack ShapeInt a gramianPosDef = HermitianPD.assureFullRank . Hermitian.gramian . Matrix.fromFull genInvertible :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Gen.MatrixInt a (HermitianP pack ShapeInt a) genInvertible pack = repack pack <$> Gen.condition Util.invertible Gen.hermitian genPositiveDefinite :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Gen.MatrixInt a (HermitianPosDefP pack ShapeInt a) genPositiveDefinite pack = repack pack . gramianPosDef <$> Gen.gramian Gen.fullRankTall genNegativeDefinite :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> Gen.MatrixInt a (HermitianNegDefP pack ShapeInt a) genNegativeDefinite pack = Hermitian.negate <$> genPositiveDefinite pack determinantPD :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => HermitianPosDefP pack ShapeInt a -> Bool determinantPD a = let d = HermitianPD.determinant a in approxReal (max (selectReal 1 1e-2) (selectReal 1e-3 1e-5 * abs d)) d (Hermitian.determinant $ Hermitian.relaxIndefinite a) inversePD :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => HermitianPosDefP pack ShapeInt a -> Bool inversePD a = Divide.approxRelArray (Hermitian.inverse $ Hermitian.relaxIndefinite a) (Hermitian.relaxIndefinite $ HermitianPD.inverse a) solvePD :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => (HermitianPosDefP pack ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool solvePD (a,b) = Divide.approxRelMatrix (Hermitian.solve (Hermitian.relaxIndefinite a) b) (HermitianPD.solve a b) solveDecomposedPD :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => (HermitianPosDefP pack ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool solveDecomposedPD (a,b) = approxArrayTol (selectReal 1e-1 1e-6) (HermitianPD.solve a b) (HermitianPD.solveDecomposed (HermitianPD.decompose a) b) eigenvaluesDeterminant :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => HermitianP pack ShapeInt a -> Bool eigenvaluesDeterminant a = approxReal (selectReal 1e-1 1e-5) (Hermitian.determinant a) (Vector.product $ Hermitian.eigenvalues a) eigensystem :: (Layout.Packing pack) => (Class.Floating a, RealOf a ~ ar, Class.Real ar) => HermitianP pack ShapeInt a -> Bool eigensystem a = approxMatrix 1e-4 a $ uncurry Hermitian.congruenceDiagonalAdjoint $ mapFst Matrix.fromFull $ Hermitian.eigensystem a checkForAll :: (Show a, QC.Testable test) => Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property checkForAll gen = Util.checkForAll (Gen.run gen 3 5) checkForAllExtra :: (Show a, Show b, QC.Testable test) => QC.Gen a -> Gen.T dim tag b -> (a -> b -> test) -> Tagged tag QC.Property checkForAllExtra = Gen.withExtra checkForAll testsVar :: (Show a, Show ar, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) => [(String, Tagged a QC.Property)] testsVar = concat $ List.transpose [Util.suffix "Packed" (testsVarPacking Layout.Packed), Util.suffix "Unpacked" (testsVarPacking Layout.Unpacked)] testsVarPacking :: (Layout.Packing pack) => (Show a, Show ar, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) => Layout.PackingSingleton pack -> [(String, Tagged a QC.Property)] testsVarPacking p = ("index", checkForAll (Indexed.genMatrixIndex $ genHermitian p) Indexed.unitDot) : ("forceOrder", checkForAllExtra genOrder ((,) <$> genHermitian p <#*|> Gen.vector) Generic.forceOrder) : ("forceOrderInverse", checkForAll (genHermitian p) Generic.forceOrderInverse) : Generic.testsDistributive (Gen.asMatrixInt $ genHermitian p) ++ ("stack", checkForAll (Gen.stack3 (genHermitian p) Gen.matrix (genHermitian p)) stack) : ("split", checkForAll (genHermitian p) split) : ("gramian", checkForAll Gen.matrix (gramian p)) : ("gramianAdjoint", checkForAll Gen.matrix (gramianAdjoint p)) : ("gramianNonAdjoint", checkForAll Gen.matrix (gramianNonAdjoint p)) : ("congruenceDiagonal", checkForAll ((,) <$> Gen.vectorReal <-*#> Gen.matrix) (congruenceDiagonal p)) : ("congruence", checkForAll ((,) <$> genHermitian p <#*#> Gen.matrix) congruence) : ("congruenceDiagonalAdjoint", checkForAll ((,) <$> Gen.matrix <#*|> Gen.vectorReal) (congruenceDiagonalAdjoint p)) : ("congruenceDiagonalGramian", checkForAll Gen.matrix (congruenceDiagonalGramian p)) : ("congruenceAdjoint", checkForAll ((,) <$> Gen.matrix <#*#> genHermitian p) congruenceAdjoint) : ("congruenceCongruenceDiagonal", checkForAllExtra genOrder ((,) <$> Gen.vectorReal <-*#> Gen.matrix) (congruenceCongruenceDiagonal p)) : ("anticommutator", checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix) (anticommutator p)) : ("anticommutatorCommutative", checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix) (anticommutatorCommutative p)) : ("anticommutatorAdjoint", checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix) (anticommutatorAdjoint p)) : ("outer", checkForAllExtra genOrder Gen.vector (outer p)) : ("sumRank1", checkForAllExtra genOrder genScaledVectors (sumRank1 p)) : ("sumRank1NonEmpty", checkForAllExtra genOrder (snd <$> genScaledVectors) (sumRank1NonEmpty p)) : ("sumRank2", checkForAllExtra genOrder genScaledVectorPairs (sumRank2 p)) : ("sumRank2NonEmpty", checkForAllExtra genOrder (snd <$> genScaledVectorPairs) (sumRank2NonEmpty p)) : ("addAdjoint", checkForAll Gen.square (addAdjoint p)) : Mosaic.testsVar Mosaic.Hermitian p ++ ("determinant", checkForAll (genHermitian p) Divide.determinant) : ("choleskyQR", checkForAll Gen.tall (choleskyQR p)) : ("solve", checkForAll ((,) <$> genInvertible p <#\#> Gen.matrix) Divide.solve) : ("solveIdentity", checkForAll ((,) <$> (repack p <$> Gen.identity `asTypeOf` Gen.hermitian) <#\#> Gen.matrix) Divide.solveIdentity) : ("inverse", checkForAll (genInvertible p) Divide.inverse) : Divide.testsVar (genInvertible p) ++ ("determinantPD", checkForAll (genPositiveDefinite p) determinantPD) : ("determinantND", checkForAll (genNegativeDefinite p) Divide.determinant) : ("inversePD", checkForAll (genPositiveDefinite p) inversePD) : ("inverseND", checkForAll (genNegativeDefinite p) Divide.inverse) : ("solvePD", checkForAll ((,) <$> genPositiveDefinite p <#\#> Gen.matrix) solvePD) : ("solveND", checkForAll ((,) <$> genNegativeDefinite p <#\#> Gen.matrix) Divide.solve) : ("solveDecomposedPD", checkForAll ((,) <$> genPositiveDefinite p <#\#> Gen.matrix) solveDecomposedPD) : Util.suffix "PD" (Divide.testsVar (genPositiveDefinite p)) ++ Util.suffix "ND" (Divide.testsVar (genNegativeDefinite p)) ++ ("eigenvaluesDeterminant", checkForAll (genHermitian p) eigenvaluesDeterminant) : ("eigensystem", checkForAll (genHermitian p) eigensystem) : []