random-fu 0.0.3.2 → 0.1.0.0

raw patch · 31 files changed

+1620/−836 lines, 31 filesdep +MonadPromptdep +erf-nativedep +mwc-randomdep −storablevectorPVP ok

version bump matches the API change (PVP)

Dependencies added: MonadPrompt, erf-native, mwc-random, tagged, vector

Dependencies removed: storablevector

API changes (from Hackage documentation)

Files

• random-fu.cabal +52/−12
• src/Data/Random.hs +64/−52
• src/Data/Random/Distribution.hs +60/−17
• src/Data/Random/Distribution/Bernoulli.hs +5/−7
• src/Data/Random/Distribution/Beta.hs +12/−11
• src/Data/Random/Distribution/Binomial.hs +27/−32
• src/Data/Random/Distribution/Categorical.hs +168/−0
• src/Data/Random/Distribution/Dirichlet.hs +30/−0
• src/Data/Random/Distribution/Discrete.hs +0/−136
• src/Data/Random/Distribution/Gamma.hs +39/−77
• src/Data/Random/Distribution/Multinomial.hs +31/−0
• src/Data/Random/Distribution/Normal.hs +28/−28
• src/Data/Random/Distribution/Poisson.hs +5/−7
• src/Data/Random/Distribution/Triangular.hs +13/−13
• src/Data/Random/Distribution/Uniform.hs +88/−42
• src/Data/Random/Distribution/Weibull.hs +16/−0
• src/Data/Random/Distribution/Ziggurat.hs +80/−56
• src/Data/Random/Internal/Find.hs +18/−14
• src/Data/Random/Internal/Primitives.hs +235/−0
• src/Data/Random/Internal/TH.hs +13/−10
• src/Data/Random/Internal/Words.hs +78/−6
• src/Data/Random/Lift.hs +6/−0
• src/Data/Random/List.hs +15/−11
• src/Data/Random/RVar.hs +141/−55
• src/Data/Random/Sample.hs +14/−2
• src/Data/Random/Source.hs +109/−56
• src/Data/Random/Source/DevRandom.hs +17/−4
• src/Data/Random/Source/MWC.hs +46/−0
• src/Data/Random/Source/PureMT.hs +72/−89
• src/Data/Random/Source/Std.hs +9/−3
• src/Data/Random/Source/StdGen.hs +129/−96

@@ -1,6 +1,6 @@ name:                   random-fu-version:                0.0.3.2-stability:              experimental+version:                0.1.0.0+stability:              provisional  cabal-version:          >= 1.2 build-type:             Simple@@ -12,34 +12,65 @@  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.+description:            Random number generation based on modeling random +                        variables in two complementary ways: first, by the+                        parameters of standard mathematical distributions and,+                        second, by an abstract type ('RVar') which can be+                        composed and manipulated monadically and sampled in+                        either monadic or \"pure\" styles.+                        +                        The primary purpose of this library is to support +                        defining and sampling a wide variety of high quality+                        random variables.  Quality is prioritized over speed,+                        but performance is an important goal too.+                        +                        In my testing, I have found it capable of speed +                        comparable to other Haskell libraries, but still+                        a fair bit slower than straight C implementations of +                        the same algorithms.+                        +                        Warning to anyone upgrading from \"< 0.1\": 'Discrete'+                        has been renamed 'Categorical', the entropy source +                        classes have been redesigned, and many things are no+                        longer exported from the root module "Data.Random"+                        (In particular, DevRandom - this is not available on +                        windows, so it will likely move to its own package +                        eventually so that client code dependencies on it will +                        be made explicit).+                        +                        The "Data.Random" module itself should now have a+                        relatively stable interface, but the other modules+                        are still subject to change.  Specifically, I am +                        considering hiding data constructors for most or all +                        of the distributions.  Flag base4 Flag base4_2     Description:        base-4.2 has an incompatible change in Data.Fixed (HasResolution)  Library+  ghc-options:          -Wall -funbox-strict-fields -fno-method-sharing   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.Discrete+                        Data.Random.Distribution.Categorical+                        Data.Random.Distribution.Dirichlet                         Data.Random.Distribution.Exponential                         Data.Random.Distribution.Gamma+                        Data.Random.Distribution.Multinomial                         Data.Random.Distribution.Normal                         Data.Random.Distribution.Poisson                         Data.Random.Distribution.Rayleigh                         Data.Random.Distribution.Triangular                         Data.Random.Distribution.Uniform+                        Data.Random.Distribution.Weibull                         Data.Random.Distribution.Ziggurat                         Data.Random.Internal.Find                         Data.Random.Internal.Fixed+                        Data.Random.Internal.Primitives                         Data.Random.Internal.TH                         Data.Random.Internal.Words                         Data.Random.Lift@@ -47,7 +78,7 @@                         Data.Random.RVar                         Data.Random.Sample                         Data.Random.Source-                        Data.Random.Source.DevRandom+                        Data.Random.Source.MWC                         Data.Random.Source.PureMT                         Data.Random.Source.Std                         Data.Random.Source.StdGen@@ -64,12 +95,21 @@        build-depends:        array,                         containers,-                        erf,                         mersenne-random-pure64,                         monad-loops >= 0.3.0.1,+                        MonadPrompt,+                        mwc-random,                         mtl,                         random,                         random-shuffle,                         stateref >= 0.3 && < 0.4,-                        storablevector,-                        template-haskell+                        tagged,+                        template-haskell,+                        vector++  if os(Windows)+    cpp-options:        -Dwindows+    build-depends:      erf-native+  else+    build-depends:      erf+    exposed-modules:    Data.Random.Source.DevRandom

