module Control.Monad.MC.GSLBase (
MC(..),
STMC,
IOMC,
evalMC,
RNG,
IORNG,
STRNG,
Seed,
mt19937,
mt19937WithState,
getRNGName,
getRNGSize,
getRNGState,
setRNGState,
uniform,
uniformInt,
normal,
exponential,
gamma,
cauchy,
levy,
levySkew,
pareto,
weibull,
logistic,
beta,
bernoulli,
poisson,
dirichlet,
multinomial,
) where
import Control.Applicative ( Applicative(..) )
import Control.Monad ( liftM )
import Control.Monad.Fix ( MonadFix(..) )
import Control.Monad.IO.Class ( MonadIO(..) )
import Control.Monad.ST ( ST, runST )
import Control.Monad.Primitive ( PrimMonad(..), unsafePrimToPrim )
import Control.Monad.Trans.Class ( MonadTrans(..) )
import Data.Typeable ( Typeable )
import Data.Word ( Word8, Word64 )
import qualified Data.Vector.Storable as VS
import qualified GSL.Random.Gen as GSL
import GSL.Random.Dist
newtype MC m a = MC { runMC :: RNG (PrimState m) -> m a }
type STMC s a = MC (ST s) a
type IOMC a = MC IO a
evalMC :: (forall s. STMC s a) -> (forall s. ST s (STRNG s)) -> a
evalMC ma mr = runST $ do
r <- mr
runMC ma r
instance (Functor m) => Functor (MC m) where
fmap f (MC m) = MC $ \r -> fmap f (m r)
instance (Applicative m) => Applicative (MC m) where
pure a = MC $ \_ -> pure a
(MC gf) <*> (MC ga) = MC $ \r -> (gf r) <*> (ga r)
instance (Monad m) => Monad (MC m) where
return a = MC $ \_ -> return a
(MC m) >>= k =
MC $ \r -> m r >>= \a ->
let (MC m') = k a
in m' r
fail msg = MC $ \_ -> fail msg
instance (MonadFix m) => MonadFix (MC m) where
mfix f = MC $ \r -> mfix $ flip (runMC . f) r
instance (MonadIO m) => MonadIO (MC m) where
liftIO io = MC $ \_ -> liftIO io
instance MonadTrans MC where
lift m = MC $ \_ -> m
newtype RNG s = RNG GSL.RNG
deriving(Eq, Show, Typeable)
type IORNG = RNG (PrimState IO)
type STRNG s = RNG (PrimState (ST s))
type Seed = Word64
getRNGName :: (PrimMonad m) => RNG (PrimState m) -> m String
getRNGName (RNG r) = unsafePrimToPrim $ GSL.getName r
getRNGSize :: (PrimMonad m) => RNG (PrimState m) -> m Int
getRNGSize (RNG r) = liftM fromIntegral $ unsafePrimToPrim $ GSL.getSize r
getRNGState :: (PrimMonad m) => RNG (PrimState m) -> m [Word8]
getRNGState (RNG r) = unsafePrimToPrim $ GSL.getState r
setRNGState :: (PrimMonad m) => RNG (PrimState m) -> [Word8] -> m ()
setRNGState (RNG r) xs = unsafePrimToPrim $ GSL.setState r xs
mt19937 :: (PrimMonad m) => Seed -> m (RNG (PrimState m))
mt19937 s = unsafePrimToPrim $ do
r <- GSL.newRNG GSL.mt19937
GSL.setSeed r s
return (RNG r)
mt19937WithState :: (PrimMonad m) => [Word8] -> m (RNG (PrimState m))
mt19937WithState xs = unsafePrimToPrim $ do
r <- GSL.newRNG GSL.mt19937
GSL.setState r xs
return (RNG r)
uniform :: (PrimMonad m) => Double -> Double -> MC m Double
uniform 0 1 = liftRan0 GSL.getUniform
uniform a b = liftRan2 getFlat a b
uniformInt :: (PrimMonad m) => Int -> MC m Int
uniformInt = liftRan1 GSL.getUniformInt
normal :: (PrimMonad m) => Double -> Double -> MC m Double
normal 0 1 = liftRan0 getUGaussianRatioMethod
normal mu 1 = liftM (mu +) $ liftRan0 getUGaussianRatioMethod
normal 0 sigma = liftRan1 getGaussianRatioMethod sigma
normal mu sigma = liftM (mu +) $ liftRan1 getGaussianRatioMethod sigma
exponential :: (PrimMonad m) => Double -> MC m Double
exponential = liftRan1 getExponential
poisson :: (PrimMonad m) => Double -> MC m Int
poisson = liftRan1 getPoisson
levy :: (PrimMonad m) => Double -> Double -> MC m Double
levy = liftRan2 getLevy
levySkew :: (PrimMonad m) => Double -> Double -> Double -> MC m Double
levySkew = liftRan3 getLevySkew
cauchy :: (PrimMonad m) => Double -> MC m Double
cauchy = liftRan1 getCauchy
beta :: (PrimMonad m) => Double -> Double -> MC m Double
beta = liftRan2 getBeta
logistic :: (PrimMonad m) => Double -> MC m Double
logistic = liftRan1 getLogistic
pareto :: (PrimMonad m) => Double -> Double -> MC m Double
pareto = liftRan2 getPareto
weibull :: (PrimMonad m) => Double -> Double -> MC m Double
weibull = liftRan2 getWeibull
gamma :: (PrimMonad m) => Double -> Double -> MC m Double
gamma = liftRan2 getGamma
multinomial :: (PrimMonad m) => Int -> VS.Vector Double -> MC m (VS.Vector Int)
multinomial = liftRan2 getMultinomial
dirichlet :: (PrimMonad m) => VS.Vector Double -> MC m (VS.Vector Double)
dirichlet = liftRan1 getDirichlet
bernoulli :: (PrimMonad m ) => Double -> MC m Bool
bernoulli p = liftM (< p) $ uniform 0 1
liftRan0 :: (PrimMonad m) => (GSL.RNG -> IO a) -> MC m a
liftRan0 f = MC $ \(RNG r) -> unsafePrimToPrim $ f r
liftRan1 :: (PrimMonad m) => (GSL.RNG -> a -> IO b) -> a -> MC m b
liftRan1 f a = MC $ \(RNG r) -> unsafePrimToPrim $ f r a
liftRan2 :: (PrimMonad m) => (GSL.RNG -> a -> b -> IO c) -> a -> b -> MC m c
liftRan2 f a b = MC $ \(RNG r) -> unsafePrimToPrim $ f r a b
liftRan3 :: (PrimMonad m) => (GSL.RNG -> a -> b -> c -> IO d) -> a -> b -> c -> MC m d
liftRan3 f a b c = MC $ \(RNG r) -> unsafePrimToPrim $ f r a b c