{-# LANGUAGE BangPatterns #-}
module QuantLib.Stochastic.Random
        ( BoxMuller
--        , createNormalGen
        , mkNormalGen
        , NormalGenerator (..)
        , InverseNormal
        , mkInverseNormal
        ) where

import           QuantLib.Math.InverseNormal
import           QuantLib.Stochastic.PureMT

class RandomGenerator a where
  create :: IO a
  next   :: a -> (Double, a)
  splitWithSeed :: Integer -> a -> (a, a)

instance RandomGenerator PureMT where
  create :: IO PureMT
create = IO PureMT
newPureMT
  next :: PureMT -> (Double, PureMT)
next   = PureMT -> (Double, PureMT)
randomDouble
  splitWithSeed :: Integer -> PureMT -> (PureMT, PureMT)
splitWithSeed = Integer -> PureMT -> (PureMT, PureMT)
splitMTwithSeed

-- | Box-Muller method
data BoxMuller a = BoxMuller {
        BoxMuller a -> Bool
bmFirst       :: Bool,
        BoxMuller a -> Double
bmSecondValue :: Double,
        BoxMuller a -> a
bmRng         :: a
        }

mkNormalGen :: RandomGenerator a => IO (BoxMuller a)
mkNormalGen :: IO (BoxMuller a)
mkNormalGen = do
        a
rng <- IO a
forall a. RandomGenerator a => IO a
create
        BoxMuller a -> IO (BoxMuller a)
forall (m :: * -> *) a. Monad m => a -> m a
return (BoxMuller a -> IO (BoxMuller a))
-> BoxMuller a -> IO (BoxMuller a)
forall a b. (a -> b) -> a -> b
$! a -> BoxMuller a
forall a. RandomGenerator a => a -> BoxMuller a
createNormalGen a
rng

-- | Creates normally distributed generator
createNormalGen :: RandomGenerator a => a -> BoxMuller a
createNormalGen :: a -> BoxMuller a
createNormalGen a
r = BoxMuller :: forall a. Bool -> Double -> a -> BoxMuller a
BoxMuller {
        bmFirst :: Bool
bmFirst         = Bool
True,
        bmSecondValue :: Double
bmSecondValue   = Double
0.0,
        bmRng :: a
bmRng           = a
r
        }

-- | Normally distributed generator
class NormalGenerator a where
        ngGetNext :: a -> (Double, a)
        ngMkNew   :: a -> IO a
        ngSplit   :: a -> (a, a)
        ngSplit   = Integer -> a -> (a, a)
forall a. NormalGenerator a => Integer -> a -> (a, a)
ngSplitWithSeed Integer
1
        ngSplitWithSeed :: Integer -> a -> (a, a)

instance RandomGenerator a => NormalGenerator (BoxMuller a) where
        ngMkNew :: BoxMuller a -> IO (BoxMuller a)
ngMkNew BoxMuller a
_       = IO (BoxMuller a)
forall a. RandomGenerator a => IO (BoxMuller a)
mkNormalGen
        ngGetNext :: BoxMuller a -> (Double, BoxMuller a)
ngGetNext = BoxMuller a -> (Double, BoxMuller a)
forall a. RandomGenerator a => BoxMuller a -> (Double, BoxMuller a)
boxMullerGetNext
        ngSplitWithSeed :: Integer -> BoxMuller a -> (BoxMuller a, BoxMuller a)
ngSplitWithSeed Integer
seed  BoxMuller a
x   = (BoxMuller a
x { bmRng :: a
bmRng = a
rng1 }, BoxMuller a
x { bmRng :: a
bmRng = a
rng2 })
          where
              (a
rng1, a
rng2) = Integer -> a -> (a, a)
forall a. RandomGenerator a => Integer -> a -> (a, a)
splitWithSeed Integer
seed (BoxMuller a -> a
forall a. BoxMuller a -> a
bmRng BoxMuller a
x)

boxMullerGetNext :: RandomGenerator a => BoxMuller a -> (Double, BoxMuller a)
boxMullerGetNext :: BoxMuller a -> (Double, BoxMuller a)
boxMullerGetNext (BoxMuller Bool
True Double
_ a
rng) = (Double
s1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ratio, BoxMuller :: forall a. Bool -> Double -> a -> BoxMuller a
BoxMuller {
                                            bmFirst :: Bool
bmFirst         = Bool
False,
                                            bmSecondValue :: Double
bmSecondValue   = Double
s2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ratio,
                                            bmRng :: a
bmRng           = a
g2
                                          })
        where
          (Double
x1, a
g1) = a -> (Double, a)
forall a. RandomGenerator a => a -> (Double, a)
next a
rng
          (Double
x2, a
g2) = a -> (Double, a)
forall a. RandomGenerator a => a -> (Double, a)
next a
g1
          (Double
r, Double
s1, Double
s2) = (Double, Double, Double)
getRs
          ratio :: Double
ratio = Double -> Double
boxMullerRatio Double
r
          getRs :: (Double, Double, Double)
getRs =
            let
              as1 :: Double
as1 = Double
2.0Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1.0
              as2 :: Double
as2 = Double
2.0Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1.0
              ar :: Double
ar = Double
s1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
s1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
s2
            in
              if Double
rDouble -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>=Double
1.0 Bool -> Bool -> Bool
|| Double
rDouble -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<=Double
0.0 then (Double, Double, Double)
getRs else (Double
ar, Double
as1, Double
as2)
boxMullerGetNext (BoxMuller Bool
False !Double
s !a
r) = (Double
s, Bool -> Double -> a -> BoxMuller a
forall a. Bool -> Double -> a -> BoxMuller a
BoxMuller Bool
True Double
s a
r)

{-# ANN boxMullerRatio "NoHerbie" #-}
boxMullerRatio :: Double -> Double
boxMullerRatio :: Double -> Double
boxMullerRatio Double
r = Double -> Double
forall a. Floating a => a -> a
sqrt (-Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
r)

-- | Normal number generation using inverse cummulative normal distribution
newtype InverseNormal a = InverseNormal a

mkInverseNormal ::  RandomGenerator a => IO (InverseNormal a)
mkInverseNormal :: IO (InverseNormal a)
mkInverseNormal = do
        a
rng <- IO a
forall a. RandomGenerator a => IO a
create
        InverseNormal a -> IO (InverseNormal a)
forall (m :: * -> *) a. Monad m => a -> m a
return (InverseNormal a -> IO (InverseNormal a))
-> InverseNormal a -> IO (InverseNormal a)
forall a b. (a -> b) -> a -> b
$! a -> InverseNormal a
forall a. a -> InverseNormal a
InverseNormal a
rng

instance RandomGenerator a => NormalGenerator (InverseNormal a) where
        ngMkNew :: InverseNormal a -> IO (InverseNormal a)
ngMkNew InverseNormal a
_       = IO (InverseNormal a)
forall a. RandomGenerator a => IO (InverseNormal a)
mkInverseNormal
        ngGetNext :: InverseNormal a -> (Double, InverseNormal a)
ngGetNext (InverseNormal a
rng)   = (Double -> Double
inverseNormal Double
x, a -> InverseNormal a
forall a. a -> InverseNormal a
InverseNormal a
newRng)
          where (Double
x, a
newRng) = a -> (Double, a)
forall a. RandomGenerator a => a -> (Double, a)
next a
rng
        ngSplitWithSeed :: Integer -> InverseNormal a -> (InverseNormal a, InverseNormal a)
ngSplitWithSeed Integer
seed (InverseNormal a
x) = (a -> InverseNormal a
forall a. a -> InverseNormal a
InverseNormal a
x1, a -> InverseNormal a
forall a. a -> InverseNormal a
InverseNormal a
x2)
          where
            (a
x1, a
x2) = Integer -> a -> (a, a)
forall a. RandomGenerator a => Integer -> a -> (a, a)
splitWithSeed Integer
seed a
x