{-# LANGUAGE TypeFamilies #-}
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,
Generator LIO -> LIO Double
generatorNormal01 :: LIO Double,
Generator LIO -> LIO Int
generatorSequenceNo :: LIO Int
}
generateUniform :: Generator LIO -> Double -> Double -> LIO Double
generateUniform = forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateUniform01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01
generateUniformInt :: Generator LIO -> Int -> Int -> LIO Int
generateUniformInt = forall (m :: * -> *). Monad m => m Double -> Int -> Int -> m Int
generateUniformInt01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01
generateTriangular :: Generator LIO -> Double -> Double -> Double -> LIO Double
generateTriangular = forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> Double -> m Double
generateTriangular01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01
generateNormal :: Generator LIO -> Double -> Double -> LIO Double
generateNormal = forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateNormal01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generatorNormal01
generateLogNormal :: Generator LIO -> Double -> Double -> LIO Double
generateLogNormal = forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateLogNormal01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generatorNormal01
generateExponential :: Generator LIO -> Double -> LIO Double
generateExponential = forall (m :: * -> *). Monad m => m Double -> Double -> m Double
generateExponential01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01
generateErlang :: Generator LIO -> Double -> Int -> LIO Double
generateErlang = forall (m :: * -> *).
Monad m =>
m Double -> Double -> Int -> m Double
generateErlang01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01
generatePoisson :: Generator LIO -> Double -> LIO Int
generatePoisson = forall (m :: * -> *). Monad m => m Double -> Double -> m Int
generatePoisson01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01
generateBinomial :: Generator LIO -> Double -> Int -> LIO Int
generateBinomial = forall (m :: * -> *). Monad m => m Double -> Double -> Int -> m Int
generateBinomial01 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 = 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 = 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 = forall (m :: * -> *).
Monad m =>
m Double -> Double -> Double -> m Double
generateWeibull01 forall b c a. (b -> c) -> (a -> b) -> a -> c
. Generator LIO -> LIO Double
generator01
generateDiscrete :: forall a. Generator LIO -> DiscretePDF a -> LIO a
generateDiscrete = forall (m :: * -> *) a. Monad m => m Double -> DiscretePDF a -> m a
generateDiscrete01 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 = forall a (m :: * -> *).
(Variate a, PrimMonad m) =>
Gen (PrimState m) -> m a
MWC.uniform forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
IO GenIO
MWC.createSystemRandom
IO Double
g' <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO (IO Double)
g
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 (forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO IO Double
g')
SimpleGeneratorWithSeed Word32
x ->
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 ->
forall (m :: * -> *).
MonadGenerator m =>
m Double -> m (Generator m)
newRandomGenerator01 LIO Double
g
newRandomGenerator :: forall g. RandomGen g => g -> LIO (Generator LIO)
newRandomGenerator g
g =
do IORef g
r <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. a -> IO (IORef a)
newIORef g
g
let g01 :: LIO Double
g01 = do g
g <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef IORef g
r
let (Double
x, g
g') = forall a g. (Random a, RandomGen g) => g -> (a, g)
random g
g
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef g
r g
g'
forall (m :: * -> *) a. Monad m => a -> m a
return Double
x
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 <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. a -> IO (IORef a)
newIORef Int
0
let gSeqNo :: LIO Int
gSeqNo =
do Int
x <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef IORef Int
gSeqNoRef
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> (a -> a) -> IO ()
modifyIORef' IORef Int
gSeqNoRef (forall a. Num a => a -> a -> a
+Int
1)
forall (m :: * -> *) a. Monad m => a -> m a
return Int
x
forall (m :: * -> *) a. Monad m => a -> m a
return GeneratorLIO { generator01 :: LIO Double
generator01 = LIO Double
g01,
generatorNormal01 :: LIO Double
generatorNormal01 = LIO Double
gNormal01,
generatorSequenceNo :: LIO Int
generatorSequenceNo = LIO Int
gSeqNo }
newNormalGenerator01 :: LIO Double
-> LIO (LIO Double)
newNormalGenerator01 :: LIO Double -> LIO (LIO Double)
newNormalGenerator01 LIO Double
g =
do IORef Double
nextRef <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. a -> IO (IORef a)
newIORef Double
0.0
IORef Bool
flagRef <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. a -> IO (IORef a)
newIORef Bool
False
IORef Double
xi1Ref <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. a -> IO (IORef a)
newIORef Double
0.0
IORef Double
xi2Ref <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. a -> IO (IORef a)
newIORef Double
0.0
IORef Double
psiRef <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. a -> IO (IORef a)
newIORef Double
0.0
let loop :: LIO ()
loop =
do Double
psi <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef IORef Double
psiRef
if (Double
psi forall a. Ord a => a -> a -> Bool
>= Double
1.0) Bool -> Bool -> Bool
|| (Double
psi 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 forall a. Num a => a -> a -> a
* Double
g1 forall a. Num a => a -> a -> a
- Double
1.0
xi2 :: Double
xi2 = Double
2.0 forall a. Num a => a -> a -> a
* Double
g2 forall a. Num a => a -> a -> a
- Double
1.0
psi :: Double
psi = Double
xi1 forall a. Num a => a -> a -> a
* Double
xi1 forall a. Num a => a -> a -> a
+ Double
xi2 forall a. Num a => a -> a -> a
* Double
xi2
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
xi1
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
xi2
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
psi
LIO ()
loop
else forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
sqrt (- Double
2.0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
psi forall a. Fractional a => a -> a -> a
/ Double
psi)
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$
do Bool
flag <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef IORef Bool
flagRef
if Bool
flag
then do forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
False
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef IORef Double
nextRef
else do forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi1Ref Double
0.0
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Double
xi2Ref Double
0.0
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Double
psiRef Double
0.0
LIO ()
loop
Double
xi1 <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef IORef Double
xi1Ref
Double
xi2 <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef IORef Double
xi2Ref
Double
psi <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef IORef Double
psiRef
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Bool
flagRef Bool
True
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> a -> IO ()
writeIORef IORef Double
nextRef forall a b. (a -> b) -> a -> b
$ Double
xi2 forall a. Num a => a -> a -> a
* Double
psi
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Double
xi1 forall a. Num a => a -> a -> a
* Double
psi