{-# LANGUAGE RankNTypes #-}

-- |
-- Module     : Simulation.Aivika.Generator
-- Copyright  : Copyright (c) 2009-2017, David Sorokin <david.sorokin@gmail.com>
-- License    : BSD3
-- Maintainer : David Sorokin <david.sorokin@gmail.com>
-- Stability  : experimental
-- Tested with: GHC 8.0.1
--
-- Below is defined a type class of the random number generator.
--
module Simulation.Aivika.Generator 
       (Generator(..),
        GeneratorType(..),
        DiscretePDF(..),
        newGenerator,
        newRandomGenerator,
        newRandomGenerator01) where

import System.Random
import qualified System.Random.MWC as MWC
import Data.IORef
import Data.Word
import Data.Vector
import Data.Functor

-- | A discrete probability density function.
type DiscretePDF a = [(a, Double)]

-- | Defines a random number generator.
data Generator =
  Generator { Generator -> Double -> Double -> IO Double
generateUniform :: Double -> Double -> IO Double,
              -- ^ Generate an uniform random number
              -- with the specified minimum and maximum.
              Generator -> Int -> Int -> IO Int
generateUniformInt :: Int -> Int -> IO Int,
              -- ^ Generate an uniform integer random number
              -- with the specified minimum and maximum.
              Generator -> Double -> Double -> Double -> IO Double
generateTriangular :: Double -> Double -> Double -> IO Double,
              -- ^ Generate a triangular random number
              -- by the specified minimum, median and maximum.
              Generator -> Double -> Double -> IO Double
generateNormal :: Double -> Double -> IO Double,
              -- ^ Generate the normal random number
              -- with the specified mean and deviation.
              Generator -> Double -> Double -> IO Double
generateLogNormal :: Double -> Double -> IO Double,
              -- ^ Generate a random number from the lognormal distribution derived
              -- from a normal distribution with the specified mean and deviation.
              Generator -> Double -> IO Double
generateExponential :: Double -> IO Double,
              -- ^ Generate the random number distributed exponentially
              -- with the specified mean (the reciprocal of the rate).
              Generator -> Double -> Int -> IO Double
generateErlang :: Double -> Int -> IO Double,
              -- ^ Generate the Erlang random number
              -- with the specified scale (the reciprocal of the rate) and integer shape.
              Generator -> Double -> IO Int
generatePoisson :: Double -> IO Int,
              -- ^ Generate the Poisson random number
              -- with the specified mean.
              Generator -> Double -> Int -> IO Int
generateBinomial :: Double -> Int -> IO Int,
              -- ^ Generate the binomial random number
              -- with the specified probability and number of trials.
              Generator -> Double -> Double -> IO Double
generateGamma :: Double -> Double -> IO Double,
              -- ^ Generate a random number from the Gamma distribution with
              -- the specified shape (kappa) and scale (theta, a reciprocal of the rate).
              --
              -- The probability density for the Gamma distribution is
              --
              -- @f x = x ** (kappa - 1) * exp (- x \/ theta) \/ theta ** kappa * Gamma kappa@
              Generator -> Double -> Double -> IO Double
generateBeta :: Double -> Double -> IO Double,
              -- ^ Generate a random number from the Beta distribution by
              -- the specified shape parameters (alpha and beta).
              --
              -- The probability density for the Beta distribution is
              --
              -- @f x = x ** (alpha - 1) * (1 - x) ** (beta - 1) \/ B alpha beta@
              Generator -> Double -> Double -> IO Double
generateWeibull :: Double -> Double -> IO Double,
              -- ^ Generate a random number from the Weibull distribution by
              -- the specified shape and scale.
              Generator -> forall a. DiscretePDF a -> IO a
generateDiscrete :: forall a. DiscretePDF a -> IO a,
              -- ^ Generate a random value from the specified discrete distribution.
              Generator -> IO Int
generateSequenceNo :: IO Int
              -- ^ Generate a sequence number which can be considered quite unique.
            }

-- | Generate the uniform random number with the specified minimum and maximum.
generateUniform01 :: IO Double
                     -- ^ the generator
                     -> Double
                     -- ^ minimum
                     -> Double
                     -- ^ maximum
                     -> IO Double
generateUniform01 :: IO Double -> Double -> Double -> IO Double
generateUniform01 IO Double
g Double
min Double
max =
  do Double
x <- IO Double
g
     Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
min Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
max Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
min)

