packages feed

random-fu (empty) → 0.0.0.2

raw patch · 24 files changed

+1255/−0 lines, 24 filesdep +basedep +bytestringdep +mersenne-random-pure64setup-changed

Dependencies added: base, bytestring, mersenne-random-pure64, monad-loops, mtl, random, stateref

Files

+ Setup.lhs view
@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain+
+ random-fu.cabal view
@@ -0,0 +1,52 @@+name:                   random-fu+version:                0.0.0.2+stability:              experimental++cabal-version:          >= 1.2+build-type:             Simple++author:                 James Cook <james.cook@usma.edu>+maintainer:             James Cook <james.cook@usma.edu>+license:                PublicDomain+homepage:               http://code.haskell.org/~mokus/random-fu++category:               Math+synopsis:               Random number generation+description:            Random number generation based on orthogonal typeclasses+                        for entropy sources and random variable distributions, all+                        served up on a monadic platter.  Aspires to be useful+                        in an idiomatic way in both \"pure\" and \"impure\" styles,+                        as well as reasonably fast.  May not yet meet the latter+                        goal, but I think the former is starting to shape up+                        nicely.++Library+  hs-source-dirs:       src+  exposed-modules:      Data.Random+                        Data.Random.Distribution+                        Data.Random.Distribution.Bernoulli+                        Data.Random.Distribution.Beta+                        Data.Random.Distribution.Binomial+                        Data.Random.Distribution.Exponential+                        Data.Random.Distribution.Discrete+                        Data.Random.Distribution.Gamma+                        Data.Random.Distribution.Normal+                        Data.Random.Distribution.Poisson+                        Data.Random.Distribution.Triangular+                        Data.Random.Distribution.Uniform+                        Data.Random.Internal.Classification+                        Data.Random.Internal.Words+                        Data.Random.RVar+                        Data.Random.Source+                        Data.Random.Source.DevRandom+                        Data.Random.Source.StdGen+                        Data.Random.Source.PureMT+                        Data.Random.Source.Std+                        +  build-depends:        base >= 3,+                        bytestring,+                        mersenne-random-pure64,+                        monad-loops >= 0.3.0.1,+                        mtl,+                        random,+                        stateref
+ src/Data/Random.hs view
@@ -0,0 +1,58 @@+{-+ -      ``Data/Random''+ -}+{-# LANGUAGE+    FlexibleContexts+  #-}++-- |Random numbers and stuff...+-- +-- Data.Random.Source exports the typeclasses for entropy sources, and+-- Data.Random.Source.* export various instances and/or functions with which+-- instances can be defined.+-- +-- Data.Random.Distribution exports the typeclasses for sampling distributions,+-- and Data.Random.Distribution.* export various specific distributions.+--+-- Data.Random.RVar exports the 'RVar' type, which is a probability distribution+-- monad that allows for concise definitions of random variables, as well as+-- a couple handy 'RVar's.++module Data.Random+    ( module Data.Random.Source+    , module Data.Random.Source.DevRandom+    , module Data.Random.Source.StdGen+    , module Data.Random.Source.PureMT+    , module Data.Random.Source.Std+    , module Data.Random.Distribution+    , module Data.Random.Distribution.Bernoulli+    , module Data.Random.Distribution.Beta+    , module Data.Random.Distribution.Binomial+    , module Data.Random.Distribution.Discrete+    , module Data.Random.Distribution.Gamma+    , module Data.Random.Distribution.Exponential+    , module Data.Random.Distribution.Normal+    , module Data.Random.Distribution.Poisson+    , module Data.Random.Distribution.Triangular+    , module Data.Random.Distribution.Uniform+    , module Data.Random.RVar+    ) where++import Data.Random.Source+import Data.Random.Source.DevRandom+import Data.Random.Source.StdGen+import Data.Random.Source.PureMT+import Data.Random.Source.Std+import Data.Random.Distribution+import Data.Random.Distribution.Bernoulli+import Data.Random.Distribution.Beta+import Data.Random.Distribution.Binomial+import Data.Random.Distribution.Discrete+import Data.Random.Distribution.Gamma+import Data.Random.Distribution.Exponential+import Data.Random.Distribution.Normal+import Data.Random.Distribution.Poisson+import Data.Random.Distribution.Triangular+import Data.Random.Distribution.Uniform+import Data.Random.RVar+
+ src/Data/Random/Distribution.hs view
@@ -0,0 +1,32 @@+{-+ -      ``Data/Random/Distribution''+ -}+{-# LANGUAGE+    MultiParamTypeClasses, FlexibleContexts+  #-}++module Data.Random.Distribution where++import {-# SOURCE #-} Data.Random.RVar+import Data.Random.Source+import Data.Random.Source.Std+import Data.Word++-- |A definition of a random variable's distribution.  From the distribution+-- an 'RVar' can be created, or the distribution can be directly sampled.+-- 'RVar' in particular is an instance of 'Distribution', and so can be 'sample'd.+--+-- Minimum instance definition: either 'rvar' or 'sampleFrom'.+class Distribution d t where+    -- |Return a random variable with this distribution.+    rvar :: d t -> RVar t+    rvar = sampleFrom StdRandom+    +    -- |Directly sample from the distribution, given a source of entropy.+    sampleFrom :: RandomSource m s => s -> d t -> m t+    sampleFrom src dist = sampleFrom src (rvar dist)++-- |Sample a distribution using the default source of entropy for the+-- monad in which the sampling occurs.+sample :: (Distribution d t, MonadRandom m) => d t -> m t+sample = sampleFrom StdRandom
+ src/Data/Random/Distribution.hs-boot view
@@ -0,0 +1,10 @@+{-+ -      ``Data/Random/Distribution''+ -}+{-# LANGUAGE+    MultiParamTypeClasses, KindSignatures+  #-}++module Data.Random.Distribution where++class Distribution (d :: * -> *) t
+ src/Data/Random/Distribution/Bernoulli.hs view
@@ -0,0 +1,47 @@+{-+ -      ``Data/Random/Distribution/Bernoulli''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts,+    UndecidableInstances+  #-}++module Data.Random.Distribution.Bernoulli where++import Data.Random.Internal.Classification++import Data.Random.Source+import Data.Random.Distribution+import Data.Random.RVar++import Data.Random.Distribution.Uniform++import Data.Int+import Data.Word++bernoulli :: (Distribution (Bernoulli b) a) => b -> RVar a+bernoulli p = sample (Bernoulli p)++boolBernoulli p = do+    x <- realFloatUniform 0 1+    return (x <= p)++generalBernoulli t f p = do+    x <- boolBernoulli p+    return (if x then t else f)++class (Classification NumericType t c) => BernoulliByClassification c t where+    bernoulliByClassification :: RealFloat a => a -> RVar t++instance (Classification NumericType t IntegralType, Num t) => BernoulliByClassification IntegralType t+    where bernoulliByClassification = generalBernoulli 0 1+instance (Classification NumericType t FractionalType, Num t) => BernoulliByClassification FractionalType t+    where bernoulliByClassification = generalBernoulli 0 1+instance (Classification NumericType t EnumType, Enum t) => BernoulliByClassification EnumType t+    where bernoulliByClassification = generalBernoulli (toEnum 0) (toEnum 1)++data Bernoulli b a = Bernoulli b++instance (BernoulliByClassification c t, RealFloat b) => Distribution (Bernoulli b) t where+    rvar (Bernoulli p) = bernoulliByClassification p
+ src/Data/Random/Distribution/Beta.hs view
@@ -0,0 +1,39 @@+{-+ -      ``Data/Random/Distribution/Beta''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts,+    UndecidableInstances+  #-}++module Data.Random.Distribution.Beta where++import Data.Random.Source+import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Gamma+import Data.Random.Distribution.Uniform++import Control.Monad++realFloatBeta :: RealFloat a => a -> a -> RVar a+realFloatBeta 1 1 = realFloatStdUniform+realFloatBeta a b = do+    x <- realFloatGamma a 1+    y <- realFloatGamma b 1+    return (x / (x + y))++realFloatBetaFromIntegral :: (Integral a, Integral b, RealFloat c) => a -> b -> RVar c+realFloatBetaFromIntegral a b =  do+    x <- realFloatErlang a+    y <- realFloatErlang b+    return (x / (x + y))++beta :: Distribution Beta a => a -> a -> RVar a+beta a b = sample (Beta a b)++data Beta a = Beta a a++instance (RealFloat a) => Distribution Beta a where+    rvar (Beta a b) = realFloatBeta a b
+ src/Data/Random/Distribution/Binomial.hs view
@@ -0,0 +1,64 @@+{-+ -      ``Data/Random/Distribution/Binomial''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts,+    UndecidableInstances+  #-}++module Data.Random.Distribution.Binomial where++import Data.Random.Internal.Classification++import Data.Random.Source+import Data.Random.Distribution+import Data.Random.RVar++import Data.Random.Distribution.Beta+import Data.Random.Distribution.Uniform++import Data.Int+import Data.Word+import Control.Monad++    -- algorithm from Knuth's TAOCP, 3rd ed., p 136+    -- specific choice of cutoff size taken from gsl source+integralBinomial :: (Integral a, RealFloat b) => a -> b -> RVar a+integralBinomial t p = bin 0 t p+    where+        bin k t p+            | t > 10    = do+                let a = 1 + t `div` 2+                    b = 1 + t - a+        +                x <- realFloatBetaFromIntegral a b+                if x >= p+                    then bin  k      (a - 1) (p / x)+                    else bin (k + a) (b - 1) ((p - x) / (1 - x))+        +            | otherwise = count k t+                where+                    count k  0    = return k+                    count k (n+1) = do+                        x <- realFloatStdUniform+                        (count $! (if x < p then k + 1 else k)) n++++binomial :: Distribution (Binomial b) a => a -> b -> RVar a+binomial t p = sample (Binomial t p)++class (Classification NumericType t c) => BinomialByClassification c t where+    binomialByClassification :: RealFloat a => t -> a -> RVar t++instance (Classification NumericType t IntegralType, Integral t) => BinomialByClassification IntegralType t+    where binomialByClassification = integralBinomial+instance (Classification NumericType t FractionalType, RealFrac t) => BinomialByClassification FractionalType t+    where binomialByClassification t p = liftM fromInteger (integralBinomial (truncate t) p)++instance (BinomialByClassification c t, RealFloat b) => Distribution (Binomial b) t where+    rvar (Binomial t p) = binomialByClassification t p++data Binomial b a = Binomial a b+
+ src/Data/Random/Distribution/Discrete.hs view
@@ -0,0 +1,35 @@+{-+ -      ``Data/Random/Distribution/Discrete''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts+  #-}++module Data.Random.Distribution.Discrete where++import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform++import Control.Monad++discrete :: Distribution (Discrete p) a => [(p,a)] -> RVar a+discrete ps = rvar (Discrete ps)++data Discrete p a = Discrete [(p, a)]++instance (Num p, Ord p, Distribution Uniform p) => Distribution (Discrete p) a where+    rvar (Discrete []) = fail "discrete distribution over empty set cannot be sampled"+    rvar (Discrete ds) = do+        let (ps, xs) = unzip ds+            cs = scanl1 (+) ps+        +        when (any (<0) ps) $ fail "negative probability in discrete distribution"+        +        u <- uniform 0 (last cs)+        return $ head+            [ x+            | (c,x) <- zip cs xs+            , c >= u+            ]
+ src/Data/Random/Distribution/Exponential.hs view
@@ -0,0 +1,28 @@+{-+ -      ``Data/Random/Distribution/Exponential''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts,+    UndecidableInstances+  #-}++module Data.Random.Distribution.Exponential where++import Data.Random.Source+import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform++data Exponential a = Exp a++realFloatExponential :: RealFloat a => a -> RVar a+realFloatExponential lambdaRecip = do+    x <- realFloatStdUniform+    return (negate (log x) * lambdaRecip)++exponential :: Distribution Exponential a => a -> RVar a+exponential = sample . Exp++instance (RealFloat a) => Distribution Exponential a where+    rvar (Exp lambdaRecip) = realFloatExponential lambdaRecip
+ src/Data/Random/Distribution/Gamma.hs view
@@ -0,0 +1,74 @@+{-+ -      ``Data/Random/Distribution/Gamma''+ -+ -  needs cleanup, verification, and automagic selection of appropriate+ -  algorithms, and proper citations.+ -+ -  should eliminate spurious 'border crossings' betweer RVars and sampleFrom+ -  perhaps Distribution class should have as its basis a function of type+ -  (d t -> RVar t)+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts,+    UndecidableInstances+  #-}++module Data.Random.Distribution.Gamma where++import Data.Random.Source+import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform+import Data.Random.Distribution.Normal++import Control.Monad++    -- translated from gsl source - seems to be best I've found by far.+    -- originally comes from Marsaglia & Tang, "A Simple Method for+    -- generating gamma variables", ACM Transactions on Mathematical+    -- Software, Vol 26, No 3 (2000), p363-372.+realFloatGamma :: RealFloat a => a -> a -> RVar a+realFloatGamma a b+    | a < 1 +    = do+        u <- realFloatStdUniform+        x <- realFloatGamma (1 + a) b+        return (x * u ** recip a)+    | otherwise+    = go+        where+            d = a - (1 / 3)+            c = recip (3 * sqrt d) -- (1 / 3) / sqrt d+            +            go = do+                x <- realFloatStdNormal+                let cx = c * x+                    v = (1 + cx) ^ 3+                    +                    x_2 = x * x+                    x_4 = x_2 * x_2+                +                if cx <= (-1)+                    then go+                    else do+                        u <- realFloatStdUniform+                        +                        if         u < 1 - 0.0331 * x_4+                            || log u < 0.5 * x_2  + d * (1 - v + log v)+                            then return (b * d * v)+                            else go++realFloatErlang :: (Integral a, RealFloat b) => a -> RVar b+realFloatErlang a = realFloatGamma (fromIntegral a) 1++gamma :: (Distribution Gamma a) => a -> a -> RVar a+gamma a b = sample (Gamma a b)++erlang :: (Distribution Gamma a, Integral b, Num a) => b -> a -> RVar a+erlang a b = sample (Gamma (fromIntegral a) b)++data Gamma a = Gamma a a++instance RealFloat a => Distribution Gamma a where+    rvar (Gamma a b) = realFloatGamma a b
+ src/Data/Random/Distribution/Normal.hs view
@@ -0,0 +1,70 @@+{-+ -      ``Data/Random/Distribution/Normal''+ -  + -  Quick and dirty implementation - eventually something faster probably + -  ought to be implemented instead.+ -}+{-# LANGUAGE+    MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,+    UndecidableInstances+  #-}++module Data.Random.Distribution.Normal where++import Data.Random.Source+import Data.Random.Distribution+import Data.Random.Distribution.Uniform+import Data.Random.RVar++import Control.Monad++-- Box-Muller method+normalPair :: (Floating a, Distribution Uniform a) => RVar (a,a)+normalPair = do+    u <- uniform 0 1+    t <- uniform 0 (2 * pi)+    let r = sqrt (-2 * log u)+        +        x = r * cos t+        y = r * sin t+    return (x,y)++-- slightly slower+knuthPolarNormalPair :: (Floating a, Ord a, Distribution Uniform a) => RVar (a,a)+knuthPolarNormalPair = do+    v1 <- uniform (-1) 1+    v2 <- uniform (-1) 1+    +    let s = v1*v1 + v2*v2+    if s >= 1+        then knuthPolarNormalPair+        else return $ if s == 0+            then (0,0)+            else let scale = sqrt (-2 * log s / s) +                  in (v1 * scale, v2 * scale)++realFloatStdNormal :: RealFloat a => RVar a+realFloatStdNormal = do+    u <- realFloatStdUniform+    t <- realFloatStdUniform+    let r = sqrt (-2 * log u)+        +        x = r * cos (t * 2 * pi)+    return x+    ++data Normal a+    = StdNormal+    | Normal a a -- mean, sd++instance (Floating a, Distribution Uniform a) => Distribution Normal a where+    rvar StdNormal = liftM fst normalPair+    rvar (Normal m s) = do+        x <- liftM fst normalPair+        return (x * s + m)++stdNormal :: Distribution Normal a => RVar a+stdNormal = rvar StdNormal++normal :: Distribution Normal a => a -> a -> RVar a+normal m s = rvar (Normal m s)
+ src/Data/Random/Distribution/Poisson.hs view
@@ -0,0 +1,67 @@+{-+ -      ``Data/Random/Distribution/Poisson''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts, UndecidableInstances+  #-}++module Data.Random.Distribution.Poisson where++import Data.Random.Source+import Data.Random.Distribution+import Data.Random.RVar++import Data.Random.Distribution.Uniform+import Data.Random.Distribution.Gamma+import Data.Random.Distribution.Binomial++import Data.Int+import Data.Word++import Control.Monad++-- from Knuth, with interpretation help from gsl sources+integralPoisson :: (Integral a, RealFloat b) => b -> RVar a+integralPoisson mu = psn 0 mu+    where+        psn k mu+            | mu > 10   = do+                let m = floor (mu * (7/8))+            +                x <- realFloatErlang m+                if x >= mu+                    then do+                        b <- integralBinomial (m - 1) (mu / x)+                        return (k + b)+                    else psn (k + m) (mu - x)+            +            | otherwise = prod 1 k+                where+                    emu = exp (-mu)+                +                    prod p k = do+                        u <- realFloatStdUniform+                        if p * u > emu+                            then prod (p * u) (k + 1)+                            else return k+++poisson :: (Distribution (Poisson b) a) => b -> RVar a+poisson mu = sample (Poisson mu)++data Poisson b a = Poisson b++instance RealFloat b => Distribution (Poisson b) Int        where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Int8       where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Int16      where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Int32      where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Int64      where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Word8      where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Word16     where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Word32     where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Word64     where rvar (Poisson mu) = integralPoisson mu+instance RealFloat b => Distribution (Poisson b) Integer    where rvar (Poisson mu) = integralPoisson mu++instance RealFloat b => Distribution (Poisson b) Float      where rvar (Poisson mu) = liftM fromIntegral (integralPoisson mu)+instance RealFloat b => Distribution (Poisson b) Double     where rvar (Poisson mu) = liftM fromIntegral (integralPoisson mu)
+ src/Data/Random/Distribution/Triangular.hs view
@@ -0,0 +1,36 @@+{-+ -      ``Data/Random/Distribution/Triangular''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances+  #-}++module Data.Random.Distribution.Triangular where++import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform++data Triangular a = Triangular+    { triLower  :: a+    , triMid    :: a+    , triUpper  :: a+    } deriving (Eq, Show)++realFloatTriangular :: (RealFloat a) => a -> a -> a -> RVar a+realFloatTriangular a b c+    | a <= b && b <= c+    = do+        let p = (c-b)/(c-a)+        u <- realFloatStdUniform+        let d   | u >= p    = a+                | otherwise = c+            x   | u >= p    = (u - p) / (1 - p)+                | otherwise = u / p+-- may prefer this: reusing u costs resolution, especially if p or 1-p is small and c-a is large.+--        x <- realFloatStdUniform+        return (b - ((1 - sqrt x) * (b-d)))++instance RealFloat a => Distribution Triangular a where+    rvar (Triangular a b c) = realFloatTriangular a b c
+ src/Data/Random/Distribution/Uniform.hs view
@@ -0,0 +1,122 @@+{-+ -      ``Data/Random/Distribution/Uniform''+ -}+{-# LANGUAGE+    MultiParamTypeClasses, FunctionalDependencies,+    FlexibleContexts, FlexibleInstances, +    UndecidableInstances+  #-}++module Data.Random.Distribution.Uniform+    ( Uniform(..)+	, UniformByClassification(..)+	, uniform+	+    , StdUniform(..)+    , StdUniformByClassification(..)+    , stdUniform+    +    , integralUniform+    , realFloatUniform+    +    , boundedStdUniform+    , boundedEnumStdUniform+    , realFloatStdUniform+    ) where++import Data.Random.Internal.Classification++import Data.Random.Source+import Data.Random.Distribution+import Data.Random.RVar++import Data.Word+import Data.Int+import Data.Bits+import Data.List++import Control.Monad.Loops++integralUniform a b+    | a > b     = compute b a+    | otherwise = compute a b+    where+        compute a b = do+            let m = 1 + toInteger b - toInteger a+            +            let bytes = bytesNeeded m+                maxXpossible = (powersOf256 !! bytes) - 1+            +            x <- iterateUntil (maxXpossible - maxXpossible `mod` m >) (nByteInteger bytes)+            return (a + fromInteger (x `mod` m))+++bytesNeeded x = case findIndex (> x) powersOf256 of+    Just x -> x+powersOf256 = iterate (256 *) 1++boundedStdUniform :: (Distribution Uniform a, Bounded a) => RVar a+boundedStdUniform = uniform minBound maxBound++boundedEnumStdUniform :: (Enum a, Bounded a) => RVar a+boundedEnumStdUniform = enumUniform minBound maxBound++-- (0,1]+realFloatStdUniform :: RealFloat a => RVar a+realFloatStdUniform | False     = return one+                    | otherwise = do+    let bitsNeeded  = floatDigits one+        (_, e) = decodeFloat one+    +    x <- nBitInteger bitsNeeded+    if x == 0+        then return 1+        else return (encodeFloat x (e-1))+    +    where one = 1++realFloatUniform :: RealFloat a => a -> a -> RVar a+realFloatUniform 0 1 = realFloatStdUniform+realFloatUniform a b = do+    x <- realFloatStdUniform+    return (a + x * (b - a))++enumUniform :: Enum a => a -> a -> RVar a+enumUniform a b = do+    x <- integralUniform (fromEnum a) (fromEnum b)+    return (toEnum x)++uniform :: Distribution Uniform a => a -> a -> RVar a+uniform a b = rvar (Uniform a b)++stdUniform :: Distribution StdUniform a => RVar a+stdUniform = rvar StdUniform++class (Classification NumericType t c) => UniformByClassification c t where+    uniformByClassification :: t -> t -> RVar t++class (Classification NumericType t c) => StdUniformByClassification c t where+    stdUniformByClassification :: RVar t++data Uniform t = Uniform !t !t+data StdUniform t = StdUniform++instance UniformByClassification c t => Distribution Uniform t+    where rvar (Uniform a b) = uniformByClassification a b++instance StdUniformByClassification c t => Distribution StdUniform t+    where rvar _ = stdUniformByClassification++instance (Classification NumericType t IntegralType, Integral t) => UniformByClassification IntegralType t+    where uniformByClassification = integralUniform+instance (Classification NumericType t FractionalType, RealFloat t) => UniformByClassification FractionalType t+    where uniformByClassification = realFloatUniform+instance (Classification NumericType t EnumType, Enum t) => UniformByClassification EnumType t+    where uniformByClassification = enumUniform++instance (Classification NumericType t IntegralType, Integral t, Bounded t) => StdUniformByClassification IntegralType t+    where stdUniformByClassification = boundedStdUniform+instance (Classification NumericType t FractionalType, RealFloat t) => StdUniformByClassification FractionalType t+    where stdUniformByClassification = realFloatStdUniform+instance (Classification NumericType t EnumType, Enum t, Bounded t) => StdUniformByClassification EnumType t+    where stdUniformByClassification = boundedStdUniform
+ src/Data/Random/Internal/Classification.hs view
@@ -0,0 +1,102 @@+{-# LANGUAGE+    MultiParamTypeClasses, FunctionalDependencies,+    EmptyDataDecls+  #-}+{-+ -      ``Data/Random/Internal/Classification''+ -}      +-- | \"Classification systems\" - for a motivating example, see+--  the implementation of the Uniform distribution.  Basically,+--  I would like to make instances like:+--  +--  > instance RealFloat a => Distribution Uniform a where ...+--  > instance Integral a =>  Distribution Uniform a where ...+--  +--  and so on.  However, this is not sound - what happens if someone+--  comes along and makes a type that's an instance of both Integral and+--  RealFloat?+--  +--  So, we introduce a classification system based on phantom types, so+--  that each type can be unambiguously declared to be \"intensionally\"+--  Integral, Floating, or whatever.+--  +--  Now, obviously it'd be nice not to clutter the Distribution typeclass+--  with extra phantom types that the end user shouldn't care about.  Hence+--  the pattern of introducing typeclasses such as "UniformByClassification"+--  +--  Now, if a new type comes along that is Integral, a single declaration+--  of the following form suffices to attach it to all such Distribution+--  instances:+--  +--  > instance Classification NumericType t IntegralType+--  +--  Not quite as automagic as the @Integral a => Distribution foo@  case,+--  but a bit closer.  Not only that, it leaves open the possibility that+--  a user may bring in a type that is \"mostly\" integral, and has an Integral+--  instance, but should be handled differently for purposes of uniform+--  random number generation.  In such a case, the user may introduce a new+--  classification of their own and provide the required instances for that+--  classification.+--  +--  All in all, although it is not yet well-tested, it has the \"feel\" of +--  a good compromise.+module Data.Random.Internal.Classification where++import Data.Int+import Data.Word+import Data.Ratio++-- |classificiation system, experimental+-- +--      * c (a phantom type) is the classification system+--      +--      * t is the type to be classified+--      +--      * tc (a phantom type) is the classification of t according to c+--+-- The functional dependency, aside from being important because the relation+-- is functional, allows the classification system to be \"discharged\" in+-- cases such as the following:+--+-- > class Classification SomeCS t c => FooByClassification t c where ...+-- > instance FooByClassification t c => Foo t where ...+-- +-- Thus the class of interest to the end user need not display anything+-- at all about the classification system, except in the superclasses of +-- the classes in the contexts of some of its instances.+class Classification c t tc | c t -> tc++-- |A simple classification system covering the cases we care+-- about when sampling distributions.  Loosely, these are the reasons we care:+-- +--   * distributions over Fractional types are handled as if the type were continuous.+-- +--   * distributions over Integral types are handled discretely.+-- +--   * distributions over Enum types (which are not Num instances) are handled +--     like Integral types, but require use of 'fromEnum' and/or 'toEnum' to work with them.+data NumericType++data IntegralType+data FractionalType+data EnumType++instance Classification NumericType Int            IntegralType+instance Classification NumericType Int8           IntegralType+instance Classification NumericType Int16          IntegralType+instance Classification NumericType Int32          IntegralType+instance Classification NumericType Int64          IntegralType+instance Classification NumericType Word8          IntegralType+instance Classification NumericType Word16         IntegralType+instance Classification NumericType Word32         IntegralType+instance Classification NumericType Word64         IntegralType+instance Classification NumericType Integer        IntegralType++instance Classification NumericType Float          FractionalType+instance Classification NumericType Double         FractionalType+instance Classification NumericType (Ratio a)      FractionalType++instance Classification NumericType Char           EnumType+instance Classification NumericType Bool           EnumType+instance Classification NumericType ()             EnumType+instance Classification NumericType Ordering       EnumType
+ src/Data/Random/Internal/Words.hs view
@@ -0,0 +1,36 @@+{-+ -      ``Data/Random/Internal/Words''+ -}++-- |A few little functions I found myself writing inline over and over again.+--+-- Note that these need to be checked to ensure proper behavior on big-endian +-- systems.  They are probably not right at the moment.+module Data.Random.Internal.Words where++import Foreign+import GHC.IOBase++import Data.Word+import Control.Monad++wordsToBytes :: [Word64] -> [Word8]+wordsToBytes = concatMap wordToBytes++wordToBytes :: Word64 -> [Word8]+wordToBytes x = unsafePerformIO . allocaBytes 8 $ \p -> do+    poke (castPtr p) x+    mapM (peekElemOff p) [0..7]++bytesToWords :: [Word8] -> [Word64]+bytesToWords = map bytesToWord . chunk 8+    where+        chunk n [] = []+        chunk n xs = case splitAt n xs of+            (ys, zs) -> ys : chunk n zs+++bytesToWord :: [Word8] -> Word64+bytesToWord bs = unsafePerformIO . allocaBytes 8 $ \p -> do+    zipWithM (pokeElemOff p) [0..7] (bs ++ repeat 0)+    peek (castPtr p)
+ src/Data/Random/RVar.hs view
@@ -0,0 +1,100 @@+{-+ -      ``Data/Random/RVar''+ -}+{-# LANGUAGE+    RankNTypes,+    MultiParamTypeClasses,+    FlexibleInstances+  #-}++-- |Random variables.  An 'RVar' is a sampleable random variable.  Because+-- probability distributions form a monad, they are quite easy to work with+-- in the standard Haskell monadic styles.  For examples, see the source for+-- any of the 'Distribution' instances - they all are defined in terms of+-- 'RVar's.+module Data.Random.RVar+    ( RVar+    +    , nByteInteger+    , nBitInteger+    ) where++import Data.Random.Distribution+import Data.Random.Source++import Data.Word+import Data.Bits++import Control.Applicative+import Control.Monad++-- |An opaque type containing a \"random variable\" - a value +-- which depends on the outcome of some random process.+newtype RVar a = RVar { runDistM :: forall m s. RandomSource m s => s -> m a }++instance Functor RVar where+    fmap = liftM++instance Monad RVar where+    return x = RVar (\_ -> return x)+    fail s   = RVar (\_ -> fail s)+    (RVar x) >>= f = RVar (\s -> do+            x <- x s+            case f x of+                RVar y -> y s+        )++instance Applicative RVar where+    pure  = return+    (<*>) = ap++instance Distribution RVar a where+    rvar = id+    sampleFrom src x = runDistM x src++instance MonadRandom RVar where+    getRandomBytes n = RVar (\s -> getRandomBytesFrom s n)+    getRandomWords n = RVar (\s -> getRandomWordsFrom s n)++-- some 'fundamental' RVars+-- this maybe ought to even be a part of the RandomSource class...+-- |A random variable evenly distributed over all unsigned integers from+-- 0 to 2^(8*n)-1, inclusive.+nByteInteger :: Int -> RVar Integer+nByteInteger n+    | n .&. 7 == 0+    = do+        xs <- getRandomWords (n `shiftR` 3)+        return $! concatWords xs+    | n > 8+    = do+        let nWords = n `shiftR` 3+            nBytes = n .&. 7+        ws <- getRandomWords nWords+        bs <- getRandomBytes nBytes+        return $! ((concatWords ws `shiftL` (nBytes `shiftL` 3)) .|. concatBytes bs)+    | otherwise+    = do+        xs <- getRandomBytes n+        return $! concatBytes xs++-- |A random variable evenly distributed over all unsigned integers from+-- 0 to 2^n-1, inclusive.+nBitInteger :: Int -> RVar Integer+nBitInteger n+    | n .&. 7 == 0+    = nByteInteger (n `shiftR` 3)+    | otherwise+    = do+        x <- nByteInteger ((n `shiftR` 3) + 1)+        return $! (x .&. (bit n - 1))++concatBytes :: (Bits a, Num a) => [Word8] -> a+concatBytes = concatBits fromIntegral++concatWords :: (Bits a, Num a) => [Word64] -> a+concatWords = concatBits fromIntegral++concatBits :: (Bits a, Bits b, Num b) => (a -> b) -> [a] -> b+concatBits f [] = 0+concatBits f (x:xs) = f x .|. (concatBits f xs `shiftL` bitSize x)
+ src/Data/Random/RVar.hs-boot view
@@ -0,0 +1,16 @@+{-+ -      ``Data/Random/RVar''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances+  #-}++module Data.Random.RVar where++import Data.Random.Source+import {-# SOURCE #-} Data.Random.Distribution++data RVar a+instance MonadRandom RVar+instance Distribution RVar a
+ src/Data/Random/Source.hs view
@@ -0,0 +1,75 @@+{-+ -      ``Data/Random/Source''+ -}+{-# LANGUAGE+    MultiParamTypeClasses, FlexibleInstances+  #-}++module Data.Random.Source+    ( MonadRandom(..)+    , RandomSource(..)+    ) where++import Data.Word+import Data.Bits+import Data.List++import Data.Random.Internal.Words++-- |A typeclass for monads with a chosen source of entropy.  For example,+-- 'RVar' is such a monad - the source from which it is (eventually) sampled+-- is the only source from which a random variable is permitted to draw, so+-- when directly requesting entropy for a random variable these functions+-- are used.+-- +-- The minimal definition is either 'getRandomBytes' or 'getRandomWords'.+class Monad m => MonadRandom m where+    -- |get the specified number of random (uniformly distributed) bytes+    getRandomBytes :: Int -> m [Word8]+    getRandomBytes n+        | n .&. 7 == 0+        = do+            let wc = n `shiftR` 3+            ws <- getRandomWords wc+            return (concatMap wordToBytes ws)+        | otherwise+        = do+            let wc = (n `shiftR` 3) + 1+            ws <- getRandomWords wc+            return . take n . concatMap wordToBytes $ ws+        +    -- |alternate basis function, providing access to larger chunks+    getRandomWords :: Int -> m [Word64]+    getRandomWords n = do+        bs <- getRandomBytes (n `shiftL` 3)+        return (bytesToWords bs)++-- |A source of entropy which can be used in the given monad.+--+-- The minimal definition is either 'getRandomBytesFrom' or 'getRandomWordsFrom'+class Monad m => RandomSource m s where+    getRandomBytesFrom :: s -> Int -> m [Word8]+    getRandomBytesFrom src n+        | n .&. 7 == 0+        = do+            let wc = n `shiftR` 3+            ws <- getRandomWordsFrom src wc+            return (concatMap wordToBytes ws)+        | otherwise+        = do+            let wc = (n `shiftR` 3) + 1+            ws <- getRandomWordsFrom src wc+            return . take n . concatMap wordToBytes $ ws+        +    +    getRandomWordsFrom :: s -> Int -> m [Word64]+    getRandomWordsFrom src n = do+        bs <- getRandomBytesFrom src (n `shiftL` 3)+        return (bytesToWords bs)++instance Monad m => RandomSource m (Int -> m [Word8]) where+    getRandomBytesFrom = id++instance Monad m => RandomSource m (Int -> m [Word64]) where+    getRandomWordsFrom = id+
+ src/Data/Random/Source/DevRandom.hs view
@@ -0,0 +1,24 @@+{-+ -      ``Data/Random/Source/DevRandom''+ -}+{-# LANGUAGE+    MultiParamTypeClasses+  #-}++module Data.Random.Source.DevRandom where++import Data.Random.Source++import GHC.IOBase (unsafePerformIO)+import Data.ByteString (hGet, unpack)+import System.IO (openBinaryFile, IOMode(..))++-- |On systems that have it, \/dev\/random is a handy-dandy ready-to-use source+-- of nonsense.+data DevRandom = DevRandom+{-# NOINLINE devRandom #-}+devRandom = unsafePerformIO (openBinaryFile "/dev/random" ReadMode)++instance RandomSource IO DevRandom where+    getRandomBytesFrom DevRandom n = fmap unpack (hGet devRandom n)+
+ src/Data/Random/Source/PureMT.hs view
@@ -0,0 +1,55 @@+{-+ -      ``Data/Random/Source/PureMT''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleContexts, FlexibleInstances+  #-}++module Data.Random.Source.PureMT where++import Data.Random.Source+import System.Random.Mersenne.Pure64++import Data.StateRef+import Data.Word++import Control.Monad.State++-- |Given a mutable reference to a 'PureMT' generator, we can make a+-- 'RandomSource' usable in any monad in which the reference can be modified.+--+-- For example, if @x :: TVar PureMT@, @getRandomWordsFromMTRef x@ can be+-- used as a 'RandomSource' in 'IO', 'STM', or any monad which is an instance+-- of 'MonadIO'.+getRandomWordsFromMTRef :: ModifyRef sr m PureMT => sr -> Int -> m [Word64]+getRandomWordsFromMTRef ref n = do+    atomicModifyRef ref (randomWords n [])+    +    where+        swap (a,b) = (b,a)+        randomWords    0  ws mt = (mt, ws)+        randomWords (n+1) ws mt = case randomWord64 mt of+            (w, mt) -> randomWords n (w:ws) mt++-- |Similarly, @getRandomWordsFromMTState x@ can be used in any \"state\"+-- monad in the mtl sense whose state is a 'PureMT' generator.+-- Additionally, the standard mtl state monads have 'MonadRandom' instances+-- which do precisely that, allowing an easy conversion of 'RVar's and+-- other 'Distribution' instances to \"pure\" random variables.+getRandomWordsFromMTState :: MonadState PureMT m => Int -> m [Word64]+getRandomWordsFromMTState n = do+    mt <- get+    let randomWords    0  ws mt = (mt, ws)+        randomWords (n+1) ws mt = case randomWord64 mt of+            (w, mt) -> randomWords n (w:ws) mt+        +        (newMt, ws) = randomWords n [] mt+    put newMt+    return ws++instance MonadRandom (State PureMT) where+    getRandomWords = getRandomWordsFromMTState++instance Monad m => MonadRandom (StateT PureMT m) where+    getRandomWords = getRandomWordsFromMTState
+ src/Data/Random/Source/Std.hs view
@@ -0,0 +1,20 @@+{-+ -      ``Data/Random/Source/Std''+ -}+{-# LANGUAGE+    MultiParamTypeClasses, FlexibleInstances+  #-}++module Data.Random.Source.Std where++import Data.Random.Source++-- |A token representing the \"standard\" entropy source in a 'MonadRandom'+-- monad.  Its sole purpose is to make the following true (when the types check):+--+-- > sampleFrom StdRandom === sample+data StdRandom = StdRandom++instance MonadRandom m => RandomSource m StdRandom where+    getRandomBytesFrom StdRandom = getRandomBytes+    getRandomWordsFrom StdRandom = getRandomWords
+ src/Data/Random/Source/StdGen.hs view
@@ -0,0 +1,88 @@+{-+ -      ``Data/Random/Source/StdGen''+ -}+{-# LANGUAGE+    MultiParamTypeClasses, FlexibleInstances, UndecidableInstances+  #-}++module Data.Random.Source.StdGen where++import Data.Random.Source+import System.Random+import Control.Monad+import Control.Monad.State+import Data.StateRef+import Data.Word++instance (ModifyRef (IORef   StdGen) m StdGen) => RandomSource m (IORef   StdGen) where+    getRandomBytesFrom = getRandomBytesFromRandomGenRef+    getRandomWordsFrom = getRandomWordsFromRandomGenRef+instance (ModifyRef (TVar    StdGen) m StdGen) => RandomSource m (TVar    StdGen) where+    getRandomBytesFrom = getRandomBytesFromRandomGenRef+    getRandomWordsFrom = getRandomWordsFromRandomGenRef+instance (ModifyRef (STRef s StdGen) m StdGen) => RandomSource m (STRef s StdGen) where+    getRandomBytesFrom = getRandomBytesFromRandomGenRef+    getRandomWordsFrom = getRandomWordsFromRandomGenRef++getRandomBytesFromStdGenIO :: Int -> IO [Word8]+getRandomBytesFromStdGenIO n = do+    ints <- replicateM n (randomRIO (0, 255))+    let bytes = map fromIntegral (ints :: [Int])+    return bytes++-- |Given a mutable reference to a 'RandomGen' generator, we can make a+-- 'RandomSource' usable in any monad in which the reference can be modified.+--+-- For example, if @x :: TVar StdGen@, @getRandomBytesFromRandomGenRef x@ can be+-- used as a 'RandomSource' in 'IO', 'STM', or any monad which is an instance+-- of 'MonadIO'.  It's generally probably better to use+-- 'getRandomWordsFromRandomGenRef' though, as this one is likely to throw+-- away a lot of perfectly good entropy.+getRandomBytesFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>+                                  sr -> Int -> m [Word8]+getRandomBytesFromRandomGenRef g n = do+    let swap (a,b) = (b,a)+    ints <- replicateM n (atomicModifyRef g (swap . randomR (0, 255)))+    let bytes = map fromIntegral (ints :: [Int])+    return bytes+    +-- |Similarly, @getRandomWordsFromRandomGenState x@ can be used in any \"state\"+-- monad in the mtl sense whose state is a 'RandomGen' generator.+-- Additionally, the standard mtl state monads have 'MonadRandom' instances+-- which do precisely that, allowing an easy conversion of 'RVar's and+-- other 'Distribution' instances to \"pure\" random variables.+getRandomBytesFromRandomGenState :: (RandomGen g, MonadState g m) =>+                                  Int -> m [Word8]+getRandomBytesFromRandomGenState n = replicateM n $ do+    g <- get+    case randomR (0,255 :: Int) g of+        (i,g) -> do+            put g+            return (fromIntegral i)++-- |See 'getRandomBytesFromRandomGenRef'+getRandomWordsFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>+                                  sr -> Int -> m [Word64]+getRandomWordsFromRandomGenRef g n = do+    let swap (a,b) = (b,a)+    ints <- replicateM n (atomicModifyRef g (swap . randomR (0, 2^64-1)))+    let bytes = map fromInteger ints+    return bytes+    +-- |See 'getRandomBytesFromRandomGenState'+getRandomWordsFromRandomGenState :: (RandomGen g, MonadState g m) =>+                                  Int -> m [Word64]+getRandomWordsFromRandomGenState n = replicateM n $ do+    g <- get+    case randomR (0,2^64-1) g of+        (i,g) -> do+            put g+            return (fromInteger i)++instance MonadRandom (State StdGen) where+    getRandomBytes = getRandomBytesFromRandomGenState+    getRandomWords = getRandomWordsFromRandomGenState++instance Monad m => MonadRandom (StateT StdGen m) where+    getRandomBytes = getRandomBytesFromRandomGenState+    getRandomWords = getRandomWordsFromRandomGenState