{-# LANGUAGE TypeFamilies #-}

-- |
-- Module     : Simulation.Aivika.Lattice.Generator
-- Copyright  : Copyright (c) 2016-2017, David Sorokin <david.sorokin@gmail.com>
-- License    : BSD3
-- Maintainer : David Sorokin <david.sorokin@gmail.com>
-- Stability  : experimental
-- Tested with: GHC 7.10.3
--
-- Here is defined a random number generator,
-- where 'LIO' is an instance of 'MonadGenerator'.
--
module Simulation.Aivika.Lattice.Generator () where

import Control.Monad
import Control.Monad.Trans

import System.Random
import qualified System.Random.MWC as MWC

import Data.IORef

import Simulation.Aivika.Trans
import Simulation.Aivika.Trans.Generator.Primitive
import Simulation.Aivika.Lattice.Internal.LIO

instance MonadGenerator LIO where

  data Generator LIO =
    GeneratorLIO { Generator LIO -> LIO Double
generator01 :: LIO Double,
                   -- ^ the generator of uniform numbers from 0 to 1
                   Generator LIO -> LIO Double
generatorNormal01 :: LIO Double,
                   -- ^ the generator of normal numbers with mean 0 and variance 1
                   Generator LIO -> LIO Int
generatorSequenceNo :: LIO Int
                   -- ^ the generator of sequence numbers
                 }

  generateUniform :: Generator LIO -> Double -> Double -> LIO Double
generateUniform = LIO Double -> Double -> Double -> LIO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateUniform01 (LIO Double -> Double -> Double -> LIO Double)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> Double
-> LIO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generateUniformInt :: Generator LIO -> Int -> Int -> LIO Int
generateUniformInt = LIO Double -> Int -> Int -> LIO Int
forall (m :: * -> *). Monad m => m Double -> Int -> Int -> m Int
generateUniformInt01 (LIO Double -> Int -> Int -> LIO Int)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Int
-> Int
-> LIO Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generateTriangular :: Generator LIO -> Double -> Double -> Double -> LIO Double
generateTriangular = LIO Double -> Double -> Double -> Double -> LIO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> Double -> m Double
generateTriangular01 (LIO Double -> Double -> Double -> Double -> LIO Double)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> Double
-> Double
-> LIO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generateNormal :: Generator LIO -> Double -> Double -> LIO Double
generateNormal = LIO Double -> Double -> Double -> LIO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateNormal01 (LIO Double -> Double -> Double -> LIO Double)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> Double
-> LIO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generatorNormal01

  generateLogNormal :: Generator LIO -> Double -> Double -> LIO Double
generateLogNormal = LIO Double -> Double -> Double -> LIO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateLogNormal01 (LIO Double -> Double -> Double -> LIO Double)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> Double
-> LIO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generatorNormal01

  generateExponential :: Generator LIO -> Double -> LIO Double
generateExponential = LIO Double -> Double -> LIO Double
forall (m :: * -> *). Monad m => m Double -> Double -> m Double
generateExponential01 (LIO Double -> Double -> LIO Double)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> LIO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generateErlang :: Generator LIO -> Double -> Int -> LIO Double
generateErlang = LIO Double -> Double -> Int -> LIO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Int -> m Double
generateErlang01 (LIO Double -> Double -> Int -> LIO Double)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> Int
-> LIO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generatePoisson :: Generator LIO -> Double -> LIO Int
generatePoisson = LIO Double -> Double -> LIO Int
forall (m :: * -> *). Monad m => m Double -> Double -> m Int
generatePoisson01 (LIO Double -> Double -> LIO Int)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> LIO Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generateBinomial :: Generator LIO -> Double -> Int -> LIO Int
generateBinomial = LIO Double -> Double -> Int -> LIO Int
forall (m :: * -> *). Monad m => m Double -> Double -> Int -> m Int
generateBinomial01 (LIO Double -> Double -> Int -> LIO Int)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> Int
-> LIO Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generateGamma :: Generator LIO -> Double -> Double -> LIO Double
generateGamma Generator LIO
g = LIO Double -> LIO Double -> Double -> Double -> LIO Double
forall (m :: * -> *).
Monad m =>
m Double -> m Double -> Double -> Double -> m Double
generateGamma01 (Generator LIO -> LIO Double
generatorNormal01 Generator LIO
g) (Generator LIO -> LIO Double
generator01 Generator LIO
g)

  generateBeta :: Generator LIO -> Double -> Double -> LIO Double
generateBeta Generator LIO
g = LIO Double -> LIO Double -> Double -> Double -> LIO Double
forall (m :: * -> *).
Monad m =>
m Double -> m Double -> Double -> Double -> m Double
generateBeta01 (Generator LIO -> LIO Double
generatorNormal01 Generator LIO
g) (Generator LIO -> LIO Double
generator01 Generator LIO
g)

  generateWeibull :: Generator LIO -> Double -> Double -> LIO Double
generateWeibull = LIO Double -> Double -> Double -> LIO Double
forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateWeibull01 (LIO Double -> Double -> Double -> LIO Double)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> Double
-> Double
-> LIO Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generateDiscrete :: Generator LIO -> DiscretePDF a -> LIO a
generateDiscrete = LIO Double -> DiscretePDF a -> LIO a
forall (m :: * -> *) a. Monad m => m Double -> DiscretePDF a -> m a
generateDiscrete01 (LIO Double -> DiscretePDF a -> LIO a)
-> (Generator LIO -> LIO Double)
-> Generator LIO
-> DiscretePDF a
-> LIO a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01

  generateSequenceNo :: Generator LIO -> LIO Int
generateSequenceNo = Generator LIO -> LIO Int
generatorSequenceNo

  newGenerator :: GeneratorType LIO -> LIO (Generator LIO)
newGenerator GeneratorType LIO
tp =
    case GeneratorType LIO
tp of
      GeneratorType LIO
SimpleGenerator ->
        do let g :: IO (IO Double)
g = Gen (PrimState IO) -> IO Double
forall a (m :: * -> *).
(Variate a, PrimMonad m) =>
Gen (PrimState m) -> m a
MWC.uniform (Gen (PrimState IO) -> IO Double)
-> IO (Gen (PrimState IO)) -> IO (IO Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                   IO (Gen (PrimState IO))
MWC.createSystemRandom
           IO Double
g' <- IO (IO Double) -> LIO (IO Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO (IO Double)
g
           LIO Double -> LIO (Generator LIO)
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 (IO Double -> LIO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO Double
g')
      SimpleGeneratorWithSeed Word32
x ->
        [Char] -> LIO (Generator LIO)
forall a. HasCallStack => [Char] -> a
error [Char]
"Unsupported generator type SimpleGeneratorWithSeed: newGenerator"
      CustomGenerator LIO (Generator LIO)
g ->
        LIO (Generator LIO)
g
      CustomGenerator01 LIO Double
g ->
        LIO Double -> LIO (Generator LIO)
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 LIO Double
g

  newRandomGenerator :: g -> LIO (Generator LIO)
newRandomGenerator g
g = 
    do IORef g
r <- IO (IORef g) -> LIO (IORef g)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef g) -> LIO (IORef g)) -> IO (IORef g) -> LIO (IORef g)
forall a b. (a -> b) -> a -> b
$ g -> IO (IORef g)
forall a. a -> IO (IORef a)
newIORef g
g
       let g01 :: LIO Double
g01 = do g
g <- IO g -> LIO g
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO g -> LIO g) -> IO g -> LIO g
forall a b. (a -> b) -> a -> b
$ 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
                    IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef g -> g -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef g