-- | Generate the uniform random number with the specified minimum and maximum.
generateUniformInt01 :: IO Double
                        -- ^ the generator
                        -> Int
                        -- ^ minimum
                        -> Int
                        -- ^ maximum
                        -> IO Int
generateUniformInt01 :: IO Double -> Int -> Int -> IO Int
generateUniformInt01 IO Double
g Int
min Int
max =
  do Double
x <- IO Double
g
     let min' :: Double
min' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
min Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5
         max' :: Double
max' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
max Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5
         z :: Int
z    = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double
min' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
max' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
min'))
         z' :: Int
z'   = if Int
z Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
min
                then Int
min
                else if Int
z Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
max
                     then Int
max
                     else Int
z
     Int -> IO Int
forall (m :: * -> *) a. Monad m => a -> m a
return Int
z'

-- | Generate the triangular random number by the specified minimum, median and maximum.
generateTriangular01 :: IO Double
                        -- ^ the generator
                        -> Double
                        -- ^ minimum
                        -> Double
                        -- ^ median
                        -> Double
                        -- ^ maximum
                        -> IO Double
generateTriangular01 :: IO Double -> Double -> Double -> Double -> IO Double
generateTriangular01 IO Double
g Double
min Double
median Double
max =
  do Double
x <- IO Double
g
     if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= (Double
median Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
min) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
max Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
min)
       then Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
min Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
median Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
min) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
max Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
min) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
       else Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
max Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
max Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
median) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
max Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
min) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x))

-- | Create a normal random number generator with mean 0 and variance 1
-- by the specified generator of uniform random numbers from 0 to 1.
newNormalGenerator01 :: IO Double
                        -- ^ the generator
                        -> IO (IO Double)
newNormalGenerator01 :: IO Double -> IO (IO Double)
newNormalGenerator01 IO Double
g =
  do IORef Double
nextRef <- Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Bool
flagRef <- Bool -> IO (IORef Bool)
forall a. a -> IO (IORef a)
newIORef Bool
False
     IORef Double
xi1Ref  <- Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Double
xi2Ref  <- Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Double
psiRef  <- Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     let loop :: IO ()
loop =
           do Double
psi <- IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
psiRef
              if (Double
psi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
1.0) Bool -> Bool -> Bool
|| (Double
psi Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0)
                then do Double
g1 <- IO Double
g
                        Double
g2 <- IO Double
g
                        let xi1 :: Double
xi1 = Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0
                            xi2 :: Double
xi2 = Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0
                            psi :: Double
psi = Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
xi2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xi2
                        IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
xi1
                        IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
xi2
                        IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
psi
                        IO ()
loop
                else IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef (Double -> IO ()) -> Double -> IO ()
