{-# LANGUAGE TypeFamilies #-}
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 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,s,vt) = Singular.decompose a
       mn = Shape.size $ Array.shape s
   in approxArrayTol 1e-3 a
        (Matrix.takeColumns mn (Matrix.generalizeWide u) #*#
         Matrix.scaleRowsReal
            (Array.mapShape (shapeInt . Shape.size) 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) :
   []