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

import qualified MathObj.Gaussian.Variance as G
import qualified Number.Root as Root

-- import qualified Algebra.Ring           as Ring

import qualified Algebra.Laws as Laws

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

import Control.Monad (liftM2, liftM3, )

import Data.Function.HT (nest, compose2, )

import NumericPrelude.Base as P
import NumericPrelude.Numeric as NP


newtype PositiveInteger = PositiveInteger Integer
   deriving Show

instance Arbitrary PositiveInteger where
   arbitrary =
      fmap (\p -> PositiveInteger $ 1 + abs p) arbitrary


{- |
For @(HoelderConjugates p q)@ it holds

> 1/p + 1/q = 1
-}
data HoelderConjugates = HoelderConjugates Rational Rational
   deriving Show

{-
instance Arbitrary HoelderConjugates where
   arbitrary = liftM2
      (\(PositiveInteger p) (PositiveInteger q) ->
         let s  = 1%p + 1%q
         in  HoelderConjugates (fromInteger p * s) (fromInteger q * s))
      arbitrary arbitrary
-}
instance Arbitrary HoelderConjugates where
   arbitrary = liftM2
      (\(PositiveInteger p) (PositiveInteger q) ->
         let s = p + q
         in  HoelderConjugates (s % p) (s % q))
      arbitrary arbitrary

{- |
For @(YoungConjugates p q r)@ it holds

> 1/p + 1/q = 1/r + 1
-}
data YoungConjugates = YoungConjugates Rational Rational Rational
   deriving Show

{-
Find positive natural numbers @a, b, c, d@ with

> a + b = c + d

and

> d >= a, d >= b, d >= c

then set

> p=d/a, q=d/b, r=d/c


a+b<=c
b+c<=a
->  2b <= 0
-}
instance Arbitrary YoungConjugates where
   arbitrary = liftM3
      (\(PositiveInteger a0) (PositiveInteger b0) (PositiveInteger c0) ->
         let guardSwap cond (x,y) =
                if cond x y then (x,y) else (y,x)
             {-
             If a+b<=c, then from b>0 it follows a<c and thus c+b>a.
             Swapping a and c is enough and we have not to consider more cases.
             -}
             (a1,c1) = guardSwap (\a c -> a+b0>c) (a0,c0)
             b1 = b0
             d1 = a1+b1-c1
             ((a2,b2),(c2,d2)) =
                guardSwap (compose2 (<=) snd)
                   (guardSwap (<=) (a1,b1),
                    guardSwap (<=) (c1,d1))
         in  YoungConjugates (d2%a2) (d2%b2) (d2%c2))
      arbitrary arbitrary arbitrary

{-
This is simpler, but may yield exponents smaller than 1.

instance Arbitrary YoungConjugates where
   arbitrary = liftM3
      (\(PositiveInteger a0) (PositiveInteger b0) (PositiveInteger c0) ->
         let {-
             If a+b<=c, then from b>0 it follows a<c and thus c+b>a.
             Swapping a and c is enough and we have not to consider more cases.
             -}
             (a1,c1) = if a0+b0<=c0 then (c0,a0) else (a0,c0)
             b1 = b0
             d1 = a1+b1-c1
         in  YoungConjugates (d1%a1) (d1%b1) (d1%c1))
      arbitrary arbitrary arbitrary
-}


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

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

testList :: [(String, IO ())]
testList =
{-
      ("convolution, dirac",
          simple $ Laws.identity (+) zero) :
-}
      ("convolution, commutative",
          simple $ Laws.commutative G.convolve) :
      ("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 via fourier",
          simple $ \x y ->
             G.fourier (G.convolve x y)
              == G.multiply (G.fourier x) (G.fourier y)) :
      ("fourier identity",
          simple $ \x -> nest 4 G.fourier x == 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))) :
      ("fourier, unitary",
          simple $ \x y ->
             G.scalarProductRoot x y
              == G.scalarProductRoot (G.fourier x) (G.fourier y)) :
      ("norm1 vs. normP 1",
          simple $ \x -> G.norm1Root x == G.normPRoot 1 x) :
      ("norm2 vs. normP 2",
          simple $ \x -> G.norm2Root x == G.normPRoot 2 x) :
{-
I would have liked to test for a monotony of norms.
Unfortunately, it does not hold.

Means contain a division by the size of the domain.
Norms do not have this division.
Means are monotonic with respect to the degree.
Norms are not.
We cannot turn the norms into means since the size of the domain
(the complete real axis) is infinitely large.
      ("norm monotony",
          simple $ \x p0 q0 ->
             let p = 1 + abs p0
                 q = 1 + abs q0
             in  case compare p q of
                    EQ -> G.normPRoot p x == G.normPRoot q x
                    LT -> G.normPRoot p x <= G.normPRoot q x
                    GT -> G.normPRoot p x >= G.normPRoot q x) :

This should also fail,
but QuickCheck does not seem to try counterexamples.
      ("infinity norm upper bound",
          simple $ \x p0 ->
             let p = 1 + abs p0
             in  G.normPRoot p x <= G.normInfRoot x) :
-}
      ("Cauchy-Schwarz inequality",
          simple $ \x y ->
             G.scalarProductRoot x y
                <= G.norm2Root x `Root.mul` G.norm2Root y) :
      ("Hoelder conjugates",
          quickCheck $ \(HoelderConjugates p q) ->
             p>=1 && q>=1 && 1/p + 1/q == 1) :
      ("Hoelder inequality with infinity norm",
          simple $ \x y ->
             G.scalarProductRoot x y
                <= G.norm1Root x `Root.mul` G.normInfRoot y) :
      ("Hoelder inequality",
          simple $ \x y (HoelderConjugates p q) ->
             G.scalarProductRoot x y
                <= G.normPRoot p x `Root.mul` G.normPRoot q y) :
      ("Young inequality with two infinity norms",
          simple $ \x y ->
             G.normInfRoot (G.convolve x y)
                <= G.norm1Root x `Root.mul` G.normInfRoot y) :
      ("Young inequality with infinity norm",
          simple $ \x y (HoelderConjugates p q) ->
             G.normInfRoot (G.convolve x y)
                <= G.normPRoot p x `Root.mul` G.normPRoot q y) :
      ("Young conjugates",
          quickCheck $ \(YoungConjugates p q r) ->
             p>=1 && q>=1 && r>=1 && 1/p + 1/q == 1/r + 1) :
      ("Young inequality",
          simple $ \x y (YoungConjugates p q r) ->
             G.normPRoot r (G.convolve x y)
                <= G.normPRoot p x `Root.mul` G.normPRoot q y) :
      []