forall a b. (a -> b) -> a -> b
$ 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
psi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
psi)
     IO Double -> IO (IO Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (IO Double -> IO (IO Double)) -> IO Double -> IO (IO Double)
forall a b. (a -> b) -> a -> b
$
       do Bool
flag <- IORef Bool -> IO Bool
forall a. IORef a -> IO a
readIORef IORef Bool
flagRef
          if Bool
flag
            then do IORef Bool -> Bool -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
False
                    IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
nextRef
            else do IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
0.0
                    IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
0.0
                    IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
0.0
                    IO ()
loop
                    Double
xi1 <- IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
xi1Ref
                    Double
xi2 <- IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
xi2Ref
                    Double
psi <- IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
psiRef
                    IORef Bool -> Bool -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
True
                    IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
nextRef (Double -> IO ()) -> Double -> IO ()
forall a b. (a -> b) -> a -> b
$ Double
xi2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
psi
                    Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
psi

-- | Return the exponential random number with the specified mean.
generateExponential01 :: IO Double
                         -- ^ the generator
                         -> Double
                         -- ^ the mean
                         -> IO Double
generateExponential01 :: IO Double -> Double -> IO Double
generateExponential01 IO Double
g Double
mu =
  do Double
x <- IO Double
g
     Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (- Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
mu)

-- | Return the Erlang random number.
generateErlang01 :: IO Double
                    -- ^ the generator
                    -> Double
                    -- ^ the scale
                    -> Int
                    -- ^ the shape
                    -> IO Double
generateErlang01 :: IO Double -> Double -> Int -> IO Double
generateErlang01 IO Double
g Double
beta Int
m =
  do Double
x <- Int -> Double -> IO Double
forall t. (Ord t, Num t) => t -> Double -> IO Double
loop Int
m Double
1
     Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (- Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
beta)
       where loop :: t -> Double -> IO Double
loop t
m Double
acc
               | t
m t -> t -> Bool
forall a. Ord a => a -> a -> Bool
< t
0     = [Char] -> IO Double
forall a. HasCallStack => [Char] -> a
error [Char]
"Negative shape: generateErlang."
               | t
m t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== t
0    = Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
acc
               | Bool
otherwise = do Double
x <- IO Double
g
                                t -> Double -> IO Double
loop (t
m t -> t -> t
forall a. Num a => a -> a -> a
- t
1) (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
acc)

-- | Generate the Poisson random number with the specified mean.
generatePoisson01 :: IO Double
                     -- ^ the generator
                     -> Double
                     -- ^ the mean
                     -> IO Int
generatePoisson01 :: IO Double -> Double -> IO Int
generatePoisson01 IO Double
g Double
mu =
  do Double
prob0 <- IO Double
g
     let loop :: Double -> Double -> t -> m t
loop Double
prob Double
prod t
acc
           | Double
prob Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
prod = t -> m t
forall (m :: * -> *) a. Monad m => a -> m a
return t
acc
           | Bool
otherwise    = Double -> Double -> t -> m t
loop
                            (Double
prob Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
prod)
                            (Double
prod Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
mu Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ t -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (t
acc t -> t -> t
forall a. Num a => a -> a -> a
+ t
1))
                            (t
acc t -> t -> t
forall a. Num a => a -> a -> a
+ t
1)
     Double -> Double -> Int -> IO Int
forall (m :: * -> *) t.
(Monad m, Integral t) =>
Double -> Double -> t -> m t
loop Double
prob0 (Double -> Double
forall a. Floating a => a -> a
exp (- Double
mu)) Int
0

-- | Generate a binomial random number with the specified probability and number of trials. 
generateBinomial01 :: IO Double
                      -- ^ the generator
                      -> Double 
                      -- ^ the probability
                      -> Int
                      -- ^ the number of trials
                      -> IO Int
generateBinomial01 :: IO Double -> Double -> Int -> IO Int
generateBinomial01 IO Double
g Double
prob Int
trials = Int -> Int -> IO Int
forall a t. (Ord a, Num a, Num t) => a -> t -> IO t
loop Int
trials Int
0 where
  loop :: a -> t -> IO t
loop a
n t
acc
    | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0     = [Char] -> IO t
forall a. HasCallStack => [Char] -> a
error [Char]
"Negative number of trials: generateBinomial."
    | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0    = t -> IO t
forall (m :: * -> *) a. Monad m => a -> m a
return t
acc
    | Bool
otherwise = do Double
x <- IO Double
g
                     if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
prob
                       then a -> t -> IO t
loop (a
n a -> a -> a
forall a. Num a => a -> a -> a
- a
1) (t
acc t -> t -> t
forall a. Num a => a -> a -> a
+ t
1)
                       else a -> t -> IO t
loop (a
n a -> a -> a
forall a. Num a => a -> a -> a
- a
1) t
acc

-- | Generate a random number from the Gamma distribution using Marsaglia and Tsang method.
generateGamma01 :: IO Double
                   -- ^ a normal random number ~ N (0,1)
                   -> IO Double
                   -- ^ an uniform random number ~ U (0, 1)
                   -> Double
                   -- ^ the shape parameter (kappa) 
                   -> Double
                   -- ^ the scale parameter (theta)
                   -> IO Double
generateGamma01 :: IO Double -> IO Double -> Double -> Double -> IO Double
generateGamma01 IO Double
gn IO Double
gu Double
kappa Double
theta
  | Double
kappa Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 = [Char] -> IO Double
forall a. HasCallStack => [Char] -> a
error [Char]
"The shape parameter (kappa) must be positive: generateGamma01"
  | Double
kappa Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1  =
    let d :: Double
d = Double
kappa Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
        c :: Double
c = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d)
        loop :: IO Double
