{-# LANGUAGE RebindableSyntax #-}
{-
We represent a Gaussian bell curve in terms of the reciprocal of its variance
and its value at the origin.

We could do some projective geometry in the exponent
in order to also have zero variance,
which corresponds to the dirac impulse.
-}
module MathObj.Gaussian.Variance where

import qualified MathObj.Polynomial as Poly
import qualified Number.Root as Root

import qualified Algebra.Transcendental as Trans
import qualified Algebra.Algebraic      as Algebraic
import qualified Algebra.Field          as Field
import qualified Algebra.Absolute       as Absolute
import qualified Algebra.Ring           as Ring
import qualified Algebra.Additive       as Additive

{-
import Algebra.Transcendental (pi, )
import Algebra.Ring ((*), (^), )
import Algebra.Additive ((+))
-}
import Test.QuickCheck (Arbitrary, arbitrary, )
import Control.Monad (liftM2, )

-- import Prelude (($))
import NumericPrelude.Numeric
import NumericPrelude.Base


{- |
Since @amp@ is the square of the actual amplitude it must be non-negative.
-}
data T a = Cons {amp, c :: a}
   deriving (Eq, Show)

instance (Absolute.C a, Arbitrary a) => Arbitrary (T a) where
   arbitrary =
      liftM2 Cons
         (fmap abs arbitrary)
         (fmap ((1+) . abs) arbitrary)


constant :: Ring.C a => T a
constant = Cons one zero

{-# INLINE evaluate #-}
evaluate :: (Trans.C a) =>
   T a -> a -> a
evaluate f x =
   sqrt (amp f) * exp (-pi * c f * x^2)

exponentPolynomial :: (Additive.C a) =>
   T a -> Poly.T a
exponentPolynomial f =
   Poly.fromCoeffs [zero, zero, c f]


norm1Root :: (Field.C a) => T a -> Root.T a
norm1Root f =
   Root.sqrt $ Root.fromNumber $ amp f / c f

norm2Root :: (Field.C a) => T a -> Root.T a
norm2Root f =
   Root.sqrt $
      Root.fromNumber (amp f)
      `Root.div`
      Root.sqrt (Root.fromNumber $ 2 * c f)

normInfRoot :: (Field.C a) => T a -> Root.T a
normInfRoot f =
   Root.sqrt $ Root.fromNumber $ amp f

normPRoot :: (Field.C a) => Rational -> T a -> Root.T a
normPRoot p f =
   Root.sqrt (Root.fromNumber (amp f))
   `Root.div`
   Root.rationalPower (recip (2*p)) (Root.fromNumber (fromRational' p * c f))


norm1 :: (Algebraic.C a) => T a -> a
norm1 f =
   sqrt $ amp f / c f

norm2 :: (Algebraic.C a) => T a -> a
norm2 f =
   sqrt $ amp f / (sqrt $ 2 * c f)

normInf :: (Algebraic.C a) => T a -> a
normInf f =
   sqrt (amp f)

normP :: (Trans.C a) => a -> T a -> a
normP p f =
   sqrt (amp f) * (p * c f) ^? (- recip (2*p))


variance :: (Trans.C a) =>
   T a -> a
variance f =
   recip $ c f * 2*pi

multiply :: (Ring.C a) =>
   T a -> T a -> T a
multiply f g =
   Cons (amp f * amp g) (c f + c g)

powerRing :: (Trans.C a) =>
   Integer -> T a -> T a
powerRing p f =
   Cons (amp f ^ p) (fromInteger p * c f)

{-
powerField does not makes sense,
since the reciprocal of a Gaussian diverges.
-}

powerAlgebraic :: (Trans.C a) =>
   Rational -> T a -> T a
powerAlgebraic p f =
   Cons (amp f ^/ p) (fromRational' p * c f)

powerTranscendental :: (Trans.C a) =>
   a -> T a -> T a
powerTranscendental p f =
   Cons (amp f ^? p) (p * c f)

{- |
> convolve x y t =
>    integrate $ \s -> x s * y(t-s)

Convergence only for @c f + c g > 0@.
-}
convolve :: (Field.C a) =>
   T a -> T a -> T a
convolve f g =
   let s = c f + c g
   in  Cons
          (amp f * amp g / s)
          (c f * c g / s)

{- |
> fourier x f =
>    integrate $ \t -> x t * cis (-2*pi*t*f)

Convergence only for @c f > 0@.
-}
fourier :: (Field.C a) =>
   T a -> T a
fourier f =
   Cons (amp f / c f) (recip $ c f)
{-
fourier (t -> exp(-(a*t)^2))
-}

dilate :: (Field.C a) => a -> T a -> T a
dilate k f =
   Cons (amp f) $ c f / k^2

shrink :: (Ring.C a) => a -> T a -> T a
shrink k f =
   Cons (amp f) $ c f * k^2

{- |
@amplify k@ scales by @abs k@!
-}
amplify :: (Ring.C a) => a -> T a -> T a
amplify k f =
   Cons (k^2 * amp f) $ c f


{- laws
fourier (convolve f g) = multiply (fourier f) (fourier g)

dilate k (dilate m f) = dilate (k*m) f

dilate k (shrink k f) = f

variance (dilate k f) = k^2 * variance f

variance (convolve f g) = variance f + variance g
-}