module Simulation.Aivika.IO.Generator () where
import Control.Monad
import Control.Monad.Trans
import System.Random
import Data.IORef
import Simulation.Aivika.Trans.Generator
import Simulation.Aivika.Trans.Generator.Primitive
import Simulation.Aivika.Trans.Template
instance (Functor m, MonadIO m, MonadTemplate m) => MonadGenerator m where
data Generator m =
Generator { generator01 :: m Double,
generatorNormal01 :: m Double
}
generateUniform = generateUniform01 . generator01
generateUniformInt = generateUniformInt01 . generator01
generateTriangular = generateTriangular01 . generator01
generateNormal = generateNormal01 . generatorNormal01
generateLogNormal = generateLogNormal01 . generatorNormal01
generateExponential = generateExponential01 . generator01
generateErlang = generateErlang01 . generator01
generatePoisson = generatePoisson01 . generator01
generateBinomial = generateBinomial01 . generator01
generateGamma g = generateGamma01 (generatorNormal01 g) (generator01 g)
generateBeta g = generateBeta01 (generatorNormal01 g) (generator01 g)
generateWeibull = generateWeibull01 . generator01
generateDiscrete = generateDiscrete01 . generator01
newGenerator tp =
case tp of
SimpleGenerator ->
liftIO newStdGen >>= newRandomGenerator
SimpleGeneratorWithSeed x ->
newRandomGenerator $ mkStdGen x
CustomGenerator g ->
g
CustomGenerator01 g ->
newRandomGenerator01 g
newRandomGenerator g =
do r <- liftIO $ newIORef g
let g01 = do g <- liftIO $ readIORef r
let (x, g') = random g
liftIO $ writeIORef r g'
return x
newRandomGenerator01 g01
newRandomGenerator01 g01 =
do gNormal01 <- newNormalGenerator01 g01
return Generator { generator01 = g01,
generatorNormal01 = gNormal01 }
newNormalGenerator01 :: MonadIO m
=> m Double
-> m (m Double)
newNormalGenerator01 g =
do nextRef <- liftIO $ newIORef 0.0
flagRef <- liftIO $ newIORef False
xi1Ref <- liftIO $ newIORef 0.0
xi2Ref <- liftIO $ newIORef 0.0
psiRef <- liftIO $ newIORef 0.0
let loop =
do psi <- liftIO $ readIORef psiRef
if (psi >= 1.0) || (psi == 0.0)
then do g1 <- g
g2 <- g
let xi1 = 2.0 * g1 1.0
xi2 = 2.0 * g2 1.0
psi = xi1 * xi1 + xi2 * xi2
liftIO $ writeIORef xi1Ref xi1
liftIO $ writeIORef xi2Ref xi2
liftIO $ writeIORef psiRef psi
loop
else liftIO $ writeIORef psiRef $ sqrt ( 2.0 * log psi / psi)
return $
do flag <- liftIO $ readIORef flagRef
if flag
then do liftIO $ writeIORef flagRef False
liftIO $ readIORef nextRef
else do liftIO $ writeIORef xi1Ref 0.0
liftIO $ writeIORef xi2Ref 0.0
liftIO $ writeIORef psiRef 0.0
loop
xi1 <- liftIO $ readIORef xi1Ref
xi2 <- liftIO $ readIORef xi2Ref
psi <- liftIO $ readIORef psiRef
liftIO $ writeIORef flagRef True
liftIO $ writeIORef nextRef $ xi2 * psi
return $ xi1 * psi