loop =
          do Double
z <- IO Double
gn
             if Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= - (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c)
               then IO Double
loop
               else do let v :: Double
v = (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
3
                       Double
u <- IO Double
gu
                       if Double -> Double
forall a. Floating a => a -> a
log Double
u Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
v
                         then IO Double
loop
                         else Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
theta
    in IO Double
loop
  | Bool
otherwise  =
    do Double
x <- IO Double -> IO Double -> Double -> Double -> IO Double
generateGamma01 IO Double
gn IO Double
gu (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
kappa) Double
theta
       Double
u <- IO Double
gu
       Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
kappa)

-- | Generate a random number from the Beta distribution.
generateBeta01 :: IO Double
                  -- ^ a normal random number ~ N (0, 1)
                  -> IO Double
                  -- ^ an uniform random number ~ U (0, 1)
                  -> Double
                  -- ^ the shape parameter alpha
                  -> Double
                  -- ^ the shape parameter beta
                  -> IO Double
generateBeta01 :: IO Double -> IO Double -> Double -> Double -> IO Double
generateBeta01 IO Double
gn IO Double
gu Double
alpha Double
beta =
  do Double
g1 <- IO Double -> IO Double -> Double -> Double -> IO Double
generateGamma01 IO Double
gn IO Double
gu Double
alpha Double
1
     Double
g2 <- IO Double -> IO Double -> Double -> Double -> IO Double
generateGamma01 IO Double
gn IO Double
gu Double
beta Double
1
     Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
g1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
g1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g2)

-- | Generate a random number from the Weibull distribution.
generateWeibull01 :: IO Double
                     -- ^ an uniform random number ~ U (0, 1)
                     -> Double
                     -- ^ shape
                     -> Double
                     -- ^ scale
                     -> IO Double
generateWeibull01 :: IO Double -> Double -> Double -> IO Double
generateWeibull01 IO Double
g Double
alpha Double
beta =
  do Double
x <- IO Double
g
     Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
beta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (- Double -> Double
forall a. Floating a => a -> a
log Double
x) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
alpha)

-- | Generate a random value from the specified discrete distribution.
generateDiscrete01 :: IO Double
                      -- ^ an uniform random number ~ U (0, 1)
                      -> DiscretePDF a
                      -- ^ a discrete probability density function
                      -> IO a
generateDiscrete01 :: IO Double -> DiscretePDF a -> IO a
generateDiscrete01 IO Double
g []   = [Char] -> IO a
forall a. HasCallStack => [Char] -> a
error [Char]
"Empty PDF: generateDiscrete01"
generateDiscrete01 IO Double
g DiscretePDF a
dpdf =
  do Double
x <- IO Double
g
     let loop :: Double -> [(p, Double)] -> p
