Singular.hs

{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
module Test.Singular (testsVar) where

import qualified Test.Generator as Gen
import qualified Test.Utility as Util
import Test.Generator ((<#\#>))
import Test.Utility
         (approxReal, approxArrayTol, approxMatrix, isUnitary, Tagged)

import qualified Numeric.LAPACK.Singular as Singular
import qualified Numeric.LAPACK.Orthogonal as Ortho
import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix as Matrix
import Numeric.LAPACK.Matrix
         (General, ShapeInt, shapeInt, (##*#), (#*#), (#*##))
import Numeric.LAPACK.Scalar (RealOf, selectReal)

import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape

import Control.Applicative ((<$>))
import Data.Semigroup ((<>))

import qualified Test.QuickCheck as QC


pseudoInverseOrtho ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   General ShapeInt ShapeInt a -> Bool
pseudoInverseOrtho a =
   let (no,invo) = Ortho.pseudoInverseRCond 1e-5 a
       (ns,invs) = Singular.pseudoInverseRCond 1e-5 a
       tol = selectReal 1e-2 1e-5
   in no==ns && approxMatrix tol invo invs

pseudoInverseProjection ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   General ShapeInt ShapeInt a -> Bool
pseudoInverseProjection a =
   let ainv = snd $ Singular.pseudoInverseRCond 1e-5 a
       tol = selectReal 1e-1 1e-5
   in approxArrayTol tol a (a <> ainv <> a) &&
      approxArrayTol tol ainv (ainv <> a <> ainv)

pseudoInverseHermitian ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   General ShapeInt ShapeInt a -> Bool
pseudoInverseHermitian a =
   let ainv = snd $ Singular.pseudoInverseRCond 1e-5 a
       tol = selectReal 1e-2 1e-5
       aainv = a <> ainv
       ainva = ainv <> a
   in approxMatrix tol aainv (Matrix.adjoint aainv) &&
      approxMatrix tol ainva (Matrix.adjoint ainva)


determinantAbsolute ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   General ShapeInt ShapeInt a -> Bool
determinantAbsolute a =
   let detOrtho = Ortho.determinantAbsolute a
       detSing = Singular.determinantAbsolute a
   in approxReal
         (selectReal 1e-3 1e-5 * max 1 (max detOrtho detSing))
         detOrtho detSing


leastSquaresMinimumNorm ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.General ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) ->
   Bool
leastSquaresMinimumNorm (a,b) =
   let (no,xo) = Ortho.leastSquaresMinimumNormRCond 1e-5 a b
       (ns,xs) = Singular.leastSquaresMinimumNormRCond 1e-5 a b
   in no==ns &&
      approxMatrix (selectReal 10 1e-3) xo xs


decompose ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.General ShapeInt ShapeInt a -> Bool
decompose a =
   let (u,sm,vt) = Singular.decompose a
       s = ArrMatrix.toVector sm
       mn = Shape.size $ Array.shape s
   in approxArrayTol 1e-3 a
        (Matrix.takeColumns mn (Matrix.generalizeWide u) #*#
         Matrix.scaleRowsReal
            (Array.reshape (shapeInt mn) s)
            (Matrix.takeRows mn (Matrix.generalizeTall vt)))
      &&
      isUnitary 1e-3 u
      &&
      isUnitary 1e-3 vt

decomposeTall ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.Tall ShapeInt ShapeInt a -> Bool
decomposeTall a =
   let (u,s,vt) = Singular.decomposeTall a
   in approxArrayTol 1e-3 a (u ##*# Matrix.scaleRowsReal s vt)
      &&
      isUnitary 1e-3 u
      &&
      isUnitary 1e-3 vt

decomposeWide ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.Wide ShapeInt ShapeInt a -> Bool
decomposeWide a =
   let (u,s,vt) = Singular.decomposeWide a
   in approxArrayTol 1e-3 a (u #*## Matrix.scaleRowsReal s vt)
      &&
      isUnitary 1e-3 u
      &&
      isUnitary 1e-3 (Matrix.transpose vt)

decomposePolar ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.Tall ShapeInt ShapeInt a -> Bool
decomposePolar a =
   let (u,h) = Singular.decomposePolar a
   in approxArrayTol 1e-3 a (u ##*# h)
      &&
      isUnitary 1e-3 u



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)

testsVar ::
   (Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>
   [(String, Tagged a QC.Property)]
testsVar =
   ("pseudoInverseOrtho",
      checkForAll Gen.matrix pseudoInverseOrtho) :
   ("pseudoInverseProjection",
      checkForAll Gen.matrix pseudoInverseProjection) :
   ("pseudoInverseHermitian",
      checkForAll Gen.matrix pseudoInverseHermitian) :
   ("determinantAbsolute",
      checkForAll Gen.matrix determinantAbsolute) :
   ("leastSquaresMinimumNorm",
      checkForAll ((,) <$> Gen.matrix <#\#> Gen.matrix) leastSquaresMinimumNorm) :
   ("decompose",
      checkForAll Gen.matrix decompose) :
   ("decomposeTall",
      checkForAll Gen.tall decomposeTall) :
   ("decomposeWide",
      checkForAll Gen.wide decomposeWide) :
   ("decomposePolar",
      checkForAll Gen.tall decomposePolar) :
   []