{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE ExistentialQuantification #-} {-# LANGUAGE GADTs #-} module Test.BandedHermitian (testsVar) where import qualified Test.Divide as Divide import qualified Test.Indexed as Indexed import qualified Test.Generator as Gen import qualified Test.Utility as Util import Test.Banded.Utility (Square(Square), genSquare, shapeBandedFromFull, natFromProxy, offDiagonalNats) import Test.Generator ((<-*#>), (<#*|>), (<-*|>), (<#*#>), (<#\#>)) import Test.Logic (Dim) import Test.Utility (approxReal, approxArray, approxRealVectorTol, approxMatrix, approxVector, genOrder, genArray, genVector, Tagged, equalListWith) import qualified Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite as BandedHermitianPD import qualified Numeric.LAPACK.Matrix.BandedHermitian as BandedHermitian import qualified Numeric.LAPACK.Matrix.Banded as Banded import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent as Extent import qualified Numeric.LAPACK.Matrix.Square as Square import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix import qualified Numeric.LAPACK.Matrix as Matrix import qualified Numeric.LAPACK.Vector as Vector import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic import Numeric.LAPACK.Matrix (ShapeInt, shapeInt, (#*#), (#*##), (-*#), (#*|)) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf, fromReal, absolute, selectReal) import qualified Numeric.Netlib.Class as Class import qualified Type.Data.Num.Unary.Proof as Proof import qualified Type.Data.Num.Unary as Unary import Type.Data.Num.Unary (unary) import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape import Foreign.Storable (Storable) import Control.Applicative (liftA2, (<$>)) import qualified Data.List.HT as ListHT import Data.Tuple.HT (mapSnd) import Data.Semigroup ((<>)) import qualified Test.QuickCheck as QC data BandedHermitian size a = forall offDiag. (Unary.Natural offDiag) => BandedHermitian (BandedHermitian.BandedHermitian offDiag size a) instance (Show size, Show a, Shape.C size, Storable a) => Show (BandedHermitian size a) where showsPrec p (BandedHermitian a) = showsPrec p a shapeBandedHermitianFromSquare :: (MatrixShape.UnaryProxy off) -> MatrixShape.Square size -> MatrixShape.BandedHermitian off size shapeBandedHermitianFromSquare k (MatrixShape.Full order extent) = MatrixShape.BandedHermitian k order $ Extent.height extent {- Non-real elements on the diagonal. -} _genBandedHermitian :: (Class.Floating a) => Gen.MatrixInt a (BandedHermitian ShapeInt a) _genBandedHermitian = flip Gen.mapGenDim Gen.squareShape $ \maxElem maxDim shape -> do k <- QC.choose (0, toInteger maxDim) Unary.reifyNatural k $ \numOff -> fmap BandedHermitian $ genArray maxElem $ shapeBandedHermitianFromSquare (unary numOff) shape genBandedHermitian :: (Dim sh, Shape.Indexed sh, Shape.Index sh ~ ix, Eq ix, Class.Floating a, RealOf a ~ ar, Class.Real ar) => Gen.Square sh a (BandedHermitian sh a) genBandedHermitian = flip Gen.mapGenDim Gen.squareShape $ \maxElem maxDim sqShape -> do k <- QC.choose (0, toInteger maxDim) Unary.reifyNatural k $ \numOff -> do let shape = shapeBandedHermitianFromSquare (unary numOff) sqShape BandedHermitian <$> (Util.genArrayIndexed shape $ \ix -> let real = case ix of MatrixShape.InsideBox r c -> r==c MatrixShape.VertOutsideBox _ _ -> False MatrixShape.HorizOutsideBox _ _ -> False in if real then fromReal <$> Util.genReal maxElem else Util.genElement maxElem) convertToFull :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian ShapeInt a -> Bool convertToFull (BandedHermitian a) = approxArray (Hermitian.toSquare $ BandedHermitian.toHermitian a) (Banded.toFull $ BandedHermitian.toBanded a) takeDiagonal :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian ShapeInt a -> Bool takeDiagonal (BandedHermitian a) = approxRealVectorTol 1e-5 (Hermitian.takeDiagonal $ BandedHermitian.toHermitian a) (BandedHermitian.takeDiagonal a) gramian :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Square ShapeInt a -> Bool gramian (Square a) = let (sub,super) = offDiagonalNats a in case (Proof.addNat sub super, Proof.addComm sub super) of (Proof.Nat, Proof.AddComm) -> approxArray (BandedHermitian.toBanded $ BandedHermitian.gramian a) (Banded.adjoint a #*# a) type StaticVector1 n = Vector (ShapeStatic.ZeroBased (Unary.Succ n)) data SumRank1 size a = forall offDiag. (Unary.Natural offDiag) => SumRank1 size [(RealOf a, (Shape.Index size, StaticVector1 offDiag a))] instance (Show size, Show (Shape.Index size), Show a, Show (RealOf a), Shape.C size, Storable a) => Show (SumRank1 size a) where showsPrec p (SumRank1 sh a) = showsPrec p (sh,a) genScaledVectors :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Gen.VectorInt a (SumRank1 ShapeInt a) genScaledVectors = flip Gen.mapGen Gen.vectorDim $ \maxElem size@(Shape.ZeroBased n) -> do k <- QC.choose (0, n-1) Unary.reifyNatural (toInteger k) $ \numOff -> fmap (SumRank1 size) $ if n==0 then return [] else QC.listOf $ liftA2 (,) (Util.genReal maxElem) $ liftA2 (,) (QC.choose (0,n-k-1)) (genVector maxElem (ShapeStatic.ZeroBased $ unary $ Unary.succ numOff)) sumRank1 :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => MatrixShape.Order -> SumRank1 ShapeInt a -> Bool sumRank1 order (SumRank1 sh xs) = approxArray (BandedHermitian.toHermitian $ BandedHermitian.sumRank1 order sh xs) (Hermitian.sumRank1 order sh $ map (mapSnd (uncurry $ displace sh)) xs) displace :: (Shape.C sh, Class.Floating a) => ShapeInt -> Int -> Vector sh a -> Vector ShapeInt a displace (Shape.ZeroBased n) k a = Array.mapShape (shapeInt . Shape.size) $ Vector.zero (shapeInt k) `Vector.append` a `Vector.append` Vector.zero (shapeInt $ max 0 $ n - k - Shape.size (Array.shape a)) multiplyIdentity :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Hermitian.Transposition -> Matrix.General ShapeInt ShapeInt a -> Bool multiplyIdentity trans m = approxArray m (BandedHermitian.multiplyFull trans (BandedHermitian.identity (Matrix.height m)) m) multiplyDiagonal :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Hermitian.Transposition -> (Vector ShapeInt ar, Matrix.General ShapeInt ShapeInt a) -> Bool multiplyDiagonal trans (d,m) = approxArray (Matrix.scaleRowsReal d m) (BandedHermitian.multiplyFull trans (BandedHermitian.diagonal d) m) multiplyFullIdentity :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian ShapeInt a -> Bool multiplyFullIdentity (BandedHermitian m) = let a = Banded.toFull $ BandedHermitian.toBanded m in approxArray a $ BandedHermitian.multiplyFull BandedHermitian.NonTransposed m $ Square.identityFrom a multiplyHermitianVector :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Hermitian.Transposition -> (BandedHermitian ShapeInt a, Vector ShapeInt a) -> Bool multiplyHermitianVector trans (BandedHermitian m, x) = approxVector (BandedHermitian.multiplyVector trans m x) (Hermitian.multiplyVector trans (BandedHermitian.toHermitian m) x) multiplyVectorDot :: (Class.Floating a, Eq a) => (Vector ShapeInt a, BandedHermitian ShapeInt a, Vector ShapeInt a) -> Bool multiplyVectorDot (x, BandedHermitian m, y) = Vector.dot x (m#*|y) == Vector.dot (x-*#m) y multiplyFullAny :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Hermitian.Transposition -> (BandedHermitian ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool multiplyFullAny trans (BandedHermitian a, b) = approxArray (BandedHermitian.multiplyFull trans a b) (Hermitian.multiplyFull trans (BandedHermitian.toHermitian a) b) multiplyFullColumns :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian.Transposition -> (BandedHermitian ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool multiplyFullColumns trans (BandedHermitian a, b) = equalListWith approxVector (Matrix.toColumns (BandedHermitian.multiplyFull trans a b)) (map (BandedHermitian.multiplyVector trans a) (Matrix.toColumns b)) multiplyFullAssoc :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian.Transposition -> (BandedHermitian ShapeInt a, Matrix.General ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool multiplyFullAssoc trans (BandedHermitian a, b, c) = approxArray (Matrix.multiply (BandedHermitian.multiplyFull trans a b) c) (BandedHermitian.multiplyFull trans a (Matrix.multiply b c)) genBandedHPD :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => Gen.MatrixInt a (BandedHermitian ShapeInt a) genBandedHPD = flip Gen.mapGenDim Gen.squareShape $ \maxElem maxDim shape -> do kl <- QC.choose (0, toInteger maxDim) ku <- QC.choose (0, toInteger maxDim) Unary.reifyNatural kl $ \subU -> Unary.reifyNatural ku $ \superU -> let sub = unary subU; subP = natFromProxy sub super = unary superU; superP = natFromProxy super in case (Proof.addNat subP superP, Proof.addComm subP superP) of (Proof.Nat, Proof.AddComm) -> fmap (BandedHermitian . BandedHermitian.gramian) $ (genArray maxElem $ shapeBandedFromFull (sub, super) shape) `QC.suchThat` (\a -> absolute (Banded.determinant a) > 0.1) determinant :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian ShapeInt a -> Bool determinant (BandedHermitian a) = let detB = BandedHermitianPD.determinant a detS = Hermitian.determinant $ BandedHermitian.toHermitian a in approxReal (selectReal 1 1e-3 * max 1 (abs detB + abs detS)) detB detS multiplySolve :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => (BandedHermitian ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool multiplySolve (BandedHermitian a, b) = approxMatrix (selectReal 10 1e-3) (a #*## BandedHermitianPD.solve a b) b solveDecomposed :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => (BandedHermitian ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool solveDecomposed (BandedHermitian a, b) = approxMatrix (selectReal 1e-3 1e-7) (BandedHermitianPD.solve a b) (BandedHermitianPD.solveDecomposed (BandedHermitianPD.decompose a) b) eigenvaluesDeterminant :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian ShapeInt a -> Bool eigenvaluesDeterminant (BandedHermitian a) = let det = BandedHermitianPD.determinant a prod = Vector.product $ BandedHermitian.eigenvalues a in approxReal ((det+prod) * selectReal 0.5 1e-6) det prod eigensystem :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian ShapeInt a -> Bool eigensystem (BandedHermitian a) = let (q,d) = BandedHermitian.eigensystem a in approxMatrix 1e-4 (Banded.toFull $ BandedHermitian.toBanded a) (q <> Matrix.scaleRowsReal d (Square.adjoint q)) eigenvaluesHermitian :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian ShapeInt a -> Bool eigenvaluesHermitian (BandedHermitian a) = approxRealVectorTol (selectReal 1e-3 1e-5) (BandedHermitian.eigenvalues a) (Hermitian.eigenvalues $ BandedHermitian.toHermitian a) eigensystemHermitian :: (Class.Floating a, RealOf a ~ ar, Class.Real ar) => BandedHermitian ShapeInt a -> QC.Property eigensystemHermitian (BandedHermitian a) = let (q0,d0) = BandedHermitian.eigensystem a (q1,d1) = Hermitian.eigensystem $ BandedHermitian.toHermitian a unit = Matrix.adjoint q0 <> q1 tol = selectReal 1e-4 1e-7 in not (or (ListHT.mapAdjacent (approxReal 0.1) (Array.toList d0))) QC.==> approxRealVectorTol tol d0 d1 && and (zipWith (\(r,c) x -> approxReal tol (absolute x) $ if r==c then 1 else 0) (Shape.indices $ ArrMatrix.shape unit) (Array.toList $ ArrMatrix.toVector unit)) checkForAll :: (Show a, QC.Testable test) => Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property checkForAll gen = Util.checkForAll (Gen.run gen 6 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 = ("index", checkForAll (Indexed.genMatrixIndexGen (\(BandedHermitian arr) -> Matrix.indices arr) genBandedHermitian) (\(mix, BandedHermitian arr) -> Indexed.unitDot (mix, arr))) : ("convertToFull", checkForAll genBandedHermitian convertToFull) : ("takeDiagonal", checkForAll genBandedHermitian takeDiagonal) : ("sumRank1", checkForAllExtra genOrder genScaledVectors sumRank1) : ("gramian", checkForAll genSquare gramian) : ("multiplyIdentity", checkForAllExtra QC.arbitraryBoundedEnum Gen.matrix multiplyIdentity) : ("multiplyDiagonal", checkForAllExtra QC.arbitraryBoundedEnum ((,) <$> Gen.vectorReal <-*#> Gen.matrix) multiplyDiagonal) : ("multiplyFullIdentity", checkForAll genBandedHermitian multiplyFullIdentity) : ("multiplyFullAny", checkForAllExtra QC.arbitraryBoundedEnum ((,) <$> genBandedHermitian <#*#> Gen.matrix) multiplyFullAny) : ("multiplyHermitianVector", checkForAllExtra QC.arbitraryBoundedEnum ((,) <$> genBandedHermitian <#*|> Gen.vector) multiplyHermitianVector) : ("multiplyVectorDot", checkForAll ((,,) <$> Gen.vector <-*#> genBandedHermitian <-*|> Gen.vector) multiplyVectorDot) : ("multiplyFullColumns", checkForAllExtra QC.arbitraryBoundedEnum ((,) <$> genBandedHermitian <#*#> Gen.matrix) multiplyFullColumns) : ("multiplyFullAssoc", checkForAllExtra QC.arbitraryBoundedEnum ((,,) <$> genBandedHermitian <#*#> Gen.matrix <#*#> Gen.matrix) multiplyFullAssoc) : ("determinant", checkForAll genBandedHPD determinant) : ("multiplySolve", checkForAll ((,) <$> genBandedHPD <#\#> Gen.matrix) multiplySolve) : ("solveDecomposed", checkForAll ((,) <$> genBandedHPD <#\#> Gen.matrix) solveDecomposed) : map (mapSnd ($ ((\(BandedHermitian m) -> Divide.SquareMatrix m) <$> genBandedHPD))) Divide.testsVarAny ++ ("eigenvaluesDeterminant", checkForAll genBandedHPD eigenvaluesDeterminant) : ("eigensystem", checkForAll genBandedHermitian eigensystem) : ("eigenvaluesHermitian", checkForAll genBandedHermitian eigenvaluesHermitian) : ("eigensystemHermitian", checkForAll genBandedHermitian eigensystemHermitian) : []