loop Double
acc [(p
a, Double
p)] = p
a
         loop Double
acc ((p
a, Double
p) : [(p, Double)]
dpdf) =
           if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
acc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
p
           then p
a
           else Double -> [(p, Double)] -> p
loop (Double
acc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
p) [(p, Double)]
dpdf
     a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> IO a) -> a -> IO a
forall a b. (a -> b) -> a -> b
$ Double -> DiscretePDF a -> a
forall p. Double -> [(p, Double)] -> p
loop Double
0 DiscretePDF a
dpdf

-- | Defines a type of the random number generator.
data GeneratorType = SimpleGenerator
                     -- ^ The simple random number generator.
                   | SimpleGeneratorWithSeed Word32
                     -- ^ The simple random number generator with the specified seed.
                   | CustomGenerator (IO Generator)
                     -- ^ The custom random number generator.
                   | CustomGenerator01 (IO Double)
                     -- ^ The custom random number generator by the specified uniform
                     -- generator of numbers from 0 to 1.

-- | Create a new random number generator by the specified type.
newGenerator :: GeneratorType -> IO Generator
newGenerator :: GeneratorType -> IO Generator
newGenerator GeneratorType
tp =
  case GeneratorType
tp of
    GeneratorType
SimpleGenerator ->
      Gen RealWorld -> IO Double
