{-# 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.
-}
module MathObj.Gaussian.Variance where

import qualified MathObj.Polynomial as Poly

import qualified Algebra.Transcendental as Trans
import qualified Algebra.Algebraic      as Algebraic
import qualified Algebra.Field          as Field
import qualified Algebra.Real           as Real
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
import PreludeBase


data T a = Cons {amp, c :: a}
   deriving (Eq, Show)

instance (Real.C a, Arbitrary a) => Arbitrary (T a) where
   arbitrary =
      liftM2 Cons
         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]


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

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

normP :: (Trans.C a, Real.C a) => a -> T a -> a
normP p f =
   sqrt (abs (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)

{- |
> convolve x y t =
>    integrate $ \s -> x s * y(t-s)
-}
convolve :: (Field.C a) =>
   T a -> T a -> T a
convolve f g =
   Cons
      (amp f * amp g / (c f + c g))
      (recip $ recip (c f) + recip (c g))

{- |
> fourier x f =
>    integrate $ \t -> x t * cis (-2*pi*t*f)
-}
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 :: (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
-}