r g
g'
                    Double -> LIO Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
x
       LIO Double -> LIO (Generator LIO)
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 LIO Double
g01

  newRandomGenerator01 :: LIO Double -> LIO (Generator LIO)
newRandomGenerator01 LIO Double
g01 =
    do LIO Double
gNormal01 <- LIO Double -> LIO (LIO Double)
newNormalGenerator01 LIO Double
g01
       IORef Int
gSeqNoRef <- IO (IORef Int) -> LIO (IORef Int)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Int) -> LIO (IORef Int))
-> IO (IORef Int) -> LIO (IORef Int)
forall a b. (a -> b) -> a -> b
$ Int -> IO (IORef Int)
forall a. a -> IO (IORef a)
newIORef Int
0
       let gSeqNo :: LIO Int
gSeqNo =
             do Int
x <- IO Int -> LIO Int
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Int -> LIO Int) -> IO Int -> LIO Int
forall a b. (a -> b) -> a -> b
$ IORef Int -> IO Int
forall a. IORef a -> IO a
readIORef IORef Int
gSeqNoRef
                IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Int -> (Int -> Int) -> IO ()
forall a. IORef a -> (a -> a) -> IO ()
modifyIORef' IORef Int
gSeqNoRef (Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                Int -> LIO Int
forall (m :: * -> *) a. Monad m => a -> m a
return Int
x
       Generator LIO -> LIO (Generator LIO)
forall (m :: * -> *) a. Monad m => a -> m a
return GeneratorLIO :: LIO Double -> LIO Double -> LIO Int -> Generator LIO
GeneratorLIO { generator01 :: LIO Double
generator01 = LIO Double
g01,
                             generatorNormal01 :: LIO Double
generatorNormal01 = LIO Double
gNormal01,
                             generatorSequenceNo :: LIO Int
generatorSequenceNo = LIO Int
gSeqNo }

-- | 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 :: LIO Double
                        -- ^ the generator
                        -> LIO (LIO Double)
newNormalGenerator01 :: LIO Double -> LIO (LIO Double)
newNormalGenerator01 LIO Double
g =
  do IORef Double
nextRef <- IO (IORef Double) -> LIO (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> LIO (IORef Double))
-> IO (IORef Double) -> LIO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Bool
flagRef <- IO (IORef Bool) -> LIO (IORef Bool)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Bool) -> LIO (IORef Bool))
-> IO (IORef Bool) -> LIO (IORef Bool)
forall a b. (a -> b) -> a -> b
$ Bool -> IO (IORef Bool)
forall a. a -> IO (IORef a)
newIORef Bool
False
     IORef Double