@@ -1,66 +1,78 @@-{-- -      ``Data/Random''- -}-{-# LANGUAGE-    FlexibleContexts-  #-}---- |Random numbers and stuff...+-- |Flexible modeling and sampling of random variables.+--+-- The central abstraction in this library is the concept of a random +-- variable.  It is not fully formalized in the standard measure-theoretic +-- language, but rather is informally defined as a \"thing you can get random +-- values out of\".  Different random variables may have different types of +-- values they can return or the same types but different probabilities for+-- each value they can return.  The random values you get out of them are+-- traditionally called \"random variates\". -- --- "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.+-- Most imperative-language random number libraries are all about obtaining +-- and manipulating random variates.  This one is about defining, manipulating +-- and sampling random variables.  Computationally, the distinction is small +-- and mostly just a matter of perspective, but from a program design +-- perspective it provides both a powerfully composable abstraction and a+-- very useful separation of concerns. -- --- "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.-+-- Abstract random variables as implemented by 'RVar' are composable.  They can+-- be defined in a monadic / \"imperative\" style that amounts to manipulating+-- variates, but with strict type-level isolation.  Concrete random variables+-- are also provided, but they do not compose as generically.  The 'Distribution'+-- type class allows concrete random variables to \"forget\" their concreteness +-- so that they can be composed.  For examples of both, see the documentation +-- for 'RVar' and 'Distribution', as well as the code for any of the concrete +-- distributions such as 'Uniform', 'Gamma', etc.+-- +-- Both abstract and concrete random variables can be sampled (despite the+-- types GHCi may list for the functions) by the functions in "Data.Random.Sample".+-- +-- Random variable sampling is done with regard to a generic basis of primitive+-- random variables defined in "Data.Random.Internal.Primitives".  This basis +-- is very low-level and the actual set of primitives is still fairly experimental,+-- which is why it is in the \"Internal\" sub-heirarchy.  User-defined variables+-- should use the existing high-level variables such as 'Uniform' and 'Normal'+-- rather than these basis variables.  "Data.Random.Source" defines classes for+-- entropy sources that provide implementations of these primitive variables. +-- Several implementations are available in the Data.Random.Source.* modules. module Data.Random-    ( module Data.Random.Sample-    , 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.Rayleigh-    , module Data.Random.Distribution.Triangular-    , module Data.Random.Distribution.Uniform-    , module Data.Random.Distribution.Ziggurat-    , module Data.Random.List-    , module Data.Random.RVar+    ( -- * Random variables+      -- ** Abstract ('RVar')+      RVar, RVarT,+      runRVar, runRVarT, runRVarTWith,++      -- ** Concrete ('Distribution')+      Distribution(..), CDF(..),+      +      -- * Sampling random variables+      Sampleable(..), sample, sampleState, sampleStateT,+      +      -- * A few very common distributions+      Uniform(..), uniform, +      StdUniform(..), stdUniform,+      Normal(..), normal, stdNormal,+      Gamma(..), gamma,+      +      -- * Entropy Sources+      MonadRandom, RandomSource, StdRandom(..),+      +      -- * Useful list-based operations+      randomElement,+      shuffle, shuffleN, shuffleNofM+           ) where  import Data.Random.Sample-import Data.Random.Source-import Data.Random.Source.DevRandom-import Data.Random.Source.StdGen-import Data.Random.Source.PureMT+import Data.Random.Source (MonadRandom, RandomSource)+import Data.Random.Source.MWC ()+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.Rayleigh-import Data.Random.Distribution.Triangular import Data.Random.Distribution.Uniform-import Data.Random.Distribution.Ziggurat+ import Data.Random.Lift () import Data.Random.List import Data.Random.RVar

@@ -1,22 +1,66 @@-{-- -      ``Data/Random/Distribution''- -} {-# LANGUAGE     MultiParamTypeClasses, FlexibleContexts   #-}- module Data.Random.Distribution where  import Data.Random.Lift import Data.Random.RVar --- |A definition of a random variable's distribution.  From the distribution--- an 'RVar' can be created, or the distribution can be directly sampled using --- 'sampleFrom' or 'sample'.+-- |A 'Distribution' is a data representation of a random variable's probability+-- structure.  For example, in "Data.Random.Distribution.Normal", the 'Normal'+-- distribution is defined as:+--+-- > data Normal a+-- >     = StdNormal+-- >     | Normal a a+-- +-- Where the two parameters of the 'Normal' data constructor are the mean and+-- standard deviation of the random variable, respectively.  To make use of+-- the 'Normal' type, one can convert it to an 'rvar' and manipulate it or+-- sample it directly:+--+-- > x <- sample (rvar (Normal 10 2))+-- > x <- sample (Normal 10 2)+-- +-- A 'Distribution' is typically more transparent than an 'RVar'+-- but less composable (precisely because of that transparency).  There are +-- several practical uses for types implementing 'Distribution':+-- +-- * Typically, a 'Distribution' will expose several parameters of a standard +-- mathematical model of a probability distribution, such as mean and std deviation for+-- the normal distribution.  Thus, they can be manipulated analytically using+-- mathematical insights about the distributions they represent.  For example,+-- a collection of bernoulli variables could be simplified into a (hopefully) smaller+-- collection of binomial variables.+-- +-- * Because they are generally just containers for parameters, they can be+-- easily serialized to persistent storage or read from user-supplied +-- configurations (eg, initialization data for a simulation).+-- +-- * If a type additionally implements the 'CDF' subclass, which extends+-- 'Distribution' with a cumulative density function, an arbitrary random+-- variable 'x' can be tested against the distribution by testing+-- @fmap (cdf dist) x@ for uniformity.+-- +-- On the other hand, most 'Distribution's will not be closed under all the+-- same operations as 'RVar' (which, being a monad, has a fully turing-complete+-- internal computational model).  The sum of two uniformly-distributed +-- variables, for example, is not uniformly distributed.  To support general +-- composition, the 'Distribution' class defines a function 'rvar' to +-- construct the more-abstract and more-composable 'RVar' representation +-- of a random variable. class Distribution d t where     -- |Return a random variable with this distribution.     rvar :: d t -> RVar t+    rvar = rvarT+    +    -- |Return a random variable with the given distribution, pre-lifted to an arbitrary 'RVarT'.+    -- Any arbitrary 'RVar' can also be converted to an 'RVarT m' for an arbitrary 'm', using+    -- either 'lift' or 'sample'.+    rvarT :: d t -> RVarT n t+    rvarT d = lift (rvar d) + class Distribution d t => CDF d t where     -- |Return the cumulative distribution function of this distribution.     -- That is, a function taking @x :: t@ to the probability that the next@@ -31,16 +75,15 @@     -- must be uniformly distributed over (0,1).  Inclusion of either endpoint is optional,     -- though the preferred range is (0,1].     -- -    -- Thus, 'cdf' for a product type should not be a joint CDF as commonly -    -- defined, as that definition violates both conditions.+    -- Note that this definition requires that  'cdf' for a product type +    -- should _not_ be a joint CDF as commonly defined, as that definition +    -- violates both conditions.     -- Instead, it should be a univariate CDF over the product type.  That is,     -- it should represent the CDF with respect to the lexicographic order-    -- of the tuple.+    -- of the product.+    -- +    -- The present specification is probably only really useful for testing+    -- conformance of a variable to its target distribution, and I am open to+    -- suggestions for more-useful specifications (especially with regard to+    -- the interaction with product types).     cdf :: d t -> t -> Double---- |Return a random variable with the given distribution, pre-lifted to an arbitrary 'RVarT'.--- Any arbitrary 'RVar' can also be converted to an 'RVarT m' for an arbitrary 'm', using--- either 'lift' or 'sample'.-rvarT :: Distribution d t => d t -> RVarT n t-rvarT d = lift (rvar d)-

@@ -16,8 +16,6 @@ import Data.Random.Distribution import Data.Random.Distribution.Uniform -import Data.Int-import Data.Word import Data.Ratio import Data.Complex @@ -35,7 +33,7 @@     return (x <= p)  boolBernoulliCDF :: (Real a) => a -> Bool -> Double-boolBernoulliCDF p True  = 1+boolBernoulliCDF _ True  = 1 boolBernoulliCDF p False = (1 - realToFrac p)  -- | @generalBernoulli t f p@ generates a random variable whose value is @t@@@ -46,10 +44,10 @@     return (if x then t else f)  generalBernoulliCDF :: CDF (Bernoulli b) Bool => (a -> a -> Bool) -> a -> a -> b -> a -> Double-generalBernoulliCDF (>=) f t p x-    | f >= t    = error "generalBernoulliCDF: f >= t"-    | x >= t    = cdf (Bernoulli p) True-    | x >= f    = cdf (Bernoulli p) False+generalBernoulliCDF gte f t p x+    | f `gte` t = error "generalBernoulliCDF: f >= t"+    | x `gte` t = cdf (Bernoulli p) True+    | x `gte` f = cdf (Bernoulli p) False     | otherwise = 0  data Bernoulli b a = Bernoulli b

@@ -4,18 +4,21 @@ {-# LANGUAGE     MultiParamTypeClasses,     FlexibleInstances, FlexibleContexts,-    UndecidableInstances+    UndecidableInstances,+    TemplateHaskell   #-}  module Data.Random.Distribution.Beta where +import Data.Random.Internal.TH+ import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Gamma import Data.Random.Distribution.Uniform -import Control.Monad-+{-# SPECIALIZE fractionalBeta :: Float  -> Float  -> RVar Float #-}+{-# SPECIALIZE fractionalBeta :: Double -> Double -> RVar Double #-} fractionalBeta :: (Fractional a, Distribution Gamma a, Distribution StdUniform a) => a -> a -> RVar a fractionalBeta 1 1 = stdUniform fractionalBeta a b = do@@ -23,16 +26,14 @@     y <- gamma b 1     return (x / (x + y)) -fractionalBetaFromIntegral :: (Fractional c, Distribution (Erlang a) c, Distribution (Erlang b) c) => a -> b -> RVar c-fractionalBetaFromIntegral a b =  do-    x <- erlang a-    y <- erlang b-    return (x / (x + y))-+{-# SPECIALIZE beta :: Float  -> Float  -> RVar Float #-}+{-# SPECIALIZE beta :: Double -> Double -> RVar Double #-} beta :: Distribution Beta a => a -> a -> RVar a beta a b = rvar (Beta a b)  data Beta a = Beta a a -instance (Fractional a, Distribution Gamma a, Distribution StdUniform a) => Distribution Beta a where-    rvar (Beta a b) = fractionalBeta a b+$( replicateInstances ''Float realFloatTypes [d|+        instance Distribution Beta Float+              where rvar (Beta a b) = fractionalBeta a b+    |])

@@ -4,7 +4,8 @@ {-# LANGUAGE     MultiParamTypeClasses,     FlexibleInstances, FlexibleContexts,-    UndecidableInstances, TemplateHaskell+    UndecidableInstances, TemplateHaskell,+    BangPatterns   #-}  module Data.Random.Distribution.Binomial where@@ -16,29 +17,19 @@ 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     -- note that although it's fast enough for large (eg, 2^10000)      -- @Integer@s, it's not accurate enough when using @Double@ as     -- the @b@ parameter.+{-# SPECIALIZE integralBinomial :: Int -> Float  -> RVar Int #-}+{-# SPECIALIZE integralBinomial :: Int -> Double -> RVar Int #-}+{-# SPECIALIZE integralBinomial :: Integer -> Float  -> RVar Integer #-}+{-# SPECIALIZE integralBinomial :: Integer -> Double -> RVar Integer #-} integralBinomial :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> b -> RVar a-integralBinomial t p = bin 0 t p+integralBinomial = bin 0     where-        -- GHC likes to discharge Beta to the Beta instance's context, which-        -- @integralBinomial@'s context doesn't (directly) satisfy.-        -- Seems like GHC could do better.  GHC could discharge it in @integralBinomial@'s-        -- context when attempting to satisfy bin, since the Beta instance covers all @b@,-        -- whereupon it would find that the contexts do, in fact, match.-        -- Of course, it's a pretty obscure case, so maybe it's better to just -        -- leave it to the coder who's doing weird stuff to tell the compiler what-        -- he/she wants.-        -- Anyway, this type signature makes GHC happy.-        bin :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> a -> b -> RVar a-        bin k t p+        bin !k !t !p             | t > 10    = do                 let a = 1 + t `div` 2                     b = 1 + t - a@@ -50,24 +41,16 @@                      | otherwise = count k t                 where-                    count k  0    = return k-                    count k (n+1) = do+                    count !k'  0    = return k'+                    count !k' (n+1) = do                         x <- stdUniform-                        (count $! (if x < p then k + 1 else k)) n-        -        -- this gives 'believable' results... but is it really any good?-        selectA cutoff t = return a---             | cutoff >= a   = return a---             | otherwise     = generalBernoulli a (a `div` cutoff) p---             -             where-                 a = 1 + t `div` 2---                 p = 0.5 :: Float -- fromIntegral cutoff / fromIntegral a :: Float-                +                        count (if x < p then k' + 1 else k') n+                    count _ _ = error "integralBinomial: negative number of trials specified"+ -- TODO: improve performance integralBinomialCDF :: (Integral a, Real b) => a -> b -> a -> Double-integralBinomialCDF n p x = sum-    [ fromIntegral (n `c` i) * p' ^^ i * (1-p') ^^ (n-i)+integralBinomialCDF t p x = sum+    [ fromIntegral (t `c` i) * p' ^^ i * (1-p') ^^ (t-i)     | i <- [0 .. x]     ]     @@ -77,12 +60,24 @@  -- would it be valid to repeat the above computation using fractional @t@? -- obviously something different would have to be done with @count@ as well...+{-# SPECIALIZE floatingBinomial :: Float  -> Float  -> RVar Float  #-}+{-# SPECIALIZE floatingBinomial :: Float  -> Double -> RVar Float  #-}+{-# SPECIALIZE floatingBinomial :: Double -> Float  -> RVar Double #-}+{-# SPECIALIZE floatingBinomial :: Double -> Double -> RVar Double #-} floatingBinomial :: (RealFrac a, Distribution (Binomial b) Integer) => a -> b -> RVar a floatingBinomial t p = fmap fromInteger (rvar (Binomial (truncate t) p))  floatingBinomialCDF :: (CDF (Binomial b) Integer, RealFrac a) => a -> b -> a -> Double floatingBinomialCDF t p x = cdf (Binomial (truncate t :: Integer) p) (floor x) +{-# SPECIALIZE binomial :: Int -> Float  -> RVar Int #-}+{-# SPECIALIZE binomial :: Int -> Double -> RVar Int #-}+{-# SPECIALIZE binomial :: Integer -> Float  -> RVar Integer #-}+{-# SPECIALIZE binomial :: Integer -> Double -> RVar Integer #-}+{-# SPECIALIZE binomial :: Float  -> Float  -> RVar Float  #-}+{-# SPECIALIZE binomial :: Float  -> Double -> RVar Float  #-}+{-# SPECIALIZE binomial :: Double -> Float  -> RVar Double #-}+{-# SPECIALIZE binomial :: Double -> Double -> RVar Double #-} binomial :: Distribution (Binomial b) a => a -> b -> RVar a binomial t p = rvar (Binomial t p) 

@@ -0,0 +1,168 @@+{-+ -      ``Data/Random/Distribution/Categorical''+ -}+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts+  #-}++module Data.Random.Distribution.Categorical where++import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform++import Control.Arrow+import Control.Monad+import Control.Applicative+import Data.Foldable (Foldable(foldMap))+import Data.Traversable (Traversable(traverse, sequenceA))++import Data.List+import Data.Function++-- |Construct a 'Categorical' distribution from a list of probabilities+-- and categories, where the probabilities all sum to 1.+categorical :: Distribution (Categorical p) a => [(p,a)] -> RVar a+categorical ps = rvar (Categorical ps)++-- | Construct a 'Categorical' distribution from a list of weighted categories,+-- where the weights do not necessarily sum to 1.+{-# INLINE weightedCategorical #-}+weightedCategorical :: (Fractional p) => [(p,a)] -> Categorical p a+weightedCategorical = normalizeCategoricalPs . Categorical++-- |Construct a 'Categorical' distribution from a list of observed outcomes.+-- Equivalent events will be grouped and counted, and the probabilities of each+-- event in the returned distribution will be proportional to the number of +-- occurrences of that event.+empirical :: (Fractional p, Ord a) => [a] -> Categorical p a+empirical xs = normalizeCategoricalPs (Categorical bins)+    where bins = [ (genericLength bin, x)+                 | bin@(x:_) <- group (sort xs)+                 ]++-- |Categorical distribution; a list of events with corresponding probabilities.+-- The sum of the probabilities must be 1, and no event should have a zero +-- or negative probability (at least, at time of sampling; very clever users+-- can do what they want with the numbers before sampling, just make sure +-- that if you're one of those clever ones, you normalize before sampling).+newtype Categorical p a = Categorical [(p, a)]+    deriving (Eq, Show)++instance (Fractional p, Ord p, Distribution StdUniform p) => Distribution (Categorical p) a where+    rvar (Categorical []) = fail "categorical distribution over empty set cannot be sampled"+    rvar (Categorical ds) = do+        let (ps, xs) = unzip ds+            cs = scanl1 (+) ps+        +        u <- stdUniform+        getEvent u cs xs+        +        where+            -- In the (hopefully) extremely rare event that, due to numerical+            -- instability, the last 'c' is less than 1 _and_ a number greater than +            -- it is drawn, simply retry the sampling.  If it comes to that, also+            -- do one last sanity check that lastC > 0, to make sure that there+            -- is some nonzero chance of termination.+            getEvent u cs0 xs0 = go 0 cs0 xs0+                where+                    go lastC [] _+                        | lastC > 0 = do {newU <- stdUniform; getEvent newU cs0 xs0}+                        | otherwise = fail "categorical distribution sampling error: total probablility not greater than zero"+                    go lastC (c:cs) (x:xs)+                        | c < lastC = fail "categorical distribution sampling error: negative probability for an event!"+                        | u > c     = go c cs xs+                        | c == c    = return x+                        | otherwise = fail "categorical distribution sampling error: NaN probability"+                    +                    go _ _ _ = error "rvar/Categorical: programming error! this case should be impossible!"++instance Functor (Categorical p) where+    fmap f (Categorical ds) = Categorical [(p, f x) | ~(p, x) <- ds]++instance Foldable (Categorical p) where+    foldMap f (Categorical ds) = foldMap (f . snd) ds++instance Traversable (Categorical p) where+    traverse f (Categorical ds) = Categorical <$> traverse (\(p,e) -> (\e' -> (p,e')) <$> f e) ds+    sequenceA  (Categorical ds) = Categorical <$> traverse (\(p,e) -> (\e' -> (p,e')) <$>   e) ds++instance Fractional p => Monad (Categorical p) where+    return x = Categorical [(1, x)]+    +    -- I'm not entirely sure whether this is a valid form of failure; see next+    -- set of comments.+    fail _ = Categorical []+    +    -- Should the normalize step be included here, or should normalization+    -- be assumed?  It seems like there is (at least) 1 valid situation where+    -- non-normal results would arise:  the distribution being modeled is +    -- "conditional" and some event arose that contradicted the assumed +    -- condition and thus was eliminated ('f' returned an empty or +    -- zero-probability consequent, possibly by 'fail'ing).+    -- +    -- It seems reasonable to continue in such circumstances, but should there+    -- be any renormalization?  If so, does it make a difference when that +    -- renormalization is done?  I'm pretty sure it does, actually.  So, the+    -- normalization will be omitted here for now, as it's easier for the+    -- user (who really better know what they mean if they're returning+    -- non-normalized probability anyway) to normalize explicitly than to+    -- undo any normalization that was done automatically.+    (Categorical xs) >>= f = {- normalizeCategoricalPs . -} Categorical $ do+        (p, x) <- xs+        +        let Categorical fx = f x+        (q, y) <- fx+        +        return (p * q, y)++instance Fractional p => Applicative (Categorical p) where+    pure = return+    (<*>) = ap++-- |Like 'fmap', but for the probabilities of a categorical distribution.+mapCategoricalPs :: (p -> q) -> Categorical p e -> Categorical q e+mapCategoricalPs f (Categorical ds) = Categorical [(f p, x) | (p, x) <- ds]++-- |Adjust all the weights of a categorical distribution so that they +-- sum to unity.+normalizeCategoricalPs :: (Fractional p) => Categorical p e -> Categorical p e+normalizeCategoricalPs orig@(Categorical ds) = +    -- For practical purposes the scale factor is strict anyway,+    -- so check if the total probability is 1 and, if so, skip +    -- the actual scaling part.+    --+    -- Along the way, discard any zero-probability events.+    if null ds || ps =~ 1+        then orig+        else Categorical+                [ (p * scale, e)+                | (p, e) <- ds+                , p /= 0+                ] +    where+        ps = foldl1' (+) (map fst ds)+        scale = recip ps+        +        -- Using same implicit-epsilon trick as in Distribution instance+        -- (see comments there)+        x =~ y  = (100 + (x-y) == 100)+++-- |Simplify a categorical distribution by combining equivalent categories (the new+-- category will have a probability equal to the sum of all the originals).+collectEvents :: (Ord e, Num p, Ord p) => Categorical p e -> Categorical p e+collectEvents = collectEventsBy compare ((sum *** head) . unzip)+        +-- |Simplify a categorical distribution by combining equivalent events (the new+-- event will have a weight equal to the sum of all the originals).+-- The comparator function is used to identify events to combine.  Once chosen,+-- the events and their weights are combined by the provided probability and+-- event aggregation function.+collectEventsBy :: (e -> e -> Ordering) -> ([(p,e)] -> (p,e))-> Categorical p e -> Categorical p e+collectEventsBy compareE combine (Categorical ds) = +    Categorical . map combine . groupEvents . sortEvents $ ds+    where+        groupEvents = groupBy (\x y -> snd x `compareE` snd y == EQ)+        sortEvents  = sortBy (compareE `on` snd)

@@ -0,0 +1,30 @@+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleInstances, FlexibleContexts,+    UndecidableInstances, GADTs+  #-}++module Data.Random.Distribution.Dirichlet where++import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Gamma++import Data.List++fractionalDirichlet :: (Fractional a, Distribution Gamma a) => [a] -> RVar [a]+fractionalDirichlet []  = return []+fractionalDirichlet [_] = return [1]+fractionalDirichlet as = do+    xs <- sequence [gamma a 1 | a <- as]+    let total = foldl1' (+) xs+    +    return (map (* recip total) xs)++dirichlet :: Distribution Dirichlet [a] => [a] -> RVar [a]+dirichlet as = rvar (Dirichlet as)++newtype Dirichlet a = Dirichlet a deriving (Eq, Show)++instance (Fractional a, Distribution Gamma a) => Distribution Dirichlet [a] where+    rvar (Dirichlet as) = fractionalDirichlet as

@@ -1,136 +0,0 @@-{-- -      ``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 Data.Random.List (randomElement)--import Control.Arrow-import Control.Monad-import Control.Applicative-import Data.Foldable (Foldable(foldMap))-import Data.Traversable (Traversable(traverse, sequenceA))--import Data.List-import Data.Function--discrete :: Distribution (Discrete p) a => [(p,a)] -> RVar a-discrete ps = rvar (Discrete ps)--empirical :: (Num p, Ord a) => [a] -> Discrete p a-empirical xs = Discrete bins-    where bins = [ (genericLength bin, x)-                 | bin@(x:_) <- group (sort xs)-                 ]--newtype Discrete p a = Discrete [(p, a)]-    deriving (Eq, Show)--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"-        -        let totalWeight = last cs-        if  totalWeight <= 0-            -- if all events are "equally impossible", just pick an arbitrary one.-            then randomElement xs-            else do-                let getU = do-                        u <- uniform 0 totalWeight-                        -- reject 0; this causes integral weights to behave sensibly (although it -                        -- potentially wastes up to 50% of all sampled integral values in the-                        -- case where totalWeight = 1) and is still valid for fractional weights-                        -- (it only prevents zero-probability events from ever occurring, which -                        -- is reasonable).-                        if u == 0 then getU-                                  else return u-                u <- getU-                return $ head-                    [ x-                    | (c,x) <- zip cs xs-                    , c >= u-                    ]--instance Functor (Discrete p) where-    fmap f (Discrete ds) = Discrete [(p, f x) | (p, x) <- ds]--instance Foldable (Discrete p) where-    foldMap f (Discrete ds) = foldMap (f . snd) ds--instance Traversable (Discrete p) where-    traverse f (Discrete ds) = Discrete <$> traverse (\(p,e) -> (\e -> (p,e)) <$> f e) ds-    sequenceA  (Discrete ds) = Discrete <$> traverse (\(p,e) -> (\e -> (p,e)) <$>   e) ds---- We want each subset of cases in fx derived from a given case --- in x to have the same relative weight as the set in x from whence they came.-instance Num p => Monad (Discrete p) where-    return x = Discrete [(1, x)]-    (Discrete x) >>= f = Discrete $ do-        (p, x) <- x-        -        let Discrete fx = f x-        (q, x) <- fx-        -        return (p * q, x)--instance Num p => Applicative (Discrete p) where-    pure = return-    (<*>) = ap---- |Like 'fmap', but for the weights of a discrete distribution.-mapDiscreteWeights :: (p -> q) -> Discrete p e -> Discrete q e-mapDiscreteWeights f (Discrete ds) = Discrete [(f p, x) | (p, x) <- ds]---- |Adjust all the weights of a discrete distribution so that they --- sum to unity.  If not possible, returns the original distribution --- unchanged.-normalizeDiscreteWeights :: (Fractional p) => Discrete p e -> Discrete p e-normalizeDiscreteWeights orig@(Discrete ds) = -    -- For practical purposes the scale factor is strict anyway,-    -- so check if it's 0 or 1 and, if so, skip the actual scaling part.-    if ws `elem` [0,1]-        then orig-        else Discrete-                [ (w * scale, e)-                | (w, e) <- ds-                ] -    where-        ws = sum (map fst ds)-        scale = recip ws---- |Simplify a discrete distribution by combining equivalent events (the new--- event will have a weight equal to the sum of all the originals).-collectDiscreteEvents :: (Ord e, Num p, Ord p) => Discrete p e -> Discrete p e-collectDiscreteEvents = collectDiscreteEventsBy compare sum head-        --- |Simplify a discrete distribution by combining equivalent events (the new--- event will have a weight equal to the sum of all the originals).--- The comparator function is used to identify events to combine.  Once chosen,--- the events and their weights are combined (independently) by the provided--- weight and event aggregation functions.-collectDiscreteEventsBy :: (e -> e -> Ordering) -> ([p] -> p) -> ([e] -> e)-> Discrete p e -> Discrete p e-collectDiscreteEventsBy compareE sumWeights mergeEvents (Discrete ds) = -    Discrete . map ((sumWeights *** mergeEvents) . unzip) . groupEvents . sortEvents $ ds-    -    where-        groupEvents = groupBy (\x y -> snd x `compareE` snd y == EQ)-        sortEvents  = sortBy (compareE `on` snd)-        -        weight (p,x)-            | p < 0     = error "negative probability in discrete distribution"-            | otherwise = p-        event ((p,x):_) = x-        -        combine (ps, xs) = (sumWeights ps, mergeEvents xs)

@@ -1,17 +1,7 @@-{-- -      ``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+    UndecidableInstances, BangPatterns   #-}  module Data.Random.Distribution.Gamma where@@ -21,74 +11,44 @@ import Data.Random.Distribution.Uniform import Data.Random.Distribution.Normal -import Control.Monad+import Data.Ratio -    -- 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 :: (Floating a, Ord a, Distribution Normal a, Distribution StdUniform a) => a -> a -> RVar a-realFloatGamma a b-    | a < 1 -    = do+-- derived from  Marsaglia & Tang, "A Simple Method for generating gamma+-- variables", ACM Transactions on Mathematical Software, Vol 26, No 3 (2000), p363-372.+{-# SPECIALIZE mtGamma :: Double -> Double -> RVar Double #-}+{-# SPECIALIZE mtGamma :: Float -> Float -> RVar Float #-}+mtGamma+    :: (Floating a, Ord a,+        Distribution StdUniform a, +        Distribution Normal a)+    => a -> a -> RVar a+mtGamma a b +    | a < 1     = do         u <- stdUniform-        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 <- stdNormal-                -                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 <- stdUniform-                        -                        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, Floating b, Ord b, Distribution Normal b, Distribution StdUniform b) => a -> RVar b-realFloatErlang a-    | a < 1 -    = fail "realFloatErlang: a < 1"-    | otherwise-    = go-        where-            d = fromIntegral a - (1 / 3)-            c = recip (3 * sqrt d) -- (1 / 3) / sqrt d+        mtGamma (1+a) $! (b * u ** recip a)+    | otherwise = go+    where+        !d = a - fromRational (1%3)+        !c = recip (sqrt (9*d))+        +        go = do+            x <- stdNormal+            let !v   = 1 + c*x             -            go = do-                x <- stdNormal-                -                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 <- stdUniform-                        -                        if         u < 1 - 0.0331 * x_4-                            || log u < 0.5 * x_2  + d * (1 - v + log v)-                            then return (d * v)-                            else go+            if v <= 0+                then go+                else do+                    u  <- stdUniform+                    let !x_2 = x*x; !x_4 = x_2*x_2+                        v3 = v*v*v+                        dv = d * v3+                    if      u < 1 - 0.0331*x_4+                     || log u < 0.5 * x_2 + d - dv + d*log v3+                        then return (b*dv)+                        else go +{-# SPECIALIZE gamma :: Float  -> Float  -> RVar Float  #-}+{-# SPECIALIZE gamma :: Double -> Double -> RVar Double #-} gamma :: (Distribution Gamma a) => a -> a -> RVar a gamma a b = rvar (Gamma a b) @@ -99,7 +59,9 @@ data Erlang a b = Erlang a  instance (Floating a, Ord a, Distribution Normal a, Distribution StdUniform a) => Distribution Gamma a where-    rvar (Gamma a b) = realFloatGamma a b+    {-# SPECIALIZE instance Distribution Gamma Double #-}+    {-# SPECIALIZE instance Distribution Gamma Float #-}+    rvar (Gamma a b) = mtGamma a b  instance (Integral a, Floating b, Ord b, Distribution Normal b, Distribution StdUniform b) => Distribution (Erlang a) b where-    rvar (Erlang a) = realFloatErlang a+    rvar (Erlang a) = mtGamma (fromIntegral a) 1

@@ -0,0 +1,31 @@+{-# LANGUAGE GADTs, MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}+module Data.Random.Distribution.Multinomial where++import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Binomial++multinomial :: Distribution (Multinomial p) [a] => [p] -> a -> RVar [a]+multinomial ps n = rvar (Multinomial ps n)++data Multinomial p a where+    Multinomial :: [p] -> a -> Multinomial p [a]++instance (Num a, Fractional p, Distribution (Binomial p) a) => Distribution (Multinomial p) [a] where+    -- TODO: implement faster version based on Categorical for small n, large (length ps)+    rvar (Multinomial ps0 t) = go t ps0 (tailSums ps0) id+        where+            go _ []     _            f = return (f [])+            go n [_]    _            f = return (f [n])+            go 0 (_:ps) (_   :psums) f = go 0 ps psums (f . (0:))+            go n (p:ps) (psum:psums) f = do+                x <- binomial n (p / psum)+                go (n-x) ps psums (f . (x:))+            +            go _ _ _ _ = error "rvar/Multinomial: programming error! this case should be impossible!"+            +            -- less wasteful version of (map sum . tails)+            tailSums [] = [0]+            tailSums (x:xs) = case tailSums xs of+                (s:rest) -> (x+s):s:rest+                _ -> error "rvar/Multinomial/tailSums: programming error! this case should be impossible!"

@@ -1,9 +1,6 @@-{-- -      ``Data/Random/Distribution/Normal''- -} {-# LANGUAGE     MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,-    UndecidableInstances, ForeignFunctionInterface+    UndecidableInstances, ForeignFunctionInterface, BangPatterns   #-}  module Data.Random.Distribution.Normal@@ -31,8 +28,9 @@ import Data.Random.Distribution.Ziggurat import Data.Random.RVar -import Control.Monad-import Foreign.Storable+import Data.Vector.Generic (Vector)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as UV  import Data.Number.Erf @@ -82,12 +80,12 @@ normalTail :: (Distribution StdUniform a, Floating a, Ord a) =>               a -> RVar a normalTail r = go-    where +    where         go = do-            u <- stdUniform-            v <- stdUniform-            let x = log u / r-                y = log v+            !u <- stdUniform+            let !x = log u / r+            !v <- stdUniform+            let !y = log v             if x*x + y+y > 0                 then go                 else return (r - x)@@ -95,8 +93,8 @@ -- |Construct a 'Ziggurat' for sampling a normal distribution, given -- @logBase 2 c@ and the 'zGetIU' implementation. normalZ ::-  (RealFloat a, Erf a, Storable a, Distribution Uniform a, Integral b) =>-  b -> RVar (Int, a) -> Ziggurat a+  (RealFloat a, Erf a, Vector v a, Distribution Uniform a, Integral b) =>+  b -> RVar (Int, a) -> Ziggurat v a normalZ p = mkZigguratRec True normalF normalFInv normalFInt normalFVol (2^p)  -- | Ziggurat target function (upper half of a non-normalized gaussian PDF)@@ -119,7 +117,7 @@ -- |A random variable sampling from the standard normal distribution -- over any 'RealFloat' type (subject to the rest of the constraints - -- it builds and uses a 'Ziggurat' internally, which requires the 'Erf'--- and 'Storable' classes).  +-- class).   --  -- Because it computes a 'Ziggurat', it is very expensive to use for -- just one evaluation, or even for multiple evaluations if not used and@@ -132,18 +130,20 @@ -- -- As far as I know, this should be safe to use in any monomorphic -- @Distribution Normal@ instance declaration.-realFloatStdNormal :: (RealFloat a, Erf a, Storable a, Distribution Uniform a) => RVar a-realFloatStdNormal = runZiggurat (normalZ p getIU)+realFloatStdNormal :: (RealFloat a, Erf a, Distribution Uniform a) => RVar a+realFloatStdNormal = runZiggurat (normalZ p getIU `asTypeOf` (undefined :: Ziggurat V.Vector a))     where +        p :: Int         p = 6                  getIU = do-            i <- getRandomByte+            i <- getRandomPrim PrimWord8             u <- uniform (-1) 1             return (fromIntegral i .&. (2^p-1), u)  -- |A random variable sampling from the standard normal distribution -- over the 'Double' type.+{-# NOINLINE doubleStdNormal #-} doubleStdNormal :: RVar Double doubleStdNormal = runZiggurat doubleStdNormalZ @@ -154,7 +154,8 @@ doubleStdNormalR = 3.852046150368388 doubleStdNormalV = 2.4567663515413507e-3 -doubleStdNormalZ :: Ziggurat Double+{-# NOINLINE doubleStdNormalZ #-}+doubleStdNormalZ :: Ziggurat UV.Vector Double doubleStdNormalZ = mkZiggurat_ True          normalF normalFInv          doubleStdNormalC doubleStdNormalR doubleStdNormalV @@ -162,23 +163,25 @@         (normalTail doubleStdNormalR)     where          getIU = do-            w <- getRandomWord+            !w <- getRandomPrim PrimWord64             let (u,i) = wordToDoubleWithExcess w-            return (fromIntegral i .&. (doubleStdNormalC-1), u+u-1)+            return $! (fromIntegral i .&. (doubleStdNormalC-1), u+u-1)  -- |A random variable sampling from the standard normal distribution -- over the 'Float' type.+{-# NOINLINE floatStdNormal #-} floatStdNormal :: RVar Float floatStdNormal = runZiggurat floatStdNormalZ --- floatStdNormalC must not be over 2^41 if using wordToFloatWithExcess+-- floatStdNormalC must not be over 2^9 if using word32ToFloatWithExcess floatStdNormalC :: Int floatStdNormalC = 512 floatStdNormalR, floatStdNormalV :: Float floatStdNormalR = 3.852046150368388 floatStdNormalV = 2.4567663515413507e-3 -floatStdNormalZ :: Ziggurat Float+{-# NOINLINE floatStdNormalZ #-}+floatStdNormalZ :: Ziggurat UV.Vector Float floatStdNormalZ = mkZiggurat_ True          normalF normalFInv          floatStdNormalC floatStdNormalR floatStdNormalV @@ -186,13 +189,10 @@         (normalTail floatStdNormalR)     where         getIU = do-            w <- getRandomWord-            let (u,i) = wordToFloatWithExcess w+            !w <- getRandomPrim PrimWord32+            let (u,i) = word32ToFloatWithExcess w             return (fromIntegral i .&. (floatStdNormalC-1), u+u-1) -normalPdf :: Real a => a -> a -> a -> Double-normalPdf m s x = recip (realToFrac s * sqrt (2*pi)) * exp (-0.5 * (realToFrac x - realToFrac m)^2 / (realToFrac s)^2)- normalCdf :: (Real a) => a -> a -> a -> Double normalCdf m s x = normcdf ((realToFrac x - realToFrac m) / realToFrac s) @@ -200,7 +200,7 @@ data Normal a     -- |The \"standard\" normal distribution - mean 0, stddev 1     = StdNormal-    -- |@Normal m s@ is a normal distribution with mean @m@ and stddev @s@.+    -- |@Normal m s@ is a normal distribution with mean @m@ and stddev @sd@.     | Normal a a -- mean, sd  instance Distribution Normal Double where

@@ -17,16 +17,14 @@ 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, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => b -> RVar a-integralPoisson mu = psn 0 mu+integralPoisson = psn 0     where         psn :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => a -> b -> RVar a-        psn k mu+        psn j mu             | mu > 10   = do                 let m = floor (mu * (7/8))             @@ -34,10 +32,10 @@                 if x >= mu                     then do                         b <- binomial (m - 1) (mu / x)-                        return (k + b)-                    else psn (k + m) (mu - x)+                        return (j + b)+                    else psn (j + m) (mu - x)             -            | otherwise = prod 1 k+            | otherwise = prod 1 j                 where                     emu = exp (-mu)                 

@@ -28,12 +28,12 @@     triUpper  :: a}     deriving (Eq, Show) --- |Compute a triangular distribution for a 'Fractional' type.  The name is--- a historical accident and may change in the future.-realFloatTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVar a-realFloatTriangular a b c-    | a <= b && b <= c-    = do+-- |Compute a triangular distribution for a 'Floating' type.+floatingTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVar a+floatingTriangular a b c+    | a > b     = floatingTriangular b a c+    | b > c     = floatingTriangular a c b+    | otherwise = do         let p = (c-b)/(c-a)         u <- stdUniform         let d   | u >= p    = a@@ -44,19 +44,19 @@ --        x <- stdUniform         return (b - ((1 - sqrt x) * (b-d))) --- |@realFloatTriangularCDF a b c@ is the CDF of @realFloatTriangular a b c@.-realFloatTriangularCDF :: RealFrac a => a -> a -> a -> a -> Double-realFloatTriangularCDF a b c x+-- |@triangularCDF a b c@ is the CDF of @realFloatTriangular a b c@.+triangularCDF :: RealFrac a => a -> a -> a -> a -> Double+triangularCDF a b c x     | x < a     = 0     | x <= b-    = realToFrac ((x - a) ^ 2 / ((c - a) * (b - a)))+    = realToFrac ((x - a)^(2 :: Int) / ((c - a) * (b - a)))     | x <= c-    = realToFrac (1 - (c - x) ^ 2 / ((c - a) * (c - b)))+    = realToFrac (1 - (c - x)^(2 :: Int) / ((c - a) * (c - b)))     | otherwise     = 1      instance (RealFloat a, Ord a, Distribution StdUniform a) => Distribution Triangular a where-    rvar (Triangular a b c) = realFloatTriangular a b c+    rvar (Triangular a b c) = floatingTriangular a b c instance (RealFrac a, Distribution Triangular a) => CDF Triangular a where-    cdf  (Triangular a b c) = realFloatTriangularCDF a b c+    cdf  (Triangular a b c) = triangularCDF a b c

@@ -5,7 +5,8 @@     MultiParamTypeClasses, FunctionalDependencies,     FlexibleContexts, FlexibleInstances,      UndecidableInstances, EmptyDataDecls,-    TemplateHaskell+    TemplateHaskell,+    BangPatterns   #-}  module Data.Random.Distribution.Uniform@@ -14,7 +15,6 @@ 	     , StdUniform(..)     , stdUniform-    , stdUniformNonneg     , stdUniformPos          , integralUniform@@ -50,30 +50,55 @@ import Control.Monad.Loops  -- |Compute a random 'Integral' value between the 2 values provided (inclusive).+{-# INLINE integralUniform #-} integralUniform :: (Integral a) => a -> a -> RVar a-integralUniform a b-    | a > b     = compute b a-    | otherwise = compute a b+integralUniform !x !y = if x < y then integralUniform' x y else integralUniform' y x++{-# SPECIALIZE integralUniform' :: Int   -> Int   -> RVar Int   #-}+{-# SPECIALIZE integralUniform' :: Int8  -> Int8  -> RVar Int8  #-}+{-# SPECIALIZE integralUniform' :: Int16 -> Int16 -> RVar Int16 #-}+{-# SPECIALIZE integralUniform' :: Int32 -> Int32 -> RVar Int32 #-}+{-# SPECIALIZE integralUniform' :: Int64 -> Int64 -> RVar Int64 #-}+{-# SPECIALIZE integralUniform' :: Word   -> Word   -> RVar Word   #-}+{-# SPECIALIZE integralUniform' :: Word8  -> Word8  -> RVar Word8  #-}+{-# SPECIALIZE integralUniform' :: Word16 -> Word16 -> RVar Word16 #-}+{-# SPECIALIZE integralUniform' :: Word32 -> Word32 -> RVar Word32 #-}+{-# SPECIALIZE integralUniform' :: Word64 -> Word64 -> RVar Word64 #-}+{-# SPECIALIZE integralUniform' :: Integer -> Integer -> RVar Integer #-}+integralUniform' :: (Integral a) => a -> a -> RVar a+integralUniform' !l !u+    | nReject == 0  = fmap shift prim+    | otherwise     = fmap shift loop     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))+        m = 1 + toInteger u - toInteger l+        (bytes, nPossible) = bytesNeeded m+        nReject = nPossible `mod` m+        +        !prim = getRandomPrim (PrimNByteInteger bytes)+        !shift = \(!z) -> l + (fromInteger $! (z `mod` m))+        +        loop = do+            z <- prim+            if z < nReject+                then loop+                else return z +integralUniformCDF :: (Integral a, Fractional b) => a -> a -> a -> b integralUniformCDF a b x     | b < a     = integralUniformCDF b a x     | x < a     = 0     | x > b     = 1     | otherwise = (fromIntegral x - fromIntegral a) / (fromIntegral b - fromIntegral a) -bytesNeeded x = case findIndex (> x) powersOf256 of-    Just x -> x-powersOf256 = iterate (256 *) 1+-- TODO: come up with a decent, fast heuristic to decide whether to return an extra+-- byte.  May involve moving calculation of nReject into this function, and then+-- accepting first if 4*nReject < nPossible or something similar.+bytesNeeded :: Integer -> (Int, Integer)+bytesNeeded x = head (dropWhile ((<= x).snd) powersOf256) +powersOf256 :: [(Int, Integer)]+powersOf256 = zip [0..] (iterate (256 *) 1)+ -- |Compute a random value for a 'Bounded' type, between 'minBound' and 'maxBound' -- (inclusive for 'Integral' or 'Enum' types, in ['minBound', 'maxBound') for Fractional types.) boundedStdUniform :: (Distribution Uniform a, Bounded a) => RVar a@@ -93,12 +118,13 @@ -- |Compute a uniform random 'Float' value in the range [0,1) floatStdUniform :: RVar Float floatStdUniform = do-    x <- getRandomWord-    return (wordToFloat x)+    x <- getRandomPrim PrimWord32+    return (word32ToFloat x)  -- |Compute a uniform random 'Double' value in the range [0,1)+{-# INLINE doubleStdUniform #-} doubleStdUniform :: RVar Double-doubleStdUniform = getRandomDouble+doubleStdUniform = getRandomPrim PrimDouble  -- |Compute a uniform random value in the range [0,1) for any 'RealFloat' type  realFloatStdUniform :: RealFloat a => RVar a@@ -137,6 +163,7 @@     return (a + x * (b - a))  -- |@doubleUniform a b@ computes a uniform random 'Double' value in the range [a,b)+{-# INLINE doubleUniform #-} doubleUniform :: Double -> Double -> RVar Double doubleUniform 0 1 = doubleStdUniform doubleUniform a b = do@@ -190,22 +217,22 @@  -- |Get a \"standard\" uniformly distributed value. -- For integral types, this means uniformly distributed over the full range--- of the type (and hence there is no support for Integer).  For fractional+-- of the type (there is no support for 'Integer').  For fractional -- types, this means uniformly distributed on the interval [0,1).+{-# SPECIALIZE stdUniform :: RVar Double #-}+{-# SPECIALIZE stdUniform :: RVar Float #-} stdUniform :: (Distribution StdUniform a) => RVar a stdUniform = rvar StdUniform --- |Like 'stdUniform', but uses 'abs' to return only positive or zero values.+-- |Like 'stdUniform', but returns only positive or zero values.  Not +-- exported because it is not truly uniform: nonzero values are twice+-- as likely as zero on signed types. stdUniformNonneg :: (Distribution StdUniform a, Num a) => RVar a-stdUniformNonneg = abs `fmap` stdUniform+stdUniformNonneg = fmap abs stdUniform  -- |Like 'stdUniform' but only returns positive values. stdUniformPos :: (Distribution StdUniform a, Num a) => RVar a-stdUniformPos = do-    x <- stdUniformNonneg-    if x /= 0-        then return x-        else stdUniformPos+stdUniformPos = iterateUntil (/= 0) stdUniformNonneg  -- |A definition of a uniform distribution over the type @t@.  See also 'uniform'. data Uniform t = @@ -229,15 +256,34 @@         instance CDF Uniform Int            where cdf  (Uniform a b) = integralUniformCDF a b     |]) --- Some integral types have specialized StdUniform rvars:-instance Distribution StdUniform Int8       where rvar ~StdUniform = fmap fromIntegral getRandomByte-instance Distribution StdUniform Word8      where rvar ~StdUniform = getRandomByte-instance Distribution StdUniform Word64     where rvar ~StdUniform = getRandomWord--- and Integer has none...-$( replicateInstances ''Int (integralTypes \\ [''Int8, ''Word8, ''Word64, ''Integer]) [d|-        instance Distribution StdUniform Int    where rvar ~StdUniform = fmap fromIntegral getRandomWord-    |])+instance Distribution StdUniform Word8      where rvarT ~StdUniform = getRandomPrim PrimWord8+instance Distribution StdUniform Word16     where rvarT ~StdUniform = getRandomPrim PrimWord16+instance Distribution StdUniform Word32     where rvarT ~StdUniform = getRandomPrim PrimWord32+instance Distribution StdUniform Word64     where rvarT ~StdUniform = getRandomPrim PrimWord64 +instance Distribution StdUniform Int8       where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord8+instance Distribution StdUniform Int16      where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord16+instance Distribution StdUniform Int32      where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord32+instance Distribution StdUniform Int64      where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord64++instance Distribution StdUniform Int where+    rvar+        | toInteger (maxBound :: Int) > toInteger (maxBound :: Int32)+        = const (fromIntegral `fmap` getRandomPrim PrimWord64)+        +        | otherwise+        = const (fromIntegral `fmap` getRandomPrim PrimWord32)++instance Distribution StdUniform Word where+    rvar+        | toInteger (maxBound :: Word) > toInteger (maxBound :: Word32)+        = const (fromIntegral `fmap` getRandomPrim PrimWord64)+        +        | otherwise+        = const (fromIntegral `fmap` getRandomPrim PrimWord32)++-- Integer has no StdUniform...+ $( replicateInstances ''Int (integralTypes \\ [''Integer]) [d|         instance CDF StdUniform Int         where cdf  ~StdUniform = boundedStdUniformCDF     |])@@ -249,7 +295,7 @@ instance CDF Uniform Double                 where cdf  (Uniform a b) = realUniformCDF a b  instance Distribution StdUniform Float      where rvar ~StdUniform = floatStdUniform-instance Distribution StdUniform Double     where rvar ~StdUniform = doubleStdUniform+instance Distribution StdUniform Double     where rvar ~StdUniform = getRandomPrim PrimDouble; rvarT ~StdUniform = getRandomPrim PrimDouble instance CDF StdUniform Float               where cdf  ~StdUniform = realStdUniformCDF instance CDF StdUniform Double              where cdf  ~StdUniform = realStdUniformCDF @@ -262,18 +308,18 @@ instance HasResolution r =>           CDF StdUniform (Fixed r)           where cdf  ~StdUniform = realStdUniformCDF -instance Distribution Uniform ()            where rvar (Uniform a b) = return ()-instance CDF Uniform ()                     where cdf  (Uniform a b) = return 1+instance Distribution Uniform ()            where rvar (Uniform _ _) = return ()+instance CDF Uniform ()                     where cdf  (Uniform _ _) = return 1 $( replicateInstances ''Char [''Char, ''Bool, ''Ordering] [d|         instance Distribution Uniform Char  where rvar (Uniform a b) = enumUniform a b         instance CDF Uniform Char           where cdf  (Uniform a b) = enumUniformCDF a b      |]) -instance Distribution StdUniform ()         where rvar ~StdUniform = return ()-instance CDF StdUniform ()                  where cdf  ~StdUniform = return 1-instance Distribution StdUniform Bool       where rvar ~StdUniform = fmap even getRandomByte-instance CDF StdUniform Bool                where cdf  ~StdUniform = boundedEnumStdUniformCDF+instance Distribution StdUniform ()         where rvarT ~StdUniform = return ()+instance CDF StdUniform ()                  where cdf   ~StdUniform = return 1+instance Distribution StdUniform Bool       where rvarT ~StdUniform = fmap even (getRandomPrim PrimWord8)+instance CDF StdUniform Bool                where cdf   ~StdUniform = boundedEnumStdUniformCDF  instance Distribution StdUniform Char       where rvar ~StdUniform = boundedEnumStdUniform instance CDF StdUniform Char                where cdf  ~StdUniform = boundedEnumStdUniformCDF

@@ -0,0 +1,16 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}+module Data.Random.Distribution.Weibull where++import Data.Random.Distribution+import Data.Random.Distribution.Uniform++data Weibull a = Weibull { weibullLambda :: !a, weibullK :: !a }+    deriving (Eq, Show)++instance (Floating a, Distribution StdUniform a) => Distribution Weibull a where+    rvarT (Weibull lambda k) = do+        u <- rvarT StdUniform+        return (lambda * (negate (log u)) ** recip k)++instance (Real a, Distribution Weibull a) => CDF Weibull a where+    cdf (Weibull lambda k) x = 1 - exp (negate ((realToFrac x / realToFrac lambda) ** realToFrac k))

@@ -1,7 +1,7 @@ {-# LANGUAGE         MultiParamTypeClasses,         FlexibleInstances, FlexibleContexts,-        RecordWildCards+        RecordWildCards, BangPatterns   #-}  -- |A generic \"ziggurat algorithm\" implementation.  Fairly rough right@@ -31,10 +31,10 @@ import Data.Random.Distribution.Uniform import Data.Random.Distribution import Data.Random.RVar-import Data.StorableVector as Vec-import Foreign.Storable--vec ! i = index vec i+import Data.Vector.Generic as Vec+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as UV+import Data.Function (fix)  -- |A data structure containing all the data that is needed -- to implement Marsaglia & Tang's \"ziggurat\" algorithm for@@ -45,7 +45,7 @@ -- be.  There are several helper functions that will build 'Ziggurat's. -- The pathologically curious may wish to read the 'runZiggurat' source. -- That is the ultimate specification of the semantics of all these fields.-data Ziggurat t = Ziggurat {+data Ziggurat v t = Ziggurat {         -- |The X locations of each bin in the distribution.  Bin 0 is the         -- 'infinite' one.         -- @@ -55,11 +55,11 @@         -- faster by not needing to specially-handle bin 0 quite as often.         -- If you really need to know why it works, see the 'runZiggurat'         -- source or \"the literature\" - it's a fairly standard setup.-        zTable_xs         :: Vector t,+        zTable_xs         :: !(v t),         -- |The ratio of each bin's Y value to the next bin's Y value-        zTable_x_ratios   :: Vector t,+        zTable_y_ratios   :: !(v t),         -- |The Y value (zFunc x) of each bin-        zTable_ys         :: Vector t,+        zTable_ys         :: !(v t),         -- |An RVar providing a random tuple consisting of:         --         --  * a bin index, uniform over [0,c) :: Int (where @c@ is the@@ -73,58 +73,63 @@         -- single random word (64 bits) can be efficiently converted to         -- a double (using 52 bits) and a bin number (using up to 12 bits),         -- for example.-        zGetIU            :: RVar (Int, t),+        zGetIU            :: !(RVar (Int, t)),                  -- |The distribution for the final \"virtual\" bin         -- (the ziggurat algorithm does not handle distributions         -- that wander off to infinity, so another distribution is needed         -- to handle the last \"bin\" that stretches to infinity)-        zTailDist         :: RVar t,+        zTailDist         :: (RVar t),                  -- |A copy of the uniform RVar generator for the base type,         -- so that @Distribution Uniform t@ is not needed when sampling         -- from a Ziggurat (makes it a bit more self-contained).-        zUniform          :: t -> t -> RVar t,+        zUniform          :: !(t -> t -> RVar t),                  -- |The (one-sided antitone) PDF, not necessarily normalized-        zFunc             :: t -> t,+        zFunc             :: !(t -> t),                  -- |A flag indicating whether the distribution should be         -- mirrored about the origin (the ziggurat algorithm it         -- its native form only samples from one-sided distributions.         -- By mirroring, we can extend it to symmetric distributions         -- such as the normal distribution)-        zMirror           :: Bool+        zMirror           :: !Bool     }  -- |Sample from the distribution encoded in a 'Ziggurat' data structure. {-# INLINE runZiggurat #-}-runZiggurat :: (Num a, Ord a, Storable a) =>-               Ziggurat a -> RVar a-runZiggurat Ziggurat{..} = go+{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Float  -> RVar Float #-}+{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Double -> RVar Double #-}+{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Float  -> RVar Float #-}+{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Double -> RVar Double #-}+runZiggurat :: (Num a, Ord a, Vector v a) =>+               Ziggurat v a -> RVar a+runZiggurat !Ziggurat{..} = go     where+        {-# NOINLINE go #-}         go = do             -- Select a bin (I) and a uniform value (U) from -1 to 1             -- (or 0 to 1 if not mirroring the distribution).             -- Let X be U scaled to the size of the selected bin.-            (i,u) <- zGetIU-            let x = u * zTable_xs ! i+            (!i,!u) <- zGetIU                          -- if the uniform value U falls in the area "clearly inside" the             -- bin, accept X immediately.             -- Otherwise, depending on the bin selected, use either the             -- tail distribution or an accept/reject test.-            if abs u < zTable_x_ratios ! i-                then return $! x+            if abs u < zTable_y_ratios ! i+                then return $! (u * zTable_xs ! i)                 else if i == 0-                    then sampleTail x-                    else sampleGreyArea i x+                    then sampleTail u+                    else sampleGreyArea i $! (u * zTable_xs ! i)                  -- when the sample falls in the "grey area" (the area between         -- the Y values of the selected bin and the bin after that one),         -- use an accept/reject method based on the target PDF.+        {-# INLINE sampleGreyArea #-}         sampleGreyArea i x = do-            v <- zUniform (zTable_ys ! (i+1)) (zTable_ys ! i)+            !v <- zUniform (zTable_ys ! (i+1)) (zTable_ys ! i)             if v < zFunc (abs x)                 then return $! x                 else go@@ -132,9 +137,10 @@         -- if the selected bin is the "infinite" one, call it quits and         -- defer to the tail distribution (mirroring if needed to ensure         -- the result has the sign already selected by zGetIU)+        {-# INLINE sampleTail #-}         sampleTail x-            | x < 0     = fmap negate zTailDist-            | otherwise = zTailDist+            | zMirror && x < 0  = fmap negate zTailDist+            | otherwise         = zTailDist   -- |Build the tables to implement the \"ziggurat algorithm\" devised by @@ -159,7 +165,12 @@ --  --  * an RVar sampling from the tail (the region where x > R) -- -mkZiggurat_ :: (RealFloat t, Storable t,+{-# INLINE mkZiggurat_ #-}+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> RVar (Int,  Float) -> RVar Float  -> Ziggurat UV.Vector Float #-}+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> RVar (Int, Double) -> RVar Double -> Ziggurat UV.Vector Double #-}+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> RVar (Int,  Float) -> RVar Float  -> Ziggurat V.Vector Float #-}+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> RVar (Int, Double) -> RVar Double -> Ziggurat V.Vector Double #-}+mkZiggurat_ :: (RealFloat t, Vector v t,                Distribution Uniform t) =>               Bool               -> (t -> t)@@ -169,18 +180,19 @@               -> t               -> RVar (Int, t)               -> RVar t-              -> Ziggurat t-mkZiggurat_ m f fInv c r v getIU tailDist = z-    where z = Ziggurat-            { zTable_xs         = zigguratTable f fInv c r v-            , zTable_x_ratios   = precomputeRatios (zTable_xs z)-            , zTable_ys         = Vec.map f (zTable_xs z)-            , zGetIU            = getIU-            , zUniform          = uniform-            , zFunc             = f-            , zTailDist         = tailDist-            , zMirror           = m-            }+              -> Ziggurat v t+mkZiggurat_ m f fInv c r v getIU tailDist = Ziggurat+    { zTable_xs         = xs+    , zTable_y_ratios   = precomputeRatios xs+    , zTable_ys         = Vec.map f xs+    , zGetIU            = getIU+    , zUniform          = uniform+    , zFunc             = f+    , zTailDist         = tailDist+    , zMirror           = m+    }+    where +        xs = zigguratTable f fInv c r v  -- |Build the tables to implement the \"ziggurat algorithm\" devised by  -- Marsaglia & Tang, attempting to automatically compute the R and V@@ -189,7 +201,7 @@ -- Arguments are the same as for 'mkZigguratRec', with an additional -- argument for the tail distribution as a function of the selected -- R value.-mkZiggurat :: (RealFloat t, Storable t,+mkZiggurat :: (RealFloat t, Vector v t,                Distribution Uniform t) =>               Bool               -> (t -> t)@@ -199,7 +211,7 @@               -> Int               -> RVar (Int, t)               -> (t -> RVar t)-              -> Ziggurat t+              -> Ziggurat v t mkZiggurat m f fInv fInt fVol c getIU tailDist =     mkZiggurat_ m f fInv c r v getIU (tailDist r)          where@@ -226,7 +238,7 @@ --  * an RVar providing the 'zGetIU' random tuple -- mkZigguratRec ::-  (RealFloat t, Storable t,+  (RealFloat t, Vector v t,    Distribution Uniform t) =>   Bool   -> (t -> t)@@ -235,14 +247,23 @@   -> t   -> Int   -> RVar (Int, t)-  -> Ziggurat t-mkZigguratRec m f fInv fInt fVol c getIU = -    mkZiggurat m f fInv fInt fVol c getIU (fix (mkTail m f fInv fInt fVol c getIU))+  -> Ziggurat v t+mkZigguratRec m f fInv fInt fVol c getIU = z         where-            fix f = f (fix f)+            z = mkZiggurat m f fInv fInt fVol c getIU (fix (mkTail m f fInv fInt fVol c getIU z)) -mkTail m f fInv fInt fVol c getIU nextTail r = do-     x <- rvar (mkZiggurat m f' fInv' fInt' fVol' c getIU nextTail)+mkTail :: +    (RealFloat a, Vector v a, Distribution Uniform a) =>+    Bool+    -> (a -> a) -> (a -> a) -> (a -> a)+    -> a+    -> Int+    -> RVar (Int, a)+    -> Ziggurat v a+    -> (a -> RVar a)+    -> (a -> RVar a)+mkTail m f fInv fInt fVol c getIU typeRep nextTail r = do+     x <- rvar (mkZiggurat m f' fInv' fInt' fVol' c getIU nextTail `asTypeOf` typeRep)      return (x + r * signum x)         where             fIntR = fInt r@@ -256,14 +277,15 @@             fVol' = fVol - fIntR          -zigguratTable :: (Fractional a, Storable a, Ord a) =>-                 (a -> a) -> (a -> a) -> Int -> a -> a -> Vector a+zigguratTable :: (Fractional a, Vector v a, Ord a) =>+                 (a -> a) -> (a -> a) -> Int -> a -> a -> v a zigguratTable f fInv c r v = case zigguratXs f fInv c r v of-    (xs, excess) -> pack xs-    where epsilon = 1e-3*v+    (xs, _excess) -> fromList xs +zigguratExcess :: (Fractional a, Ord a) => (a -> a) -> (a -> a) -> Int -> a -> a -> a zigguratExcess f fInv c r v = snd (zigguratXs f fInv c r v) +zigguratXs :: (Fractional a, Ord a) => (a -> a) -> (a -> a) -> Int -> a -> a -> ([a], a) zigguratXs f fInv c r v = (xs, excess)     where         xs = Prelude.map x [0..c] -- sample c x@@ -273,6 +295,7 @@         x 1 = r         x i | i == c = 0         x (i+1) = next i+        x _ = error "zigguratXs: programming error! this case should be impossible!"                  next i = let x_i = xs!!i                   in if x_i <= 0 then -1 else fInv (ys!!i + (v / x_i))@@ -280,7 +303,8 @@         excess = xs!!(c-1) * (f 0 - ys !! (c-1)) - v   -precomputeRatios zTable_xs = sample (c-1) $ \i -> zTable_xs!(i+1) / zTable_xs!i+precomputeRatios :: (Vector v a, Fractional a) => v a -> v a+precomputeRatios zTable_xs = generate (c-1) $ \i -> zTable_xs!(i+1) / zTable_xs!i     where         c = Vec.length zTable_xs @@ -302,7 +326,7 @@ -- Result: (R,V) findBin0 :: (RealFloat b) =>      Int -> (b -> b) -> (b -> b) -> (b -> b) -> b -> (b, b)-findBin0 cInt f fInv fInt fVol = (r,v r)+findBin0 cInt f fInv fInt fVol = (rMin,v rMin)     where         c = fromIntegral cInt         v r = r * f r + fVol - fInt r@@ -310,11 +334,11 @@         -- initial R guess:         r0 = findMin (\r -> v r <= fVol / c)         -- find a better R:-        r = findMinFrom r0 1 $ \r -> +        rMin = findMinFrom r0 1 $ \r ->              let e = exc r               in e >= 0 && not (isNaN e)                  exc x = zigguratExcess f fInv cInt x (v x) -instance (Num t, Ord t, Storable t) => Distribution Ziggurat t where+instance (Num t, Ord t, Vector v t) => Distribution (Ziggurat v) t where     rvar = runZiggurat

@@ -9,21 +9,25 @@ findMax :: (Fractional a, Ord a) => (a -> Bool) -> a findMax p = negate (findMin (p.negate)) -findMaxFrom :: (Fractional a, Ord a) => a -> a -> (a -> Bool) -> a-findMaxFrom z 0 p = findMaxFrom z 1 p-findMaxFrom z step1 p = findMinFrom z (negate step1) p- -- |Given an upward-closed predicate on an ordered Fractional type, -- find the smallest value satisfying the predicate. findMin :: (Fractional a, Ord a) => (a -> Bool) -> a findMin = findMinFrom 0 1 +-- |Given an upward-closed predicate on an ordered Fractional type,+-- find the smallest value satisfying the predicate.  Starts at the+-- specified point with the specified stepsize, performs an exponential+-- search out from there until it finds an interval bracketing the+-- change-point of the predicate, and then performs a bisection search+-- to isolate the change point.  Note that infinitely-divisible domains +-- such as 'Rational' cannot be searched by this function because it does+-- not terminate until it reaches a point where further subdivision of the+-- interval has no effect. findMinFrom :: (Fractional a, Ord a) => a -> a -> (a -> Bool) -> a-findMinFrom z 0 p = findMinFrom z 1 p-findMinFrom z step1 p-    | p z   = descend (z-step1) z-    | otherwise-    = fixZero (ascend z (z+step1))+findMinFrom z0 0 p = findMinFrom z0 1 p+findMinFrom z0 step1 p+    | p z0      = descend (z0-step1) z0+    | otherwise = fixZero (ascend z0 (z0+step1))     where         -- eliminate negative zero, which, in many domains, is technically         -- a feasible answer@@ -35,13 +39,13 @@         -- 0 <= l < x         ascend l x              | p x       = bisect l x-            | otherwise = ascend x $! 2*x-z+            | otherwise = ascend x $! 2*x-z0                  -- preconditions:         -- p h         -- x < h <= 0         descend x h -            | p x       = (descend $! 2*x-z) x+            | p x       = (descend $! 2*x-z0) x             | otherwise = bisect x h                  -- preconditions:@@ -49,11 +53,11 @@         -- p h         -- l <= h         bisect l h -            | l >= h    = h-            | l >= mid || mid >= h+            | l /< h    = h+            | l /< mid || mid /< h             = if p mid then mid else h             | p mid     = bisect l mid             | otherwise = bisect mid h             where -                a >= b = not (a < b)+                a /< b = not (a < b)                 mid = (l+h)*0.5

@@ -0,0 +1,235 @@+{-# LANGUAGE GADTs, RankNTypes, DeriveDataTypeable #-}+-- |This is an experimental interface to support an extensible set of primitives,+-- where a RandomSource will be able to support whatever subset of them they want+-- and have well-founded defaults generated automatically for any unsupported+-- primitives.+--+-- The purpose, in case it's not clear, is to decouple the implementations of+-- entropy sources from any particular set of primitives, so that implementors+-- of random variates can make use of a large number of primitives, supported+-- on all entropy sources, while the burden on entropy-source implementors+-- is only to provide one or two basic primitives of their choice.+-- +-- One challenge I foresee with this interface is optimization - different +-- compilers or even different versions of GHC may treat this interface +-- radically differently, making it very hard to achieve reliable performance+-- on all platforms.  It may even be that no compiler optimizes sufficiently+-- to make the flexibility this system provides worth the overhead.  I hope+-- this is not the case, but if it turns out to be a major problem, this+-- system may disappear or be modified in significant ways.+module Data.Random.Internal.Primitives (Prim(..), decomposePrimWhere) where++import Data.Random.Internal.Words+import Data.Word+import Data.Bits+import Data.Typeable++import Control.Monad.Prompt++-- |A 'Prompt' GADT describing a request for a primitive random variate.+-- Random variable definitions will request their entropy via these prompts,+-- and entropy sources will satisfy some or all of them.  The 'decomposePrimWhere'+-- function extends an entropy source's incomplete definition to a complete +-- definition, essentially defining a very flexible implementation-defaulting+-- system.+-- +-- Some possible future additions:+--    PrimFloat :: Prim Float+--    PrimInt :: Prim Int+--    PrimPair :: Prim a -> Prim b -> Prim (a :*: b)+--    PrimNormal :: Prim Double+--    PrimChoice :: [(Double :*: a)] -> Prim a+--+-- Unfortunately, I cannot get Haddock to accept my comments about the +-- data constructors, but hopefully they should be reasonably self-explanatory.+data Prim a where+    -- An unsigned byte, uniformly distributed from 0 to 0xff+    PrimWord8           :: Prim Word8+    -- An unsigned 16-bit word, uniformly distributed from 0 to 0xffff+    PrimWord16          :: Prim Word16+    -- An unsigned 32-bit word, uniformly distributed from 0 to 0xffffffff+    PrimWord32          :: Prim Word32+    -- An unsigned 64-bit word, uniformly distributed from 0 to 0xffffffffffffffff+    PrimWord64          :: Prim Word64+    -- A double-precision float U, uniformly distributed 0 <= U < 1+    PrimDouble          :: Prim Double+    -- A uniformly distributed 'Integer' 0 <= U < 2^(8*n)+    PrimNByteInteger    :: !Int -> Prim Integer+    deriving (Typeable)++instance Show (Prim a) where+    showsPrec _p PrimWord8               = showString "PrimWord8"+    showsPrec _p PrimWord16              = showString "PrimWord16"+    showsPrec _p PrimWord32              = showString "PrimWord32"+    showsPrec _p PrimWord64              = showString "PrimWord64"+    showsPrec _p PrimDouble              = showString "PrimDouble"+    showsPrec  p (PrimNByteInteger n)    = showParen (p > 10) (showString "PrimNByteInteger " . showsPrec 11 n)++-- |This is essentially a suite of interrelated default implementations,+-- each definition making use of only \"supported\" primitives.  It _really_+-- ought to be inlined to the point where the @supported@ predicate+-- is able to be inlined into it and eliminated.  +-- +-- When inlined sufficiently, it should in theory be optimized down to the+-- static set of "best" definitions for each required primitive in terms of +-- only supported primitives.+-- +-- Hopefully, when not inlined, it does not impose too much overhead.+{-# INLINE decomposePrimWhere #-}+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word8   -> Prompt Prim Word8   #-}+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word16  -> Prompt Prim Word16  #-}+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word32  -> Prompt Prim Word32  #-}+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word64  -> Prompt Prim Word64  #-}+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Double  -> Prompt Prim Double  #-}+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Integer -> Prompt Prim Integer #-}+decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim a -> Prompt Prim a+decomposePrimWhere supported requested = decomp requested+    where+        {-# INLINE decomp #-}++        {-# SPECIALIZE decomp :: Prim Word8   -> Prompt Prim Word8   #-}+        {-# SPECIALIZE decomp :: Prim Word16  -> Prompt Prim Word16  #-}+        {-# SPECIALIZE decomp :: Prim Word32  -> Prompt Prim Word32  #-}+        {-# SPECIALIZE decomp :: Prim Word64  -> Prompt Prim Word64  #-}+        {-# SPECIALIZE decomp :: Prim Double  -> Prompt Prim Double  #-}+        {-# SPECIALIZE decomp :: Prim Integer -> Prompt Prim Integer #-}+        -- First, all supported prims should just be evaluated directly.+        decomp :: Prim a -> Prompt Prim a+        decomp prim+            | supported prim = prompt prim+        -- beyond this point, all definitions must be in terms of+        -- 'prompt's referring to other supported primitives or +        -- 'decomp's referring to other primitives in a well-founded way+        +        decomp PrimWord8+            | supported PrimWord16 = do+                w <- prompt PrimWord16+                return (fromIntegral w)+            | supported PrimWord32 = do+                w <- prompt PrimWord32+                return (fromIntegral w)+            | supported PrimWord64 = do+                w <- prompt PrimWord64+                return (fromIntegral w)+            | supported PrimDouble = do+                d <- prompt PrimDouble+                return (truncate (d * 256))+            | supported (PrimNByteInteger 1) = do+                i <- prompt (PrimNByteInteger 1)+                return (fromInteger i)+        +        decomp PrimWord16+            | supported PrimWord8 = do+                b0 <- prompt PrimWord8+                b1 <- prompt PrimWord8+                return (buildWord16 b0 b1)+            | supported PrimWord32 = do+                w <- prompt PrimWord32+                return (fromIntegral w)+            | supported PrimWord64 = do+                w <- prompt PrimWord64+                return (fromIntegral w)+            | supported PrimDouble = do+                d <- prompt PrimDouble+                return (truncate (d * 65536))+            | supported (PrimNByteInteger 2) = do+                i <- prompt (PrimNByteInteger 2)+                return (fromInteger i)+        +        decomp PrimWord32+            | supported PrimWord16 = do+                w0 <- prompt PrimWord16+                w1 <- prompt PrimWord16+                +                return (buildWord32' w0 w1)+            | supported PrimWord8 = do+                b0 <- prompt PrimWord8+                b1 <- prompt PrimWord8+                b2 <- prompt PrimWord8+                b3 <- prompt PrimWord8+                +                return (buildWord32 b0 b1 b2 b3)+            | supported PrimWord64 = do+                w <- prompt PrimWord64+                return (fromIntegral w)+            | supported PrimDouble = do+                d <- prompt PrimDouble+                return (truncate (d * 4294967296))+            | supported (PrimNByteInteger 4) = do+                i <- prompt (PrimNByteInteger 4)+                return (fromInteger i)+        +        decomp PrimWord64+            | supported PrimWord32 = do+                w0 <- prompt PrimWord32+                w1 <- prompt PrimWord32+                +                return (buildWord64'' w0 w1)+            | supported PrimWord16 = do+                w0 <- prompt PrimWord16+                w1 <- prompt PrimWord16+                w2 <- prompt PrimWord16+                w3 <- prompt PrimWord16+                +                return (buildWord64' w0 w1 w2 w3)+            | supported PrimWord8 = do+                b0 <- prompt PrimWord8+                b1 <- prompt PrimWord8+                b2 <- prompt PrimWord8+                b3 <- prompt PrimWord8+                b4 <- prompt PrimWord8+                b5 <- prompt PrimWord8+                b6 <- prompt PrimWord8+                b7 <- prompt PrimWord8+                +                return (buildWord64 b0 b1 b2 b3 b4 b5 b6 b7)+            | supported PrimDouble = do+                -- Need 2 doubles, because a uniform [0,1) double only has+                -- about 52 bits of reliable entropy+                d0 <- prompt PrimDouble+                d1 <- prompt PrimDouble+                +                let w0 = truncate (d0 * 4294967296)+                    w1 = truncate (d1 * 4294967296)+                +                return (w0 .|. (w1 `shiftL` 32))+            | supported (PrimNByteInteger 8) = do+                i <- prompt (PrimNByteInteger 8)+                return (fromInteger i)+        +        decomp PrimDouble = do+            word <- decomp PrimWord64+            return (wordToDouble word)+        +        decomp (PrimNByteInteger 1) = do+            x <- decomp PrimWord8+            return $! toInteger x+        decomp (PrimNByteInteger 2) = do+            x <- decomp PrimWord16+            return $! toInteger x+        decomp (PrimNByteInteger 4) = do+            x <- decomp PrimWord32+            return $! toInteger x+        decomp (PrimNByteInteger 8) = do+            x <- decomp PrimWord64+            return $! toInteger x+        decomp (PrimNByteInteger (n+8))  = do+            x <- decomp PrimWord64+            y <- decomp (PrimNByteInteger n)+            return $! (toInteger x `shiftL` (n `shiftL` 3)) .|. y+        decomp (PrimNByteInteger (n+4))  = do+            x <- decomp PrimWord32+            y <- decomp (PrimNByteInteger n)+            return $! (toInteger x `shiftL` (n `shiftL` 3)) .|. y+        decomp (PrimNByteInteger (n+2))  = do+            x <- decomp PrimWord16+            y <- decomp (PrimNByteInteger n)+            return $! (toInteger x `shiftL` (n `shiftL` 3)) .|. y+-- REDUNDANT CASE+--        decomp (PrimNByteInteger (n+1))  = do+--            x <- decomp PrimWord8+--            y <- decomp (PrimNByteInteger n)+--            return $! (toInteger x `shiftL` (n `shiftL` 3)) .|. y+        decomp (PrimNByteInteger _) = return 0+        +        decomp _ = error ("decomposePrimWhere: no supported primitive to satisfy " ++ show requested)

@@ -26,13 +26,14 @@  import Data.Word import Data.Int+import Control.Monad  -- |Names of standard 'Integral' types integralTypes :: [Name] integralTypes = -    [ ''Int,   ''Integer-    , ''Int8,  ''Int16,  ''Int32,  ''Int64-    , ''Word8, ''Word16, ''Word32, ''Word64+    [ ''Integer+    , ''Int,  ''Int8,  ''Int16,  ''Int32,  ''Int64+    , ''Word, ''Word8, ''Word16, ''Word32, ''Word64     ]  -- |Names of standard 'RealFloat' types@@ -69,11 +70,13 @@ -- This code takes those 2 instance declarations and creates identical ones for -- every type named in 'integralTypes'. replicateInstances :: (Monad m, Data t) => Name -> [Name] -> m [t] -> m [t]-replicateInstances standin types decls = do-    decls <- decls-    sequence-        [ everywhereM (mkM (return . replaceName standin t)) dec-        | t <- types-        , dec <- decls-        ]+replicateInstances standin types getDecls = liftM concat $ sequence+    [ do+        decls <- getDecls+        sequence+            [ everywhereM (mkM (return . replaceName standin t)) dec+            | dec <- decls+            ]+    | t <- types+    ] 

@@ -6,19 +6,55 @@ module Data.Random.Internal.Words where  import Foreign-import GHC.IOBase -import Data.Bits-import Data.Word+-- TODO: add a build flag for endianness-invariance, or just find a way+-- to make sure these operations all do the right thing without costing +-- anything extra at runtime -{-# INLINE buildWord #-}+{-# INLINE buildWord16 #-} -- |Build a word out of 8 bytes.  No promises are made regarding the order -- in which the bytes are stuffed.  Note that this means that a 'RandomSource' -- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc., -- may return different random values on different platforms when started  -- with the same seed, depending on the platform's endianness.-buildWord :: Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word64-buildWord b0 b1 b2 b3 b4 b5 b6 b7+buildWord16 :: Word8 -> Word8 -> Word16+buildWord16 b0 b1+    = unsafePerformIO . allocaBytes 2 $ \p -> do+        pokeByteOff p 0 b0+        pokeByteOff p 1 b1+        peek (castPtr p)++{-# INLINE buildWord32 #-}+-- |Build a word out of 8 bytes.  No promises are made regarding the order+-- in which the bytes are stuffed.  Note that this means that a 'RandomSource'+-- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc.,+-- may return different random values on different platforms when started +-- with the same seed, depending on the platform's endianness.+buildWord32 :: Word8 -> Word8 -> Word8 -> Word8 -> Word32+buildWord32 b0 b1 b2 b3+    = unsafePerformIO . allocaBytes 4 $ \p -> do+        pokeByteOff p 0 b0+        pokeByteOff p 1 b1+        pokeByteOff p 2 b2+        pokeByteOff p 3 b3+        peek (castPtr p)++{-# INLINE buildWord32' #-}+buildWord32' :: Word16 -> Word16 -> Word32+buildWord32' w0 w1+    = unsafePerformIO . allocaBytes 4 $ \p -> do+        pokeByteOff p 0 w0+        pokeByteOff p 2 w1+        peek (castPtr p)++{-# INLINE buildWord64 #-}+-- |Build a word out of 8 bytes.  No promises are made regarding the order+-- in which the bytes are stuffed.  Note that this means that a 'RandomSource'+-- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc.,+-- may return different random values on different platforms when started +-- with the same seed, depending on the platform's endianness.+buildWord64 :: Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word64+buildWord64 b0 b1 b2 b3 b4 b5 b6 b7     = unsafePerformIO . allocaBytes 8 $ \p -> do         pokeByteOff p 0 b0         pokeByteOff p 1 b1@@ -30,6 +66,36 @@         pokeByteOff p 7 b7         peek (castPtr p) +{-# INLINE buildWord64' #-}+buildWord64' :: Word16 -> Word16 -> Word16 -> Word16 -> Word64+buildWord64' w0 w1 w2 w3+    = unsafePerformIO . allocaBytes 8 $ \p -> do+        pokeByteOff p 0 w0+        pokeByteOff p 2 w1+        pokeByteOff p 4 w2+        pokeByteOff p 6 w3+        peek (castPtr p)++{-# INLINE buildWord64'' #-}+buildWord64'' :: Word32 -> Word32 -> Word64+buildWord64'' w0 w1+    = unsafePerformIO . allocaBytes 8 $ \p -> do+        pokeByteOff p 0 w0+        pokeByteOff p 4 w1+        peek (castPtr p)++{-# INLINE word32ToFloat #-}+-- |Pack the low 23 bits from a 'Word32' into a 'Float' in the range [0,1).+-- Used to convert a 'stdUniform' 'Word32' to a 'stdUniform' 'Double'.+word32ToFloat :: Word32 -> Float+word32ToFloat x = (encodeFloat $! toInteger (x .&. 0x007fffff {- 2^23-1 -} )) $ (-23)++{-# INLINE word32ToFloatWithExcess #-}+-- |Same as word32ToFloat, but also return the unused bits (as the 9+-- least significant bits of a 'Word32')+word32ToFloatWithExcess :: Word32 -> (Float, Word32)+word32ToFloatWithExcess x = (word32ToFloat x, x `shiftR` 23)+ {-# INLINE wordToFloat #-} -- |Pack the low 23 bits from a 'Word64' into a 'Float' in the range [0,1). -- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'.@@ -47,6 +113,12 @@ -- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'. wordToDouble :: Word64 -> Double wordToDouble x = (encodeFloat $! toInteger (x .&. 0x000fffffffffffff {- 2^52-1 -})) $ (-52)++{-# INLINE word32ToDouble #-}+-- |Pack a 'Word32' into a 'Double' in the range [0,1).  Note that a Double's +-- mantissa is 52 bits, so this does not fill all of them.+word32ToDouble :: Word32 -> Double+word32ToDouble x = (encodeFloat $! toInteger x) $ (-32)  {-# INLINE wordToDoubleWithExcess #-} -- |Same as wordToDouble, but also return the unused bits (as the 12

@@ -13,6 +13,12 @@ -- For instances where 'm' and 'n' have 'return'/'pure' defined, -- these instances must satisfy -- @lift (return x) == return x@.+-- +-- This form of 'lift' has an extremely general type and is used primarily to+-- support 'sample'.  Its excessive generality is the main reason it's not+-- exported from "Data.Random".  'RVarT' is, however, an instance of +-- 'T.MonadTrans', which in most cases is the preferred way+-- to do the lifting. class Lift m n where     lift :: m a -> n a 

@@ -24,16 +24,20 @@          return (SRS.shuffle xs (reverse is)) --- | A random variable that shuffles a list of a known length. Useful for --- shuffling large lists when the length is known in advance.--- Avoids needing to traverse the list to discover its length.  Each ordering--- has equal probability.------ Throws an error the list is not exactly as long as claimed.+-- | A random variable that shuffles a list of a known length (or a list+-- prefix of the specified length). Useful for shuffling large lists when +-- the length is known in advance.  Avoids needing to traverse the list to+-- discover its length.  Each ordering has equal probability. shuffleN :: Int -> [a] -> RVar [a]-shuffleN 0 xs = return []-shuffleN m@(n+1) xs = do-    is <- sequence [uniform 0 i | i <- [n,n-1..1]]-    return (SRS.shuffle xs is)+shuffleN n xs = shuffleNofM n n xs -    +-- | A random variable that selects N arbitrary elements of a list of known length M.+shuffleNofM :: Int -> Int -> [a] -> RVar [a]+shuffleNofM 0 _ _ = return []+shuffleNofM n m xs+    | n > m     = error "shuffleNofM: n > m"+    | otherwise = do+        is <- sequence [uniform 0 i | i <- take n [m-1, m-2 ..1]]+        return (take n $ SRS.shuffle (take m xs) is)+shuffleNofM _ _ _ = error "shuffleNofM: negative length specified"+

@@ -18,68 +18,183 @@     , runRVar     , RVarT     , runRVarT-    , nByteInteger-    , nBitInteger+    , runRVarTWith     ) where  +import Data.Random.Internal.Primitives import Data.Random.Source import Data.Random.Lift as L -import Data.Bits- import qualified Control.Monad.Trans as T import Control.Applicative-import Control.Monad.Reader import Control.Monad.Identity+import Control.Monad.Prompt (PromptT, runPromptT, prompt) --- |An opaque type containing a \"random variable\" - a value --- which depends on the outcome of some random process.+-- |An opaque type modeling a \"random variable\" - a value +-- which depends on the outcome of some random event.  'RVar's +-- can be conveniently defined by an imperative-looking style:+-- +-- > normalPair =  do+-- >     u <- stdUniform+-- >     t <- stdUniform+-- >     let r = sqrt (-2 * log u)+-- >         theta = (2 * pi) * t+-- >         +-- >         x = r * cos theta+-- >         y = r * sin theta+-- >     return (x,y)+-- +-- OR by a more applicative style:+-- +-- > logNormal = exp <$> stdNormal+--+-- Once defined (in any style), there are a couple ways to sample 'RVar's:+-- +-- * In a monad, using a 'RandomSource':+-- +-- > sampleFrom DevRandom (uniform 1 100) :: IO Int+-- +-- * As a pure function transforming a functional RNG:+-- +-- > sampleState (uniform 1 100) :: StdGen -> (Int, StdGen) type RVar = RVarT Identity --- | single combined container allowing all the relevant --- dictionaries (plus the RandomSource item itself) to be passed--- with one pointer.-data RVarDict n m where-    RVarDict :: (Lift n m, Monad m, RandomSource m s) => s -> RVarDict n m-+-- |\"Run\" an 'RVar' - samples the random variable from the provided+-- source of entropy. runRVar :: RandomSource m s => RVar a -> s -> m a runRVar = runRVarT  -- |A random variable with access to operations in an underlying monad.  Useful -- examples include any form of state for implementing random processes with hysteresis, -- or writer monads for implementing tracing of complicated algorithms.-newtype RVarT n a = RVarT { unRVarT :: forall m r. (a -> m r) -> RVarDict n m -> m r }+-- +-- For example, a simple random walk can be implemented as an 'RVarT' 'IO' value:+--+-- > rwalkIO :: IO (RVarT IO Double)+-- > rwalkIO d = do+-- >     lastVal <- newIORef 0+-- >     +-- >     let x = do+-- >             prev    <- lift (readIORef lastVal)+-- >             change  <- rvarT StdNormal+-- >             +-- >             let new = prev + change+-- >             lift (writeIORef lastVal new)+-- >             return new+-- >         +-- >     return x+--+-- To run the random walk, it must first be initialized, and then it can be sampled as usual:+--+-- > do+-- >     rw <- rwalkIO+-- >     x <- sampleFrom DevURandom rw+-- >     y <- sampleFrom DevURandom rw+-- >     ...+--+-- The same random-walk process as above can be implemented using MTL types+-- as follows (using @import Control.Monad.Trans as MTL@):+-- +-- > rwalkState :: RVarT (State Double) Double+-- > rwalkState = do+-- >     prev <- MTL.lift get+-- >     change  <- rvarT StdNormal+-- >     +-- >     let new = prev + change+-- >     MTL.lift (put new)+-- >     return new+-- +-- Invocation is straightforward (although a bit noisy) if you're used +-- to MTL, but there is a gotcha lurking here: @sample@ and 'runRVarT' +-- inherit the extreme generality of 'lift', so there will almost always+-- need to be an explicit type signature lurking somewhere in any client +-- code making use of 'RVarT' with MTL types.  In this example, the +-- inferred type of @start@ would be too general to be practical, so the+-- signature for @rwalk@  explicitly fixes it to 'Double'.  Alternatively, +-- in this case @sample@ could be replaced with+-- @\\x -> runRVarTWith MTL.lift x StdRandom@.+-- +-- > rwalk :: Int -> Double -> StdGen -> ([Double], StdGen)+-- > rwalk count start gen = evalState (runStateT (sample (replicateM count rwalkState)) gen) start+newtype RVarT m a = RVarT { unRVarT :: PromptT Prim m a } --- | \"Runs\" the monad.+-- | \"Runs\" an 'RVarT', sampling the random variable it defines.+-- +-- The 'Lift' context allows random variables to be defined using a minimal+-- underlying functor ('Identity' is sufficient for \"conventional\" random+-- variables) and then sampled in any monad into which the underlying functor +-- can be embedded (which, for 'Identity', is all monads).+-- +-- The lifting is very important - without it, every 'RVar' would have+-- to either be given access to the full capability of the monad in which it+-- will eventually be sampled (which, incidentally, would also have to be +-- monomorphic so you couldn't sample one 'RVar' in more than one monad)+-- or functions manipulating 'RVar's would have to use higher-ranked +-- types to enforce the same kind of isolation and polymorphism.+-- +-- For non-standard liftings or those where you would rather not introduce a+-- 'Lift' instance, see 'runRVarTWith'.+{-# INLINE runRVarT #-} runRVarT :: (Lift n m, RandomSource m s) => RVarT n a -> s -> m a-runRVarT (RVarT m) (src) = m return (RVarDict src)+runRVarT (RVarT m) src = runPromptT return bindP bindN m+    where+        bindP prim cont = getRandomPrimFrom src prim >>= cont+        bindN nExp cont = lift nExp >>= cont +-- |Like 'runRVarT' but allowing a user-specified lift operation.  This +-- operation must obey the \"monad transformer\" laws:+--+-- > lift . return = return+-- > lift (x >>= f) = (lift x) >>= (lift . f)+--+-- One example of a useful non-standard lifting would be one that takes @State s@ to+-- another monad with a different state representation (such as @IO@ with the+-- state mapped to an @IORef@):+--+-- > embedState :: (Monad m) => m s -> (s -> m ()) -> State s a -> m a+-- > embedState get put = \m -> do+-- >     s <- get+-- >     (res,s) <- return (runState m s)+-- >     put s+-- >     return res+{-# INLINE runRVarTWith #-}+runRVarTWith :: (RandomSource m s) => (forall t. n t -> m t) -> RVarT n a -> s -> m a+runRVarTWith liftN (RVarT m) src = runPromptT return bindP bindN m+    where+        bindP prim cont = getRandomPrimFrom src prim >>= cont+        bindN nExp cont = liftN nExp >>= cont+ instance Functor (RVarT n) where     fmap = liftM  instance Monad (RVarT n) where-    return x = RVarT $ \k _ -> k x-    fail s   = RVarT $ \_ (RVarDict _) -> fail s-    (RVarT m) >>= k = RVarT $ \c s -> m (\a -> unRVarT (k a) c s) s+    return x = RVarT (return $! x)+    fail s   = RVarT (fail s)+    (RVarT m) >>= k = RVarT (m >>= \x -> x `seq` unRVarT (k x))  instance Applicative (RVarT n) where     pure  = return     (<*>) = ap  instance T.MonadTrans RVarT where-    lift m = RVarT $ \k r@(RVarDict _) -> L.lift m >>= \a -> k a+    lift m = RVarT (T.lift m)  instance Lift (RVarT Identity) (RVarT m) where-    lift (RVarT m) = RVarT $ \k (RVarDict src) -> m k (RVarDict src)+    lift (RVarT m) = RVarT (runPromptT return bindP bindN m)+        where+            bindP prim  cont = prompt prim >>= cont+            bindN idExp cont = cont (runIdentity idExp) -instance MonadIO m => MonadIO (RVarT m) where-    liftIO = T.lift . liftIO+instance T.MonadIO m => T.MonadIO (RVarT m) where+    liftIO = T.lift . T.liftIO  instance MonadRandom (RVarT n) where-    getRandomByte   = RVarT $ \k (RVarDict s) -> getRandomByteFrom   s >>= \a -> k a-    getRandomWord   = RVarT $ \k (RVarDict s) -> getRandomWordFrom   s >>= \a -> k a-    getRandomDouble = RVarT $ \k (RVarDict s) -> getRandomDoubleFrom s >>= \a -> k a+    supportedPrims _ _ = True+    {-# INLINE getSupportedRandomPrim #-}+    getSupportedRandomPrim p    = RVarT (prompt p)+    {-# INLINE getRandomPrim #-}+    getRandomPrim p = RVarT (prompt p)  -- I would really like to be able to do this, but I can't because of the -- blasted Eq and Show in Num's class context...@@ -91,32 +206,3 @@ --     signum = liftA signum --     abs = liftA abs --     fromInteger = pure . fromInteger---- some 'fundamental' RVarTs--- this maybe ought to even be a part of the RandomSource class...-{-# INLINE nByteInteger #-}--- |A random variable evenly distributed over all unsigned integers from--- 0 to 2^(8*n)-1, inclusive.-nByteInteger :: Int -> RVarT m Integer-nByteInteger 1 = do-    x <- getRandomByte-    return $! toInteger x-nByteInteger 8 = do-    x <- getRandomWord-    return $! toInteger x-nByteInteger n = nBitInteger (n `shiftL` 3)--{-# INLINE nBitInteger #-}--- |A random variable evenly distributed over all unsigned integers from--- 0 to 2^n-1, inclusive.-nBitInteger :: Int -> RVarT m Integer-nBitInteger 8  = do-    x <- getRandomByte-    return $! toInteger x-nBitInteger (n+64) = do-    x <- getRandomWord-    y <- nBitInteger n-    return $! (toInteger x `shiftL` n) .|. y-nBitInteger n = do-        x <- getRandomWord-        return $! toInteger (x `shiftR` (64-n))

@@ -9,6 +9,7 @@  module Data.Random.Sample where +import Control.Monad.State  import Data.Random.Distribution import Data.Random.Lift import Data.Random.RVar@@ -17,7 +18,8 @@  -- |A typeclass allowing 'Distribution's and 'RVar's to be sampled.  Both may -- also be sampled via 'runRVar' or 'runRVarT', but I find it psychologically--- pleasing to be able to sample both using this function.+-- pleasing to be able to sample both using this function, as they are two+-- separate abstractions for one base concept: a random variable. class Sampleable d m t where     -- |Directly sample from a distribution or random variable, using the given source of entropy.     sampleFrom :: RandomSource m s => s -> d t -> m t@@ -29,7 +31,17 @@ instance Lift m n => Sampleable (RVarT m) n t where     sampleFrom src x = runRVarT x src --- |Sample a distribution using the default source of entropy for the+-- |Sample a random variable using the default source of entropy for the -- monad in which the sampling occurs. sample :: (Sampleable d m t, MonadRandom m) => d t -> m t sample = sampleFrom StdRandom++-- |Sample a random variable in a \"functional\" style.  Typical instantiations+-- of @s@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.+sampleState :: (Sampleable d (State s) t, MonadRandom (State s)) => d t -> s -> (t, s)+sampleState thing = runState (sample thing)++-- |Sample a random variable in a \"semi-functional\" style.  Typical instantiations+-- of @s@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.+sampleStateT :: (Sampleable d (StateT s m) t, MonadRandom (StateT s m)) => d t -> s -> m (t, s)+sampleStateT thing = runStateT (sample thing)

@@ -2,18 +2,20 @@  -      ``Data/Random/Source''  -} {-# LANGUAGE-    MultiParamTypeClasses, FlexibleInstances+    MultiParamTypeClasses, FlexibleInstances, GADTs   #-}  module Data.Random.Source     ( MonadRandom(..)     , RandomSource(..)+    , Prim(..)     ) where  import Data.Word-import Control.Monad+import Control.Monad.Prompt+import Data.Tagged -import Data.Random.Internal.Words+import Data.Random.Internal.Primitives  -- |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@@ -21,68 +23,119 @@ -- when directly requesting entropy for a random variable these functions -- are used. -- --- The minimal definition is either 'getRandomByte' or 'getRandomWord'.--- 'getRandomDouble' is defaulted in terms of 'getRandomWord'.+-- The minimal definition is 'supportedPrims' and 'getSupportedRandomPrim'+-- with cases for those primitives where 'supportedPrims' returns 'True'.+--+-- It is recommended (despite the warnings it generates) that, even when+-- all primitives are supported, a final wildcard case of 'supportedPrims' is+-- specified, as:+-- +-- > supportedPrims _ _ = False+--+-- The overlapping pattern warnings can be suppressed (without suppressing +-- other, genuine, overlapping-pattern warnings) by the GHC flag+-- @-fno-warn-simple-patterns@.  This is not actually the documented behavior+-- of that flag as far as I can find in 3 google-minutes, but it works with+-- GHC 6.12.1 anyway, and that's good enough for me.+--+-- Note that it is very important that at least 'supportedPrims' (and preferably+-- 'getSupportedRandomPrim' as well) gets inlined into the default implementation+-- of 'getRandomPrim'.  If your 'supportedPrims' is more than about 2 or 3+-- cases, add an INLINE pragma so that it can be optimized out of 'getRandomPrim'. class Monad m => MonadRandom m where-    -- |Get a random uniformly-distributed byte.-    getRandomByte :: m Word8-    getRandomByte = do-        word <- getRandomWord-        return (fromIntegral word)-    -    -- |Get a random 'Word64' uniformly-distributed over the full range of the type.-    getRandomWord :: m Word64-    getRandomWord = do-        b0 <- getRandomByte-        b1 <- getRandomByte-        b2 <- getRandomByte-        b3 <- getRandomByte-        b4 <- getRandomByte-        b5 <- getRandomByte-        b6 <- getRandomByte-        b7 <- getRandomByte-        -        return (buildWord b0 b1 b2 b3 b4 b5 b6 b7)+    -- |Predicate indicating whether a given primitive is supported by the+    -- instance.  The first parameter is a phantom used to select the instance.+    supportedPrims :: m () -> Prim t -> Bool     -    -- |Get a random 'Double' uniformly-distributed over the interval [0,1)-    getRandomDouble :: m Double-    getRandomDouble = do-        word <- getRandomWord-        return (wordToDouble word)+    -- |Generate a random value corresponding to the specified primitive.  Will+    -- not be called unless supportedPrims returns true for that primitive.+    getSupportedRandomPrim :: Prim t -> m t +    -- This could just be a function, but placing it in a dictionary gives+    -- GHC a place to optimize it separately for each instance, which is +    -- kinda the whole point of the 'Prim' machinery:+    -- +    -- |Generate a random value corresponding to the specified primitive.  The+    -- default implementation makes use of 'supportedPrims' and 'getSupportedRandomPrim'+    -- to construct any required Prim out of the supported ones.+    {-# NOINLINE getRandomPrim #-}+    getRandomPrim :: Prim t -> m t+    getRandomPrim prim = val+        where+            val = runPromptM getSupportedRandomPrim (decomposePrimWhere (supportedPrims mPhantom) prim)+            mPhantom = error "supportedPrims tried to evaluate a phantom parameter" `asTypeOf` (val >> return ())+ -- |A source of entropy which can be used in the given monad. ----- The minimal definition is either 'getRandomByteFrom' or 'getRandomWordFrom'.--- 'getRandomDoubleFrom' is defaulted in terms of 'getRandomWordFrom'+-- The minimal definition is 'supportedPrimsFrom' and 'getSupportedRandomPrimFrom'+-- with cases for those primitives where 'supportedPrimsFrom' returns 'True'.+-- +-- Note that it is very important that at least 'supportedPrimsFrom' (and preferably+-- 'getSupportedRandomPrimFrom' as well) gets inlined into the default implementation+-- of 'getRandomPrimFrom'.  If your 'supportedPrimsFrom' is more than about 2 or 3+-- cases, add an INLINE pragma so that it can be optimized out of 'getRandomPrimFrom'.+-- +-- See also 'MonadRandom'. class Monad m => RandomSource m s where-    -- |Get a random uniformly-distributed byte.-    getRandomByteFrom :: s -> m Word8-    getRandomByteFrom src = do-        word <- getRandomWordFrom src-        return (fromIntegral word)+    -- |Predicate indicating whether a given primitive is supported by the+    -- instance.  The tag on the first parameter is a phantom used only to+    -- select the instance, but the value itself may be inspected.+    supportedPrimsFrom :: Tagged (m ()) s -> Prim t -> Bool     -    -- |Get a random 'Word64' uniformly-distributed over the full range of the type.-    getRandomWordFrom :: s -> m Word64-    getRandomWordFrom src = do-        b0 <- getRandomByteFrom src-        b1 <- getRandomByteFrom src-        b2 <- getRandomByteFrom src-        b3 <- getRandomByteFrom src-        b4 <- getRandomByteFrom src-        b5 <- getRandomByteFrom src-        b6 <- getRandomByteFrom src-        b7 <- getRandomByteFrom src-        -        return (buildWord b0 b1 b2 b3 b4 b5 b6 b7)+    -- |Generate a random value corresponding to the specified primitive+    getSupportedRandomPrimFrom :: s -> Prim t -> m t     -    -- |Get a random 'Double' uniformly-distributed over the interval [0,1)-    getRandomDoubleFrom :: s -> m Double-    getRandomDoubleFrom src = do-        word <- getRandomWordFrom src-        return (wordToDouble word)+    +    -- This could just be a function, but placing it in a dictionary gives+    -- GHC a place to optimize it separately for each instance, which is +    -- kinda the whole point of the 'Prim' machinery:+    -- +    -- |Generate a random value corresponding to the specified primitive.  The+    -- default implementation makes use of 'supportedPrimsFrom' and+    -- 'getSupportedRandomPrimFrom' to construct any required Prim out of +    -- the supported ones.+    {-# NOINLINE getRandomPrimFrom #-}+    getRandomPrimFrom :: s -> Prim t -> m t+    getRandomPrimFrom src prim = val+        where+            val = runPromptM (getSupportedRandomPrimFrom src) (decomposePrimWhere supported prim)+            supported :: Prim t -> Bool+            supported = supportedPrimsFrom (tagIt (val >> return ()) src)+            +            tagIt :: a -> b -> Tagged a b+            tagIt _ it = Tagged it  instance Monad m => RandomSource m (m Word8) where-    getRandomByteFrom = id+    supportedPrimsFrom _ PrimWord8 = True+    supportedPrimsFrom _ _ = False+    +    getSupportedRandomPrimFrom f PrimWord8 = f+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word8): unsupported prim requested: " ++ show p) +instance Monad m => RandomSource m (m Word16) where+    supportedPrimsFrom _ PrimWord16 = True+    supportedPrimsFrom _ _ = False+    +    getSupportedRandomPrimFrom f PrimWord16 = f+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word16): unsupported prim requested: " ++ show p)++instance Monad m => RandomSource m (m Word32) where+    supportedPrimsFrom _ PrimWord32 = True+    supportedPrimsFrom _ _ = False+    +    getSupportedRandomPrimFrom f PrimWord32 = f+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word32): unsupported prim requested: " ++ show p)+ instance Monad m => RandomSource m (m Word64) where-    getRandomWordFrom = id+    supportedPrimsFrom _ PrimWord64 = True+    supportedPrimsFrom _ _ = False+    +    getSupportedRandomPrimFrom f PrimWord64 = f+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word64): unsupported prim requested: " ++ show p)++instance Monad m => RandomSource m (m Double) where+    supportedPrimsFrom _ PrimDouble = True+    supportedPrimsFrom _ _ = False+    +    getSupportedRandomPrimFrom f PrimDouble = f+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Double): unsupported prim requested: " ++ show p)

@@ -2,7 +2,7 @@  -      ``Data/Random/Source/DevRandom''  -} {-# LANGUAGE-    MultiParamTypeClasses+    MultiParamTypeClasses, GADTs   #-}  module Data.Random.Source.DevRandom @@ -11,7 +11,7 @@  import Data.Random.Source -import System.IO (openBinaryFile, hGetBuf, IOMode(..))+import System.IO (openBinaryFile, hGetBuf, Handle, IOMode(..)) import Foreign  -- |On systems that have it, \/dev\/random is a handy-dandy ready-to-use source@@ -22,19 +22,32 @@ -- purposes other than cryptography, \/dev\/urandom is preferable because when it -- runs out of real entropy it'll still churn out pseudorandom data. data DevRandom = DevRandom | DevURandom+    deriving (Eq, Show)  {-# NOINLINE devRandom  #-}+devRandom :: Handle devRandom  = unsafePerformIO (openBinaryFile "/dev/random"  ReadMode) {-# NOINLINE devURandom #-}+devURandom :: Handle devURandom = unsafePerformIO (openBinaryFile "/dev/urandom" ReadMode) +dev :: DevRandom -> Handle dev DevRandom  = devRandom dev DevURandom = devURandom  instance RandomSource IO DevRandom where-    getRandomByteFrom src  = allocaBytes 1 $ \buf -> do+    supportedPrimsFrom _ PrimWord8          = True+    supportedPrimsFrom _ PrimWord32         = True+    supportedPrimsFrom _ PrimWord64         = True+    supportedPrimsFrom _ _ = False+    +    getSupportedRandomPrimFrom src PrimWord8  = allocaBytes 1 $ \buf -> do         1 <- hGetBuf (dev src) buf  1         peek buf-    getRandomWordFrom src  = allocaBytes 8 $ \buf -> do+    getSupportedRandomPrimFrom src PrimWord32  = allocaBytes 1 $ \buf -> do+        4 <- hGetBuf (dev src) buf  4+        peek (castPtr buf)+    getSupportedRandomPrimFrom src PrimWord64  = allocaBytes 8 $ \buf -> do         8 <- hGetBuf (dev src) buf  8         peek (castPtr buf)+    getSupportedRandomPrimFrom src prim = error ("getSupportedRandomPrimFrom/" ++ show src ++ ": unsupported prim requested: " ++ show prim)

@@ -0,0 +1,46 @@+{-# LANGUAGE+        MultiParamTypeClasses,+        FlexibleInstances,+        GADTs+  #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}+-- |This module defines the following instances:+-- +-- > instance RandomSource (ST s) (Gen s)+-- > instance RandomSource IO (Gen RealWorld)+module Data.Random.Source.MWC where++import Data.Random.Internal.Primitives+import Data.Random.Internal.Words+import Data.Random.Source+import System.Random.MWC+import Control.Monad.ST++instance RandomSource (ST s) (Gen s) where+    {-# INLINE supportedPrimsFrom #-}+    supportedPrimsFrom _ PrimWord8  = True+    supportedPrimsFrom _ PrimWord16 = True+    supportedPrimsFrom _ PrimWord32 = True+    supportedPrimsFrom _ PrimWord64 = True+    supportedPrimsFrom _ PrimDouble = True+    supportedPrimsFrom _ _ = False+    +    {-# INLINE getSupportedRandomPrimFrom #-}+    getSupportedRandomPrimFrom gen PrimWord8    = uniform gen+    getSupportedRandomPrimFrom gen PrimWord16   = uniform gen+    getSupportedRandomPrimFrom gen PrimWord32   = uniform gen+    getSupportedRandomPrimFrom gen PrimWord64   = uniform gen+    getSupportedRandomPrimFrom gen PrimDouble   = fmap wordToDouble (uniform gen)+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/Gen s: unsupported prim requested: " ++ show p)++instance RandomSource IO (Gen RealWorld) where+    {-# INLINE supportedPrimsFrom #-}+    supportedPrimsFrom _ PrimWord8  = True+    supportedPrimsFrom _ PrimWord16 = True+    supportedPrimsFrom _ PrimWord32 = True+    supportedPrimsFrom _ PrimWord64 = True+    supportedPrimsFrom _ PrimDouble = True+    supportedPrimsFrom _ _ = False+    +    {-# INLINE getSupportedRandomPrimFrom #-}+    getSupportedRandomPrimFrom gen prim = stToIO (getSupportedRandomPrimFrom gen prim)

@@ -1,8 +1,11 @@ {-# LANGUAGE+    BangPatterns,     MultiParamTypeClasses,     FlexibleContexts, FlexibleInstances,-    UndecidableInstances+    UndecidableInstances,+    GADTs   #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}  -- |This module provides functions useful for implementing new 'MonadRandom' -- and 'RandomSource' instances for state-abstractions containing 'PureMT'@@ -11,52 +14,38 @@ -- cases. module Data.Random.Source.PureMT where -import Data.Random.Internal.Words+import Data.Random.Internal.Primitives import Data.Random.Source import System.Random.Mersenne.Pure64  import Data.StateRef-import Data.Word +import Control.Monad.Prompt import Control.Monad.State import qualified Control.Monad.ST.Strict as S import qualified Control.Monad.State.Strict as S --- |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@, @getRandomWordFromMTRef x@ can be--- used as a 'RandomSource' in 'IO', 'STM', or any monad which is an instance--- of 'MonadIO'.  These functions can also be used to implement additional--- 'RandomSource' instances for mutable references to 'PureMT' states.-getRandomWordFromMTRef :: ModifyRef sr m PureMT => sr -> m Word64-getRandomWordFromMTRef ref = do-    atomicModifyReference ref (swap . randomWord64)-    -    where-        swap (a,b) = (b,a)--getRandomByteFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> m Word8-getRandomByteFromMTRef ref = do-    x <- atomicModifyReference ref (swap . randomInt)-    return (fromIntegral x)-    -    where-        swap (a,b) = (b,a)---- for whatever reason, my simple wordToDouble is faster than whatever--- the mersenne random library is using, at least in the version I have.--- if this changes, switch to the commented version.--- Same thing below, in getRandomDoubleFromMTState.-getRandomDoubleFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> m Double-getRandomDoubleFromMTRef src = liftM wordToDouble (getRandomWordFromMTRef src)--- getRandomDoubleFromMTRef ref = do---     atomicModifyReference ref (swap . randomDouble)---     ---     where---         swap (a,b) = (b,a)+-- |Given a mutable reference to a 'PureMT' generator, we can implement+-- 'RandomSource' for in any monad in which the reference can be modified.+getRandomPrimFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> Prim a -> m a+getRandomPrimFromMTRef ref prim+    | supported prim = getThing (genPrim prim)+    | otherwise = runPromptM (getRandomPrimFromMTRef ref) (decomposePrimWhere supported prim)+    where +        supported :: Prim a -> Bool+        supported PrimWord64 = True+        supported PrimDouble = True+        supported _          = False+        +        genPrim :: Prim a -> (PureMT -> (a, PureMT))+        genPrim PrimWord64 = randomWord64+        genPrim PrimDouble = randomDouble+        genPrim p = error ("getRandomPrimFromMTRef: genPrim called for unsupported prim " ++ show p)+        +        getThing thing = atomicModifyReference ref $ \(!oldMT) -> case thing oldMT of (!w, !newMT) -> (newMT, w)+             --- |Similarly, @getRandomWordFromMTState x@ can be used in any \"state\"+-- |Similarly, @getRandomPrimFromMTState 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@@ -64,71 +53,65 @@ -- @runState . sample :: Distribution d t => d t -> PureMT -> (t, PureMT)@. -- 'PureMT' in the type there can be replaced by 'StdGen' or anything else  -- satisfying @MonadRandom (State s) => s@).-getRandomWordFromMTState :: MonadState PureMT m => m Word64-getRandomWordFromMTState = do-    mt <- get-    let (ws, newMt) = randomWord64 mt-    put newMt-    return ws--getRandomByteFromMTState :: MonadState PureMT m => m Word8-getRandomByteFromMTState = do-    mt <- get-    let (ws, newMt) = randomInt mt-    put newMt-    return (fromIntegral ws)---getRandomDoubleFromMTState :: MonadState PureMT m => m Double-getRandomDoubleFromMTState = liftM wordToDouble getRandomWordFromMTState--- getRandomDoubleFromMTState = do---     mt <- get---     let (x, newMt) = randomDouble mt---     put newMt---     return x+getRandomPrimFromMTState :: MonadState PureMT m => Prim a -> m a+getRandomPrimFromMTState prim+    | supported prim = getThing (genPrim prim)+    | otherwise = runPromptM getRandomPrimFromMTState (decomposePrimWhere supported prim)+    where+        supported :: Prim a -> Bool+        supported PrimWord64 = True+        supported PrimDouble = True+        supported _          = False+        +        genPrim :: Prim a -> (PureMT -> (a, PureMT))+        genPrim PrimWord64 = randomWord64+        genPrim PrimDouble = randomDouble+        genPrim p = error ("getRandomPrimFromMTRef: genPrim called for unsupported prim " ++ show p)+        +        getThing thing = do+            !mt <- get+            let (!ws, !newMt) = thing mt+            put newMt+            return ws  instance MonadRandom (State PureMT) where-    getRandomByte   = getRandomByteFromMTState-    getRandomWord   = getRandomWordFromMTState-    getRandomDouble = getRandomDoubleFromMTState+    supportedPrims _ _ = True+    getSupportedRandomPrim = getRandomPrimFromMTState  instance MonadRandom (S.State PureMT) where-    getRandomByte   = getRandomByteFromMTState-    getRandomWord   = getRandomWordFromMTState-    getRandomDouble = getRandomDoubleFromMTState+    supportedPrims _ _ = True+    getSupportedRandomPrim = getRandomPrimFromMTState  instance (Monad m1, ModifyRef (Ref m2 PureMT) m1 PureMT) => RandomSource m1 (Ref m2 PureMT) where-    getRandomByteFrom   = getRandomByteFromMTRef-    getRandomWordFrom   = getRandomWordFromMTRef-    getRandomDoubleFrom = getRandomDoubleFromMTRef-+    supportedPrimsFrom _ _ = True+    getSupportedRandomPrimFrom = getRandomPrimFromMTRef+     instance Monad m => MonadRandom (StateT PureMT m) where-    getRandomByte   = getRandomByteFromMTState-    getRandomWord   = getRandomWordFromMTState-    getRandomDouble = getRandomDoubleFromMTState+    supportedPrims _ _ = True+    getSupportedRandomPrim = getRandomPrimFromMTState  instance Monad m => MonadRandom (S.StateT PureMT m) where-    getRandomByte   = getRandomByteFromMTState-    getRandomWord   = getRandomWordFromMTState-    getRandomDouble = getRandomDoubleFromMTState+    supportedPrims _ _ = True+    getSupportedRandomPrim = getRandomPrimFromMTState  instance (Monad m, ModifyRef (IORef PureMT) m PureMT) => RandomSource m (IORef PureMT) where     {-# SPECIALIZE instance RandomSource IO (IORef PureMT)#-}-    getRandomByteFrom   = getRandomByteFromMTRef-    getRandomWordFrom   = getRandomWordFromMTRef-    getRandomDoubleFrom = getRandomDoubleFromMTRef-+    supportedPrimsFrom _ _ = True+    getSupportedRandomPrimFrom = getRandomPrimFromMTRef+     instance (Monad m, ModifyRef (STRef s PureMT) m PureMT) => RandomSource m (STRef s PureMT) where     {-# SPECIALIZE instance RandomSource (ST s) (STRef s PureMT) #-}     {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s PureMT) #-}-    getRandomByteFrom   = getRandomByteFromMTRef-    getRandomWordFrom   = getRandomWordFromMTRef-    getRandomDoubleFrom = getRandomDoubleFromMTRef--instance (Monad m, ModifyRef (TVar PureMT) m PureMT) => RandomSource m (TVar PureMT) where-    {-# SPECIALIZE instance RandomSource IO  (TVar PureMT) #-}-    {-# SPECIALIZE instance RandomSource STM (TVar PureMT) #-}-    getRandomByteFrom   = getRandomByteFromMTRef-    getRandomWordFrom   = getRandomWordFromMTRef-    getRandomDoubleFrom = getRandomDoubleFromMTRef+    supportedPrimsFrom _ _ = True+    getSupportedRandomPrimFrom = getRandomPrimFromMTRef +-- Note that this instance is probably a Bad Idea.  STM allows random variables+-- to interact in spooky quantum-esque ways - One transaction can 'retry' until+-- it gets a \"random\" answer it likes, which causes it to selectively consume +-- entropy, biasing the supply from which other random variables will draw.+-- instance (Monad m, ModifyRef (TVar PureMT) m PureMT) => RandomSource m (TVar PureMT) where+--     {-# SPECIALIZE instance RandomSource IO  (TVar PureMT) #-}+--     {-# SPECIALIZE instance RandomSource STM (TVar PureMT) #-}+--     supportedPrimsFrom _ _ = True+--     getSupportedRandomPrimFrom = getRandomPrimFromMTRef+    

@@ -8,6 +8,7 @@ module Data.Random.Source.Std where  import Data.Random.Source+import Data.Tagged  -- |A token representing the \"standard\" entropy source in a 'MonadRandom' -- monad.  Its sole purpose is to make the following true (when the types check):@@ -16,6 +17,11 @@ data StdRandom = StdRandom  instance MonadRandom m => RandomSource m StdRandom where-    getRandomByteFrom   StdRandom = getRandomByte-    getRandomWordFrom   StdRandom = getRandomWord-    getRandomDoubleFrom StdRandom = getRandomDouble+    {-SPECIALIZE instance MonadRandom m => RandomSource m StdRandom -}+    supportedPrimsFrom w = supportedPrims (mkWit w)+        where+            mkWit :: Tagged a b -> a+            mkWit = error "supportedPrims tried to evaluate its phantom parameter"+    getSupportedRandomPrimFrom   StdRandom = getSupportedRandomPrim+    +    getRandomPrimFrom StdRandom = getRandomPrim

@@ -1,9 +1,8 @@-{-- -      ``Data/Random/Source/StdGen''- -} {-# LANGUAGE-    MultiParamTypeClasses, FlexibleInstances, UndecidableInstances+    MultiParamTypeClasses, FlexibleInstances, UndecidableInstances, GADTs,+    BangPatterns, RankNTypes   #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}  -- |This module provides functions useful for implementing new 'MonadRandom' -- and 'RandomSource' instances for state-abstractions containing 'StdGen'@@ -13,131 +12,165 @@ module Data.Random.Source.StdGen where  import Data.Random.Internal.Words+import Data.Random.Internal.Primitives import Data.Random.Source import System.Random+import Control.Monad.Prompt import Control.Monad.State import qualified Control.Monad.ST.Strict as S import qualified Control.Monad.State.Strict as S import Data.StateRef import Data.Word + instance (Monad m1, ModifyRef (Ref m2 StdGen) m1 StdGen) => RandomSource m1 (Ref m2 StdGen) where-    getRandomByteFrom   = getRandomByteFromRandomGenRef-    getRandomWordFrom   = getRandomWordFromRandomGenRef-    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef+    supportedPrimsFrom _ _ = True+    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef  instance (Monad m, ModifyRef (IORef   StdGen) m StdGen) => RandomSource m (IORef   StdGen) where     {-# SPECIALIZE instance RandomSource IO (IORef StdGen) #-}-    getRandomByteFrom   = getRandomByteFromRandomGenRef-    getRandomWordFrom   = getRandomWordFromRandomGenRef-    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef-instance (Monad m, ModifyRef (TVar    StdGen) m StdGen) => RandomSource m (TVar    StdGen) where-    {-# SPECIALIZE instance RandomSource IO  (TVar StdGen) #-}-    {-# SPECIALIZE instance RandomSource STM (TVar StdGen) #-}-    getRandomByteFrom   = getRandomByteFromRandomGenRef-    getRandomWordFrom   = getRandomWordFromRandomGenRef-    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef+    supportedPrimsFrom _ _ = True+    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef++-- Note that this instance is probably a Bad Idea.  STM allows random variables+-- to interact in spooky quantum-esque ways - One transaction can 'retry' until+-- it gets a \"random\" answer it likes, which causes it to selectively consume +-- entropy, biasing the supply from which other random variables will draw.+-- instance (Monad m, ModifyRef (TVar    StdGen) m StdGen) => RandomSource m (TVar    StdGen) where+--     {-# SPECIALIZE instance RandomSource IO  (TVar StdGen) #-}+--     {-# SPECIALIZE instance RandomSource STM (TVar StdGen) #-}+--     supportedPrimsFrom _ _ = True+--     getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef+ instance (Monad m, ModifyRef (STRef s StdGen) m StdGen) => RandomSource m (STRef s StdGen) where     {-# SPECIALIZE instance RandomSource (ST s) (STRef s StdGen) #-}     {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s StdGen) #-}-    getRandomByteFrom   = getRandomByteFromRandomGenRef-    getRandomWordFrom   = getRandomWordFromRandomGenRef-    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef--getRandomByteFromStdGenIO :: IO Word8-getRandomByteFromStdGenIO = do-    int <- randomRIO (0, 255) :: IO Int-    return (fromIntegral int)+    supportedPrimsFrom _ _ = True+    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef -getRandomWordFromStdGenIO :: IO Word64-getRandomWordFromStdGenIO = do-    int <- randomRIO (0, 0xffffffffffffffff)-    return (fromInteger int)+getRandomPrimFromStdGenIO :: Prim a -> IO a+getRandomPrimFromStdGenIO prim+    | supported prim = genPrim prim+    | otherwise = runPromptM getRandomPrimFromStdGenIO (decomposePrimWhere supported prim)+    where +        {-# INLINE supported #-}+        supported :: Prim a -> Bool+        supported PrimWord8             = True+        supported PrimWord16            = True+        supported PrimWord32            = True+        supported PrimWord64            = True+        supported PrimDouble            = True+        supported (PrimNByteInteger _)  = True+        supported _                     = False+        +        -- based on reading the source of the "random" library's implementation, I do+        -- not believe that the randomRIO (0,1) implementation for Double is capable of producing+        -- the value 0.  Therefore, I'm not using it.  If this is an incorrect reading on+        -- my part, or if this changes, then feel free to change the implementation.+        -- Same goes for the other getRandomDouble... functions here. --- based on reading the source of the "random" library's implementation, I do--- not believe that the randomRIO (0,1) implementation for Double is capable of producing--- the value 0.  Therefore, I'm not using it.  If this is an incorrect reading on--- my part, or if this changes, then feel free to use the commented version.--- Same goes for the other getRandomDouble... functions here.-getRandomDoubleFromStdGenIO :: IO Double-getRandomDoubleFromStdGenIO = liftM wordToDouble getRandomWordFromStdGenIO--- getRandomDoubleFromStdGenIO = randomRIO (0, 1)+        {-# INLINE genPrim #-}+        genPrim :: Prim a -> IO a+        genPrim PrimWord8            = fmap fromIntegral                  (randomRIO (0, 0xff) :: IO Int)+        genPrim PrimWord16           = fmap fromIntegral                  (randomRIO (0, 0xffff) :: IO Int)+        genPrim PrimWord32           = fmap fromInteger                   (randomRIO (0, 0xffffffff))+        genPrim PrimWord64           = fmap fromInteger                   (randomRIO (0, 0xffffffffffffffff))+        genPrim PrimDouble           = fmap (wordToDouble . fromInteger)  (randomRIO (0, 0xffffffffffffffff))+        genPrim (PrimNByteInteger n) = randomRIO (0, iterate (*256) 1 !! n)+        genPrim p = error ("getRandomPrimFromStdGenIO: genPrim called for unsupported prim " ++ show p)  -- |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@, @getRandomByteFromRandomGenRef 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--- 'getRandomWordFromRandomGenRef' though, as this one is likely to throw--- away a lot of perfectly good entropy.  Better still is to use these 3 functions--- together to create a 'RandomSource' instance for the reference you're using,--- if one does not already exist.-getRandomByteFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>-                                  sr -> m Word8-getRandomByteFromRandomGenRef g = atomicModifyReference g (swap . randomR (0,255))+getRandomPrimFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>+                                  sr -> Prim a -> m a+getRandomPrimFromRandomGenRef ref prim+    | supported prim = genPrim prim getThing+    | otherwise = runPromptM (getRandomPrimFromRandomGenRef ref) (decomposePrimWhere supported prim)     where -        swap :: (Int, a) -> (a, Word8)-        swap (a,b) = (b,fromIntegral a)+        {-# INLINE supported #-}+        supported :: Prim a -> Bool+        supported PrimWord8             = True+        supported PrimWord16            = True+        supported PrimWord32            = True+        supported PrimWord64            = True+        supported PrimDouble            = True+        supported (PrimNByteInteger _)  = True+        supported _                     = False+        +        {-# INLINE genPrim #-}+        genPrim :: (RandomGen g) => Prim a -> (forall b. (g -> (b, g)) -> (b -> a) -> c) -> c+        genPrim PrimWord8            f = f (randomR (0, 0xff))                (fromIntegral :: Int -> Word8)+        genPrim PrimWord16           f = f (randomR (0, 0xffff))              (fromIntegral :: Int -> Word16)+        genPrim PrimWord32           f = f (randomR (0, 0xffffffff))          (fromInteger)+        genPrim PrimWord64           f = f (randomR (0, 0xffffffffffffffff))  (fromInteger)+        genPrim PrimDouble           f = f (randomR (0, 0x000fffffffffffff))  (flip encodeFloat (-52))+        genPrim (PrimNByteInteger n) f = f (randomR (0, iterate (*256) 1 !! n)) (id :: Integer -> Integer)+        genPrim p _ = error ("getRandomPrimFromRandomGenRef: genPrim called for unsupported prim " ++ show p)+        +        {-# INLINE getThing #-}+        getThing thing f = atomicModifyReference ref $ \(!oldMT) -> case thing oldMT of (!w, !newMT) -> (newMT, f w) -getRandomWordFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>-                                  sr -> m Word64-getRandomWordFromRandomGenRef g = atomicModifyReference g (swap . randomR (0,0xffffffffffffffff))-    where swap (a,b) = (b,fromInteger a) -getRandomDoubleFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>-                                  sr -> m Double-getRandomDoubleFromRandomGenRef g = liftM wordToDouble (getRandomWordFromRandomGenRef g)--- getRandomDoubleFromRandomGenRef g = atomicModifyRef g (swap . randomR (0,1))---     where swap (a,b) = (b,a)- -- |Similarly, @getRandomWordFromRandomGenState 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.-getRandomByteFromRandomGenState :: (RandomGen g, MonadState g m) => m Word8-getRandomByteFromRandomGenState = do-    g <- get-    case randomR (0, 255 :: Int) g of-        (i,g) -> do-            put g-            return (fromIntegral i)--getRandomWordFromRandomGenState :: (RandomGen g, MonadState g m) => m Word64-getRandomWordFromRandomGenState = do-    g <- get-    case randomR (0, 0xffffffffffffffff) g of-        (i,g) -> do-            put g-            return (fromInteger i)--getRandomDoubleFromRandomGenState :: (RandomGen g, MonadState g m) => m Double-getRandomDoubleFromRandomGenState = liftM wordToDouble getRandomWordFromRandomGenState--- getRandomDoubleFromRandomGenState = do---     g <- get---     case randomR (0, 1) g of---         (x,g) -> do---             put g---             return x-+{-# SPECIALIZE getRandomPrimFromRandomGenState :: Prim a -> State StdGen a #-}+{-# SPECIALIZE getRandomPrimFromRandomGenState :: Monad m => Prim a -> StateT StdGen m a #-}+getRandomPrimFromRandomGenState :: (RandomGen g, MonadState g m) => Prim a -> m a+getRandomPrimFromRandomGenState prim+    | supported prim = genSupported prim+    | otherwise = runPromptM genSupported (decomposePrimWhere supported prim)+    where +        {-# INLINE genSupported #-}+        genSupported prim = genPrim prim getThing+        +        {-# INLINE supported #-}+        supported :: Prim a -> Bool+        supported PrimWord8             = True+        supported PrimWord16            = True+        supported PrimWord32            = True+        supported PrimWord64            = True+        supported PrimDouble            = True+        supported (PrimNByteInteger _)  = True+        supported _                     = False+        +        {-# INLINE genPrim #-}+        genPrim :: (RandomGen g) => Prim a -> (forall b. (g -> (b, g)) -> (b -> a) -> c) -> c+        genPrim PrimWord8            f = f (randomR (0, 0xff))                (fromIntegral :: Int -> Word8)+        genPrim PrimWord16           f = f (randomR (0, 0xffff))              (fromIntegral :: Int -> Word16)+        genPrim PrimWord32           f = f (randomR (0, 0xffffffff))          (fromInteger)+        genPrim PrimWord64           f = f (randomR (0, 0xffffffffffffffff))  (fromInteger)+        genPrim PrimDouble           f = f (randomR (0, 0x000fffffffffffff))  (flip encodeFloat (-52))+          {- not using the Random Double instance for 2 reasons.  1st, it only generates 32 bits of entropy, when +             a [0,1) Double has room for 52.  Second, it appears there's a bug where it can actually generate a +             negative number in the case where randomIvalInteger returns minBound::Int32. -}+--        genPrim PrimDouble f = f (randomR (0, 1.0))  (id)+        genPrim (PrimNByteInteger n) f = f (randomR (0, iterate (*256) 1 !! n)) id+        genPrim p _ = error ("getRandomPrimFromRandomGenState: genPrim called for unsupported prim " ++ show p)+        +        {-# INLINE getThing #-}+        getThing thing f = do+            !oldGen <- get+            case thing oldGen of+                (!i,!newGen) -> do+                    put newGen+                    return (f $! i)  instance MonadRandom (State StdGen) where-    getRandomByte   = getRandomByteFromRandomGenState-    getRandomWord   = getRandomWordFromRandomGenState-    getRandomDouble = getRandomDoubleFromRandomGenState+    supportedPrims _ _ = True+    getSupportedRandomPrim = getRandomPrimFromRandomGenState  instance Monad m => MonadRandom (StateT StdGen m) where-    getRandomByte   = getRandomByteFromRandomGenState-    getRandomWord   = getRandomWordFromRandomGenState-    getRandomDouble = getRandomDoubleFromRandomGenState+    supportedPrims _ _ = True+    getSupportedRandomPrim = getRandomPrimFromRandomGenState  instance MonadRandom (S.State StdGen) where-    getRandomByte   = getRandomByteFromRandomGenState-    getRandomWord   = getRandomWordFromRandomGenState-    getRandomDouble = getRandomDoubleFromRandomGenState+    supportedPrims _ _ = True+    getSupportedRandomPrim = getRandomPrimFromRandomGenState  instance Monad m => MonadRandom (S.StateT StdGen m) where-    getRandomByte   = getRandomByteFromRandomGenState-    getRandomWord   = getRandomWordFromRandomGenState-    getRandomDouble = getRandomDoubleFromRandomGenState+    supportedPrims _ _ = True+    getSupportedRandomPrim = getRandomPrimFromRandomGenState+