lapack-0.4: test/Test/Hermitian.hs
{-# 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) :
("addDistributive",
checkForAll
(Generic.genDistribution (genHermitian p))
Generic.addDistributive) :
("subDistributive",
checkForAll
(Generic.genDistribution (genHermitian p))
Generic.subDistributive) :
("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) :
[]