{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE Rank2Types #-} module Test.Triangular (testsVar) where import qualified Test.Divide as Divide import qualified Test.Multiply as Multiply import qualified Test.Generic as Generic import qualified Test.Indexed as Indexed import qualified Test.Generator as Gen import qualified Test.Utility as Util import Test.Generator ((<-*#>), (<#*|>), (<#*#>), (<#\#>)) import Test.Utility (approx, approxArray, approxArrayTol, approxMatrix, approxVector, equalArray, Tagged, (!|||), (!===)) import qualified Numeric.LAPACK.Matrix.Triangular as Triangular import qualified Numeric.LAPACK.Matrix.Square as Square import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix import qualified Numeric.LAPACK.Matrix as Matrix import qualified Numeric.LAPACK.Vector as Vector import Numeric.LAPACK.Matrix.Triangular (Triangular) 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 Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable ((!)) import Data.Array.Comfort.Shape ((:+:)) import Control.Applicative ((<$>)) import Data.Tuple.HT (mapFst, uncurry3) import Data.Semigroup ((<>)) import qualified Test.QuickCheck as QC -- cf. Test.Generic.addDistributive addDistributive :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Eq lo, Eq up, Class.Floating a, RealOf a ~ ar, Class.Real ar) => ((Triangular lo diag up ShapeInt a, Triangular lo diag up ShapeInt a), Vector ShapeInt a) -> Bool addDistributive ((a,b),x) = approxVector ((Triangular.strictNonUnitDiagonal a #+# Triangular.strictNonUnitDiagonal b) #*| x) (a#*|x |+| b#*|x) subDistributive :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Eq lo, Eq up, Class.Floating a, RealOf a ~ ar, Class.Real ar) => ((Triangular lo diag up ShapeInt a, Triangular lo diag up ShapeInt a), Vector ShapeInt a) -> Bool subDistributive ((a,b),x) = approxVector ((Triangular.strictNonUnitDiagonal a #-# Triangular.strictNonUnitDiagonal b) #*| x) (a#*|x |-| b#*|x) expandTriangle :: (MatrixShape.Content lo, MatrixShape.TriDiag diag, MatrixShape.Content up, Shape.C sh, Class.Floating a) => Triangular lo diag up sh a -> General sh sh a expandTriangle = Matrix.fromFull . Triangular.toSquare transposedZero :: (Class.Floating a, Shape.C width, Shape.C height) => General height width a -> General width height a transposedZero = ArrMatrix.zero . ArrMatrix.shape . Matrix.transpose stackDiagonal :: (MatrixShape.TriDiag diag, Class.Floating a) => (Triangular.FlexDiagonal diag ShapeInt a, Triangular.FlexDiagonal diag ShapeInt a) -> Bool stackDiagonal (a,c) = let ac = expandTriangle $ Triangular.stackDiagonal a c b = Matrix.zero $ MatrixShape.general MatrixShape.RowMajor (Matrix.height a) (Matrix.height c) in equalArray ac $ (expandTriangle a ||| b === Matrix.transpose b ||| expandTriangle c) stackLower :: (MatrixShape.TriDiag diag, Class.Floating a) => (Triangular.FlexLower diag ShapeInt a, General ShapeInt ShapeInt a, Triangular.FlexLower diag ShapeInt a) -> Bool stackLower (a,b,c) = let abc = expandTriangle $ Triangular.stackLower a b c in equalArray abc $ (expandTriangle a !||| transposedZero b === b !||| expandTriangle c) stackUpper :: (MatrixShape.TriDiag diag, Class.Floating a) => (Triangular.FlexUpper diag ShapeInt a, General ShapeInt ShapeInt a, Triangular.FlexUpper diag ShapeInt a) -> Bool stackUpper (a,b,c) = let abc = Matrix.fromFull $ Triangular.toSquare $ Triangular.stackUpper a b c in equalArray abc $ (expandTriangle a ||| b !=== transposedZero b ||| expandTriangle c) stackSymmetric :: (MatrixShape.TriDiag diag, Class.Floating a) => (Triangular.FlexSymmetric diag ShapeInt a, General ShapeInt ShapeInt a, Triangular.FlexSymmetric diag ShapeInt a) -> Bool stackSymmetric (a,b,c) = let abc = expandTriangle $ Triangular.stackSymmetric a b c in equalArray abc $ (expandTriangle a ||| b !=== Matrix.transpose b ||| expandTriangle c) splitDiagonal :: (MatrixShape.TriDiag diag, Eq diag, Class.Floating a) => Triangular.FlexDiagonal diag (ShapeInt:+:ShapeInt) a -> Bool splitDiagonal abc = equalArray abc $ uncurry Triangular.stackDiagonal $ Triangular.splitDiagonal abc splitLower :: (MatrixShape.TriDiag diag, Eq diag, Class.Floating a) => Triangular.FlexLower diag (ShapeInt:+:ShapeInt) a -> Bool splitLower abc = equalArray abc $ uncurry3 Triangular.stackLower $ Triangular.splitLower abc splitUpper :: (MatrixShape.TriDiag diag, Eq diag, Class.Floating a) => Triangular.FlexUpper diag (ShapeInt:+:ShapeInt) a -> Bool splitUpper abc = equalArray abc $ uncurry3 Triangular.stackUpper $ Triangular.splitUpper abc splitSymmetric :: (MatrixShape.TriDiag diag, Eq diag, Class.Floating a) => Triangular.FlexSymmetric diag (ShapeInt:+:ShapeInt) a -> Bool splitSymmetric abc = equalArray abc $ uncurry3 Triangular.stackSymmetric $ Triangular.splitSymmetric abc multiplyIdentityVector :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, Vector ShapeInt a) -> Bool multiplyIdentityVector (eye,a) = approxVector a (Triangular.multiplyVector eye a) multiplyIdentityFull :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, General ShapeInt ShapeInt a) -> Bool multiplyIdentityFull (eye,a) = approxArray a (Triangular.multiplyFull eye a) multiplyIdentity :: (MatrixShape.DiagUpLo lo up, MatrixShape.TriDiag diag, Eq lo, Eq diag, Eq up, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, Triangular lo diag up ShapeInt a) -> Bool multiplyIdentity (eye,a) = approxArray a (eye <> a) multiplyVector :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, Vector ShapeInt a) -> Bool multiplyVector (a,x) = approxVector (Triangular.toSquare a #*| x) (Triangular.multiplyVector a x) multiply :: (MatrixShape.DiagUpLo lo up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, Triangular lo diag up ShapeInt a) -> Bool multiply (a,b) = approxArray (Triangular.toSquare a <> Triangular.toSquare b) (Triangular.toSquare $ a <> b) multiplyFull :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, General ShapeInt ShapeInt a) -> Bool multiplyFull (a,b) = approxArray (Triangular.toSquare a #*## b) (Triangular.multiplyFull a b) multiplyVectorLeft :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Vector ShapeInt a, Triangular lo diag up ShapeInt a) -> Bool multiplyVectorLeft (x,a) = approxVector (x -*# Triangular.toSquare a) (x -*# a) multiplyVectorRight :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, Vector ShapeInt a) -> Bool multiplyVectorRight (a,x) = approxVector (Triangular.toSquare a #*| x) (a #*| x) multiplyLeft :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (General ShapeInt ShapeInt a, Triangular lo diag up ShapeInt a) -> Bool multiplyLeft (a,b) = approxMatrix 1e-5 (a ##*# Triangular.toSquare b) (a ##*# b) multiplyRight :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, General ShapeInt ShapeInt a) -> Bool multiplyRight (a,b) = approxArray (Triangular.toSquare a #*## b) (a #*## b) genInvertible :: (MatrixShape.Content up, MatrixShape.Content lo, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => GenTriangular lo diag up ShapeInt a genInvertible = Gen.condition Util.invertible Gen.triangular inverse :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => Triangular lo diag up ShapeInt a -> Bool inverse a = approxArrayTol (selectReal 1 1e-5) (Triangular.toSquare $ Triangular.inverse a) (Square.inverse $ Triangular.toSquare a) solve :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool solve (a, b) = approxMatrix (selectReal 1 1e-5) (Triangular.solve a b) (Square.solve (Triangular.toSquare a) b) solveIdentity :: (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => (Triangular lo diag up ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool solveIdentity (eye, a) = approxMatrix (selectReal 1e-3 1e-5) a (Triangular.solve eye a) eigenvaluesDeterminant :: (MatrixShape.DiagUpLo lo up, MatrixShape.TriDiag diag, Class.Floating a, RealOf a ~ ar, Class.Real ar) => Triangular lo diag up ShapeInt a -> Bool eigenvaluesDeterminant a = approx (selectReal 1e-1 1e-5) (Triangular.determinant a) (Vector.product $ Triangular.eigenvalues a) genDiagonalizable :: (MatrixShape.Content lo, MatrixShape.Content up, Class.Floating a, RealOf a ~ ar, Class.Real ar) => GenTriangular lo MatrixShape.NonUnit up ShapeInt a genDiagonalizable = flip Gen.mapGen Gen.triangularShape $ \maxElem shape -> do d <- Util.genDistinct 3 10 $ MatrixShape.triangularSize shape Util.genArrayExtraDiag maxElem shape (\r -> return (d!r)) eigensystem :: (MatrixShape.DiagUpLo lo up, Eq lo, Eq up, Class.Floating a, RealOf a ~ ar, Class.Real ar) => MatrixShape.Order -> Triangular lo MatrixShape.NonUnit up ShapeInt a -> Bool eigensystem order a = let (vr,d,vl) = Triangular.eigensystem a scal = Triangular.takeDiagonal $ vl <> vr in approxMatrix (selectReal 1e-3 1e-5) a (vr <> Triangular.diagonal order (Vector.divide d scal) <> vl) 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 type GenTriangular lo diag up sh a = Gen.Square sh a (Triangular lo diag up sh a) addSuperName :: String -> [(String, a)] -> [(String, a)] addSuperName superName = map (mapFst ((superName++) . ("."++))) checkAnyFlexDiag :: (MatrixShape.TriDiag diag) => String -> (forall lo up. (MatrixShape.Content lo, MatrixShape.Content up, Eq lo, Eq up, Show lo, Show up) => GenTriangular lo diag up sh a -> Tagged a QC.Property) -> (forall lo up. (MatrixShape.Content lo, MatrixShape.Content up, Eq lo, Eq up, Show lo, Show up) => GenTriangular lo diag up sh a) -> [(String, Tagged a QC.Property)] checkAnyFlexDiag name checker gen = (checkDiagUpLoFlexDiag name checker gen ++) $ addSuperName name $ ("Symmetric", checker (Triangular.asSymmetric <$> gen)) : [] checkAny :: String -> (forall lo up diag. (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Eq lo, Eq up, Show lo, Show up, Show diag) => GenTriangular lo diag up sh a -> Tagged a QC.Property) -> (forall lo up diag. (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Eq lo, Eq up, Show lo, Show up, Show diag) => GenTriangular lo diag up sh a) -> [(String, Tagged a QC.Property)] checkAny name checker gen = checkAnyFlexDiag (name++".Unit") checker (Triangular.requireUnitDiagonal <$> gen) ++ checkAnyFlexDiag (name++".NonUnit") checker (Triangular.requireNonUnitDiagonal <$> gen) checkDiagUpLoFlexDiag :: (MatrixShape.TriDiag diag) => String -> (forall lo up. (MatrixShape.DiagUpLo lo up, Eq lo, Eq up, Show lo, Show up) => GenTriangular lo diag up sh a -> Tagged a QC.Property) -> (forall lo up. (MatrixShape.DiagUpLo lo up, Eq lo, Eq up, Show lo, Show up) => GenTriangular lo diag up sh a) -> [(String, Tagged a QC.Property)] checkDiagUpLoFlexDiag name checker gen = addSuperName name $ ("Diagonal", checker (Triangular.asDiagonal <$> gen)) : ("Lower", checker (Triangular.asLower <$> gen)) : ("Upper", checker (Triangular.asUpper <$> gen)) : [] checkDiagUpLo :: String -> (forall lo up diag. (MatrixShape.DiagUpLo lo up, MatrixShape.TriDiag diag, Eq lo, Eq diag, Eq up, Show lo, Show diag, Show up) => GenTriangular lo diag up sh a -> Tagged a QC.Property) -> (forall lo up diag. (MatrixShape.DiagUpLo lo up, MatrixShape.TriDiag diag, Eq lo, Eq diag, Eq up, Show lo, Show diag, Show up) => GenTriangular lo diag up sh a) -> [(String, Tagged a QC.Property)] checkDiagUpLo name checker gen = checkDiagUpLoFlexDiag (name++".Unit") checker (Triangular.requireUnitDiagonal <$> gen) ++ checkDiagUpLoFlexDiag (name++".NonUnit") checker (Triangular.requireNonUnitDiagonal <$> gen) newtype Power diag sh a lo up = Power {getPower :: GenTriangular lo diag up sh a -> Tagged a QC.Property} restrictDiagUpLo :: (MatrixShape.DiagUpLo lo0 up0, MatrixShape.TriDiag diag0, Eq lo0, Eq diag0, Eq up0, Show lo0, Show diag0, Show up0) => (forall lo up diag. (Triangular.PowerContentDiag lo diag up, Eq lo, Eq diag, Eq up, Show lo, Show diag, Show up) => GenTriangular lo diag up sh a -> Tagged a QC.Property) -> GenTriangular lo0 diag0 up0 sh a -> Tagged a QC.Property restrictDiagUpLo f = getPower $ MatrixShape.switchDiagUpLo (Power f) (Power f) (Power f) checkDiagUpLoSym :: String -> (forall lo up diag. (Triangular.PowerContentDiag lo diag up, Eq lo, Eq diag, Eq up, Show lo, Show diag, Show up) => GenTriangular lo diag up ShapeInt a -> Tagged a QC.Property) -> (forall lo up diag. (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag, Eq lo, Eq diag, Eq up, Show lo, Show diag, Show up) => GenTriangular lo diag up ShapeInt a) -> [(String, Tagged a QC.Property)] checkDiagUpLoSym name checker gen = (checkDiagUpLo name (restrictDiagUpLo checker) gen ++) $ addSuperName name $ ("Symmetric", checker (Triangular.asSymmetric <$> gen)) : [] checkFlexDiag :: String -> (forall diag. (MatrixShape.TriDiag diag, Eq diag, Show diag) => GenTriangular lo diag up sh a -> Tagged a QC.Property) -> (forall diag. (MatrixShape.TriDiag diag, Eq diag, Show diag) => GenTriangular lo diag up sh a) -> [(String, Tagged a QC.Property)] checkFlexDiag name checker gen = (name++".Unit", checker (Triangular.requireUnitDiagonal <$> gen)) : (name++".NonUnit", checker (Triangular.requireNonUnitDiagonal <$> gen)) : [] testsVar :: (Show a, Show ar, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) => [(String, Tagged a QC.Property)] testsVar = checkAny "index" (\gen -> checkForAll (Indexed.genMatrixIndex gen) Indexed.unitDot) Gen.triangular ++ checkFlexDiag "stackDiagonal" (\gen -> checkForAll (Gen.stackDiagonal gen gen) stackDiagonal) Gen.triangular ++ checkFlexDiag "stackLower" (\gen -> checkForAll (Gen.stack3 gen (Matrix.transpose <$> Gen.matrix) gen) stackLower) Gen.triangular ++ checkFlexDiag "stackUpper" (\gen -> checkForAll (Gen.stack3 gen Gen.matrix gen) stackUpper) Gen.triangular ++ checkFlexDiag "stackSymmetric" (\gen -> checkForAll (Gen.stack3 gen Gen.matrix gen) stackSymmetric) Gen.triangular ++ checkFlexDiag "splitDiagonal" (\gen -> checkForAll gen splitDiagonal) Gen.triangular ++ checkFlexDiag "splitLower" (\gen -> checkForAll gen splitLower) Gen.triangular ++ checkFlexDiag "splitUpper" (\gen -> checkForAll gen splitUpper) Gen.triangular ++ checkFlexDiag "splitSymmetric" (\gen -> checkForAll gen splitSymmetric) Gen.triangular ++ checkAny "forceOrder" (\gen -> checkForAllExtra Util.genOrder ((,) <$> gen <#*|> Gen.vector) Generic.forceOrder) Gen.triangular ++ checkAny "addDistributive" (\gen -> checkForAll (Generic.genDistribution gen) addDistributive) Gen.triangular ++ checkAny "subDistributive" (\gen -> checkForAll (Generic.genDistribution gen) subDistributive) Gen.triangular ++ checkAny "multiplyIdentityVector" (\gen -> checkForAll ((,) <$> gen <#*|> Gen.vector) multiplyIdentityVector) (Triangular.relaxUnitDiagonal <$> Gen.identity) ++ checkAny "multiplyIdentityFull" (\gen -> checkForAll ((,) <$> gen <#*#> Gen.matrix) multiplyIdentityFull) (Triangular.relaxUnitDiagonal <$> Gen.identity) ++ checkDiagUpLo "multiplyIdentity" (\gen -> checkForAll ((,) <$> gen <#*#> Gen.triangular) multiplyIdentity) (Triangular.relaxUnitDiagonal <$> Gen.identity) ++ checkAny "multiplyVector" (\gen -> checkForAll ((,) <$> gen <#*|> Gen.vector) multiplyVector) Gen.triangular ++ checkAny "multiplyFull" (\gen -> checkForAll ((,) <$> gen <#*#> Gen.matrix) multiplyFull) Gen.triangular ++ checkAny "multiplyVectorLeft" (\gen -> checkForAll ((,) <$> Gen.vector <-*#> gen) multiplyVectorLeft) Gen.triangular ++ checkAny "multiplyVectorRight" (\gen -> checkForAll ((,) <$> gen <#*|> Gen.vector) multiplyVectorRight) Gen.triangular ++ checkAny "multiplyLeft" (\gen -> checkForAll ((,) <$> Gen.matrix <#*#> gen) multiplyLeft) Gen.triangular ++ checkAny "multiplyRight" (\gen -> checkForAll ((,) <$> gen <#*#> Gen.matrix) multiplyRight) Gen.triangular ++ checkDiagUpLo "multiply" (\gen -> checkForAll ((,) <$> gen <#*#> gen) multiply) Gen.triangular ++ checkDiagUpLoSym "multiplySquare" (\gen -> checkForAll gen Multiply.multiplySquare) Gen.triangular ++ checkDiagUpLoSym "squareSquare" (\gen -> checkForAll gen Multiply.squareSquare) Gen.triangular ++ checkDiagUpLoSym "power" (\gen -> checkForAllExtra (QC.choose (0,10::Int)) gen Multiply.power) Gen.triangular ++ checkAny "determinant" (\gen -> checkForAll gen Divide.determinant) Gen.triangular ++ checkAny "solve" (\gen -> checkForAll ((,) <$> gen <#\#> Gen.matrix) solve) genInvertible ++ checkAny "solveIdentity" (\gen -> checkForAll ((,) <$> gen <#\#> Gen.matrix) solveIdentity) (Triangular.relaxUnitDiagonal <$> Gen.identity) ++ checkAny "inverse" (\gen -> checkForAll gen inverse) genInvertible ++ concatMap (\(name,test) -> checkAny name (test . fmap Divide.SquareMatrix) genInvertible) Divide.testsVarAny ++ checkDiagUpLo "eigenvaluesDeterminant" (\gen -> checkForAll gen eigenvaluesDeterminant) Gen.triangular ++ checkDiagUpLoFlexDiag "eigensystem" (\gen -> checkForAllExtra Util.genOrder gen eigensystem) genDiagonalizable ++ []