{-# LANGUAGE NoImplicitPrelude #-}
{-
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.

The Gaussians form a nice multiplicative commutative monoid.
Maybe we should have such a structure.
It would also be useful for the Root data type
and a new Exponential data type.
-}
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

{- |
eigenfunction of 'fourier'
-}
unit :: Ring.C a => T a
unit = Cons one one

{-# 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]


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

scalarProductRoot :: (Field.C a) => T a -> T a -> Root.T a
scalarProductRoot f g =
   integrateRoot (multiply f g)


norm1Root :: (Field.C a) => T a -> Root.T a
norm1Root = integrateRoot

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))


-- ToDo: implement NormedSpace.Sum et.al.
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
-}