xi1Ref  <- IO (IORef Double) -> LIO (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> LIO (IORef Double))
-> IO (IORef Double) -> LIO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Double
xi2Ref  <- IO (IORef Double) -> LIO (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> LIO (IORef Double))
-> IO (IORef Double) -> LIO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     IORef Double
psiRef  <- IO (IORef Double) -> LIO (IORef Double)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (IORef Double) -> LIO (IORef Double))
-> IO (IORef Double) -> LIO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef Double
0.0
     let loop :: LIO ()
loop =
           do Double
psi <- IO Double -> LIO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> LIO Double) -> IO Double -> LIO Double
forall a b. (a -> b) -> a -> b
$ 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 <- LIO Double
g
                        Double
g2 <- LIO 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
                        IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
xi1
                        IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
xi2
                        IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
psi
                        LIO ()
loop
                else IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ 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)
     LIO Double -> LIO (LIO Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (LIO Double -> LIO (LIO Double)) -> LIO Double -> LIO (LIO Double)
forall a b. (a -> b) -> a -> b
$
       do Bool
flag <- IO Bool -> LIO Bool
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Bool -> LIO Bool) -> IO Bool -> LIO Bool
forall a b. (a -> b) -> a -> b
$ IORef Bool -> IO Bool
forall a. IORef a -> IO a
readIORef IORef Bool
flagRef
          if Bool
flag
            then do IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Bool -> Bool -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
False
                    IO Double -> LIO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> LIO Double) -> IO Double -> LIO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
nextRef
            else do IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
0.0
                    IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
0.0
                    IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Double -> Double -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
0.0
                    LIO ()
loop
                    Double
xi1 <- IO Double -> LIO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> LIO Double) -> IO Double -> LIO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
xi1Ref
                    Double
xi2 <- IO Double -> LIO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> LIO Double) -> IO Double -> LIO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
xi2Ref
                    Double
psi <- IO Double -> LIO Double
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Double -> LIO Double) -> IO Double -> LIO Double
forall a b. (a -> b) -> a -> b
$ IORef Double -> IO Double
forall a. IORef a -> IO a
readIORef IORef Double
psiRef
                    IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ IORef Bool -> Bool -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
True
                    IO () -> LIO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> LIO ()) -> IO () -> LIO ()
forall a b. (a -> b) -> a -> b
$ 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 -> LIO Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> LIO Double) -> Double -> LIO Double
forall a b. (a -> b) -> a -> b
$ Double
xi1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
psi