forall a (m :: * -> *).
(Variate a, PrimMonad m) =>
Gen (PrimState m) -> m a
MWC.uniform (Gen RealWorld -> IO Double)
-> IO (Gen RealWorld) -> IO (IO Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> IO (Gen RealWorld)
IO GenIO
MWC.createSystemRandom IO (IO Double) -> (IO Double -> IO Generator) -> IO Generator
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= IO Double -> IO Generator
newRandomGenerator01
    SimpleGeneratorWithSeed Word32
x ->
      Gen RealWorld -> IO Double
forall a (m :: * -> *).
(Variate a, PrimMonad m) =>
Gen (PrimState m) -> m a
MWC.uniform (Gen RealWorld -> IO Double)
-> IO (Gen RealWorld) -> IO (IO Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Vector Word32 -> IO GenIO
forall (m :: * -> *) (v :: * -> *).
(PrimMonad m, Vector v Word32) =>
v Word32 -> m (Gen (PrimState m))
MWC.initialize (Word32 -> Vector Word32
forall a. a -> Vector a
singleton Word32
x) IO (IO Double) -> (IO Double -> IO Generator) -> IO Generator
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= IO Double -> IO Generator
newRandomGenerator01
    CustomGenerator IO Generator
g ->
      IO Generator
g
    CustomGenerator01 IO Double
g ->
      IO Double -> IO Generator
newRandomGenerator01 IO Double
g

-- | Create a new random generator by the specified standard generator.
newRandomGenerator :: RandomGen g => g -> IO Generator
newRandomGenerator :: g -> IO Generator
newRandomGenerator g
g =
  do IORef g
r <- g -> IO (IORef g)
forall a. a -> IO (IORef a)
newIORef g
g
     let g1 :: IO Double
g1 = do g
g <- IORef g -> IO g
forall a. IORef a -> IO a
readIORef IORef g
r
                 let (Double
x, g
g') = g -> (Double, g)
forall a g. (Random a, RandomGen g) => g -> (a, g)
random g
g
                 IORef g -> g -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef g
r g
g'
                 Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
x
     IO Double -> IO Generator
newRandomGenerator01 IO Double
g1

-- | Create a new random generator by the specified uniform generator of numbers from 0 to 1.
newRandomGenerator01 :: IO Double -> IO Generator
newRandomGenerator01 :: IO Double -> IO Generator
newRandomGenerator01 IO Double
g =
  do let g1 :: IO Double
g1 = IO Double
g
     IO Double
g2 <- IO Double -> IO (IO Double)
newNormalGenerator01 IO Double
g1
     let g3 :: Double -> Double -> IO Double
g3 Double
mu Double
nu = do { Double
x <- IO Double
g2; Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double
mu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
nu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x }
         g4 :: Double -> Double -> IO Double
g4 Double
mu Double
nu = do { Double
x <- IO Double
g2; Double -> IO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
exp (Double
mu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
nu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) }
     IORef Int
seqNoRef <- Int -> IO (IORef Int)
forall a. a -> IO (IORef a)
newIORef Int
0
     let seqNo :: IO Int
seqNo = do { Int
x <- IORef Int -> IO Int
forall a. IORef a -> IO a
readIORef IORef Int
seqNoRef; IORef Int -> (Int -> Int) -> IO ()
forall a. IORef a -> (a -> a) -> IO ()
modifyIORef' IORef Int
seqNoRef (Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1); Int -> IO Int
forall (m :: * -> *) a. Monad m => a -> m a
return Int
x }
     Generator -> IO Generator
forall (m :: * -> *) a. Monad m => a -> m a
return Generator :: (Double -> Double -> IO Double)
-> (Int -> Int -> IO Int)
-> (Double -> Double -> Double -> IO Double)
-> (Double -> Double -> IO Double)
-> (Double -> Double -> IO Double)
-> (Double -> IO Double)
-> (Double -> Int -> IO Double)
-> (Double -> IO Int)
-> (Double -> Int -> IO Int)
-> (Double -> Double -> IO Double)
-> (Double -> Double -> IO Double)
-> (Double -> Double -> IO Double)
-> (forall a. DiscretePDF a -> IO a)
-> IO Int
-> Generator
Generator { generateUniform :: Double -> Double -> IO Double
generateUniform = IO Double -> Double -> Double -> IO Double
generateUniform01 IO Double
g1,
                        generateUniformInt :: Int -> Int -> IO Int
generateUniformInt = IO Double -> Int -> Int -> IO Int
generateUniformInt01 IO Double
g1,
                        generateTriangular :: Double -> Double -> Double -> IO Double
generateTriangular = IO Double -> Double -> Double -> Double -> IO Double
generateTriangular01 IO Double
g1,
                        generateNormal :: Double -> Double -> IO Double
generateNormal = Double -> Double -> IO Double
g3,
                        generateLogNormal :: Double -> Double -> IO Double
generateLogNormal = Double -> Double -> IO Double
g4,
                        generateExponential :: Double -> IO Double
generateExponential = IO Double -> Double -> IO Double
generateExponential01 IO Double
g1,
                        generateErlang :: Double -> Int -> IO Double
generateErlang = IO Double -> Double -> Int -> IO Double
generateErlang01 IO Double
g1,
                        generatePoisson :: Double -> IO Int
generatePoisson = IO Double -> Double -> IO Int
generatePoisson01 IO Double
g1,
                        generateBinomial :: Double -> Int -> IO Int
generateBinomial = IO Double -> Double -> Int -> IO Int
generateBinomial01 IO Double
g1,
                        generateGamma :: Double -> Double -> IO Double
generateGamma = IO Double -> IO Double -> Double -> Double -> IO Double
generateGamma01 IO Double
g2 IO Double
g1,
                        generateBeta :: Double -> Double -> IO Double
generateBeta = IO Double -> IO Double -> Double -> Double -> IO Double
generateBeta01 IO Double
g2 IO Double
g1,
                        generateWeibull :: Double -> Double -> IO Double
generateWeibull = IO Double -> Double -> Double -> IO Double
generateWeibull01 IO Double
g1,
                        generateDiscrete :: forall a. DiscretePDF a -> IO a
generateDiscrete = IO Double -> DiscretePDF a -> IO a
forall a. IO Double -> DiscretePDF a -> IO a
generateDiscrete01 IO Double
g1,
                        generateSequenceNo :: IO Int
generateSequenceNo = IO Int
seqNo }