{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
module Test.MathObj.Gaussian.Polynomial where

import qualified MathObj.Gaussian.Polynomial as G
import qualified MathObj.Gaussian.Bell as B

import qualified MathObj.Polynomial as Poly

-- import qualified Algebra.Ring           as Ring

import qualified Algebra.Laws as Laws

import qualified Number.Complex as Complex

import Test.NumericPrelude.Utility (testUnit)
import Test.QuickCheck (Testable, quickCheck, (==>))
import qualified Test.HUnit as HUnit

import qualified Number.NonNegative as NonNeg
import Data.Function.HT (nest, )
import Data.Tuple.HT (mapSnd, )

-- import Debug.Trace (trace, )

import PreludeBase as P
import NumericPrelude as NP


simple ::
   (Testable t) =>
   (G.T Rational -> t) -> IO ()
simple f =
   quickCheck (\x -> f (x :: G.T Rational))

tests :: HUnit.Test
tests =
   HUnit.TestLabel "polynomial" $
   HUnit.TestList $
   map testUnit $
   testList

testList :: [(String, IO ())]
testList =
{-
      ("convolution, dirac",
          simple $ Laws.identity (+) zero) :
-}
      ("convolution, commutative",
          simple $ Laws.commutative G.convolve) :
--          simple $ \x -> Laws.commutative G.convolve (trace (show x) x)) :
      ("convolution, associative",
          simple $ Laws.associative G.convolve) :
      ("multiplication, one",
          simple $ Laws.identity G.multiply G.constant) :
      ("multiplication, commutative",
          simple $ Laws.commutative G.multiply) :
      ("multiplication, associative",
          simple $ Laws.associative G.multiply) :
      ("convolution, multplication, fourier",
          simple $ \x y ->
             G.fourier (G.convolve x y)
              == G.multiply (G.fourier x) (G.fourier y)) :
      ("fourier reverse",
          simple $ \x -> nest 2 G.fourier x == G.reverse x) :
      ("reverse identity",
          simple $ \x -> nest 2 G.reverse x == x) :
      ("fourier eigenfunction differential",
          quickCheck $ \m ->
             m <= 15 ==>
                let n = NonNeg.toNumber m
                    x = G.eigenfunctionDifferential n :: G.T Rational
                    k = Complex.conjugate Complex.imaginaryUnit ^ fromIntegral n
                in  G.fourier x  ==  G.scaleComplex k x) :
      ("fourier eigenfunction iterative",
          quickCheck $ \m ->
             m <= 15 ==>
                let n = NonNeg.toNumber m
                    x = G.eigenfunctionIterative n :: G.T Rational
                    k = Complex.conjugate Complex.imaginaryUnit ^ fromIntegral n
                in  G.fourier x  ==  G.scaleComplex k x) :
{- this does not hold, both functions compute different eigenbases
      ("fourier eigenfunction diff vs. iterative",
          quickCheck $ \n ->
             n <= 15 ==>
                G.eigenfunctionDifferential n ==
                (G.eigenfunctionIterative n :: G.T Rational)) :
-}
      ("translate additive",
          simple $ \x a b ->
             G.translate a (G.translate b x) == G.translate (a+b) x) :
      ("translateComplex additive",
          simple $ \x a b ->
             G.translateComplex a (G.translateComplex b x) == G.translateComplex (a+b) x) :
      ("translateComplex real",
          simple $ \x a ->
             G.translateComplex (Complex.fromReal a) x == G.translate a x) :
      ("modulate additive",
          simple $ \x a b ->
             G.modulate a (G.modulate b x) == G.modulate (a+b) x) :
      ("commute translate modulate",
          simple $ \x a b ->
             G.modulate b (G.translate a x)
              == G.turn (a*b) (G.translate a (G.modulate b x))) :
      ("fourier translate",
          simple $ \x a ->
             G.fourier (G.translate a x)
              == G.modulate a (G.fourier x)) :
      ("dilate multiplicative",
          simple $ \x a b -> a>0 && b>0 ==>
             G.dilate a (G.dilate b x) == G.dilate (a*b) x) :
      ("dilate by shrink",
          simple $ \x a -> a>0 ==>
             G.shrink a x == G.dilate (recip a) x) :
      ("fourier dilate",
          simple $ \x a -> a>0 ==>
             G.fourier (G.dilate a x) == G.amplify a (G.shrink a (G.fourier x))) :
      ("integrate differentiate",
          simple $ \x ->
             G.integrate (G.differentiate x) == (zero, x)) :
      ("fourier differentiate",
          simple $ \x ->
             G.fourier (G.differentiate x) ==
              let y = G.fourier x
              in  y{G.polynomial =
                      Poly.fromCoeffs [0, 0 Complex.+: 2] * G.polynomial y}) :
      ("approximate by bells, translate",
          simple $ \x unit d -> unit/=0 ==>
             G.approximateByBells unit (G.translateComplex d x) ==
             map (mapSnd (B.translateComplex d)) (G.approximateByBells unit x)) :
      ("approximate by bells, dilate",
          simple $ \x unit d -> unit/=0 && d/=0 ==>
             G.approximateByBells unit (G.dilate d x) ==
             map (mapSnd (B.dilate d)) (G.approximateByBells (unit/d) x)) :
      ("approximate by bells, shrink",
          simple $ \x unit d -> unit/=0 && d/=0 ==>
             G.approximateByBells unit (G.shrink d x) ==
             map (mapSnd (B.shrink d)) (G.approximateByBells (unit*d) x)) :
      ("approximate by bells, different implementations",
          quickCheck $ \unit d s p -> unit/=0 ==>
             G.approximateByBellsAtOnce unit d s (p::Poly.T (Complex.T Rational)) ==
             G.approximateByBellsByTranslation unit d s p) :
      []