-----------------------------------------------------------------------------
-- |
-- Module  :  ForSyDe.Shallow.Utility.Gaussian
-- Copyright   :  (c) ForSyDe Group, KTH 2007-2008
-- License     :  BSD-style (see the file LICENSE)
-- 
-- Maintainer  :  forsyde-dev@ict.kth.se
-- Stability   :  experimental
-- Portability :  portable
--
-- We follow the Box-Muller method to generate white gaussian noise, 
-- described at: <http://www.dspguru.com/howto/tech/wgn.htm>
-----------------------------------------------------------------------------
module ForSyDe.Shallow.Utility.Gaussian (
      pGaussianNoise
    )
where

import ForSyDe.Shallow.MoC.Untimed
import ForSyDe.Shallow.Core.Signal

import qualified System.Random

-- |To generate an infinite Signal of Gaussian values
pGaussianNoise:: Double -- Mean value of the Gaussian noise
    -> Double   -- Variance of the Gaussian noise
    -> Int  -- The seed
    -> Signal Double -- Output gaussian noise signal
pGaussianNoise :: Double -> Double -> Int -> Signal Double
pGaussianNoise Double
mean Double
variance = Int -> ([Double] -> [Double]) -> Signal Double -> Signal Double
forall a b. Int -> ([a] -> [b]) -> Signal a -> Signal b
mapU Int
2 [Double] -> [Double]
gaussianXY (Signal Double -> Signal Double)
-> (Int -> Signal Double) -> Int -> Signal Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Signal Double
pUnitNormXY
  where
    gaussianXY :: [Double] -> [Double]
gaussianXY [Double
x, Double
y] = [Double
mean Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt(Double
variance) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x,
         Double
mean Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt(Double
variance) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y]
    gaussianXY [Double]
_  = [Char] -> [Double]
forall a. HasCallStack => [Char] -> a
error [Char]
"gaussianXY: unexpected pattern."

-- |To get a uniform random variable in the range [0, 1]
uniform :: (Fractional a, System.Random.RandomGen g, System.Random.Random a) => 
  g -> (a, g)
uniform :: g -> (a, g)
uniform g
rGen = (a, a) -> g -> (a, g)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
System.Random.randomR (a
0.0,a
1.0) g
rGen

-- |To generate an infinite signal of unit normal random variables,
-- with the specified seed
pUnitNormXY :: Int     -- The seed
     -> Signal Double  -- The infinite ouput signal
pUnitNormXY :: Int -> Signal Double
pUnitNormXY = Int -> ([Double] -> [Double]) -> Signal Double -> Signal Double
forall a b. Int -> ([a] -> [b]) -> Signal a -> Signal b
mapU Int
3 [Double] -> [Double]
forall a. Floating a => [a] -> [a]
unitNormXY (Signal Double -> Signal Double)
-> (Int -> Signal Double) -> Int -> Signal Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Signal Double
forall a. [a] -> Signal a
signal ([Double] -> Signal Double)
-> (Int -> [Double]) -> Int -> Signal Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. StdGen -> [Double]
svGenerator (StdGen -> [Double]) -> (Int -> StdGen) -> Int -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> StdGen
System.Random.mkStdGen
  where
    unitNormXY :: [a] -> [a]
unitNormXY [a
s, a
v1, a
v2] = [a -> a
forall a. Floating a => a -> a
sqrt(-a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
log(a
s) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s) a -> a -> a
forall a. Num a => a -> a -> a
* a
v1,
          a -> a
forall a. Floating a => a -> a
sqrt(-a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
log(a
s) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s) a -> a -> a
forall a. Num a => a -> a -> a
* a
v2]
    unitNormXY [a]
_ = [Char] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"pUnitNormXY: Unexpected pattern."


-- |To generate the s, v1, v2 value
svGenerator :: System.Random.StdGen -> [Double]
svGenerator :: StdGen -> [Double]
svGenerator StdGen
s
    | Double
sVal Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>=Double
1 = [][Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ StdGen -> [Double]
svGenerator StdGen
newStdG
    | Bool
otherwise = [Double]
svVal [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ StdGen -> [Double]
svGenerator StdGen
newStdG
  where
    svGen1 :: ([Double], StdGen)
svGen1 = StdGen -> ([Double], StdGen)
svHelper  StdGen
s
    svVal :: [Double]
svVal = ([Double], StdGen) -> [Double]
forall a b. (a, b) -> a
fst ([Double], StdGen)
svGen1
    sVal :: Double
sVal = [Double] -> Double
forall a. [a] -> a
head [Double]
svVal
    newStdG :: StdGen
newStdG = ([Double], StdGen) -> StdGen
forall a b. (a, b) -> b
snd ([Double], StdGen)
svGen1
    svHelper :: System.Random.StdGen -> ([Double], System.Random.StdGen)
    svHelper :: StdGen -> ([Double], StdGen)
svHelper StdGen
stdG = ([Double
s, Double
v1, Double
v2], StdGen
sNew2)
      where
        (Double
u1, StdGen
sNew1) = StdGen -> (Double, StdGen)
forall a g. (Fractional a, RandomGen g, Random a) => g -> (a, g)
uniform StdGen
stdG
        (Double
u2, StdGen
sNew2) = StdGen -> (Double, StdGen)
forall a g. (Fractional a, RandomGen g, Random a) => g -> (a, g)
uniform StdGen
sNew1
        v1 :: Double
v1=Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u1 Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1
        v2 :: Double
v2=Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u2 Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1
        s :: Double
s = Double
v1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
v1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
v2