mwc-random 0.11.0.0 → 0.12.0.0

raw patch · 8 files changed

+433/−45 lines, 8 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• README.markdown +0/−23
• System/Random/MWC.hs +71/−14
• System/Random/MWC/CondensedTable.hs +269/−0
• System/Random/MWC/Distributions.hs +2/−2
• benchmarks/Benchmark.hs +31/−1
• mwc-random.cabal +2/−1
• test/visual.R +43/−0
• test/visual.hs +15/−4

@@ -5,29 +5,6 @@ time-efficient way.  -# Performance--This library has been carefully optimised for high performance.  To-obtain the best runtime efficiency, it is imperative to compile-libraries and applications that use this library using a high level of-optimisation.--Suggested GHC options:--    -O -fvia-C -funbox-strict-fields--To illustrate, here are the times (in seconds) to generate and sum 250-million random Word32 values, on a laptop with a 2.4GHz Core2 Duo-P8600 processor, running Fedora 11 and GHC 6.10.3:--    no flags   200+-    -O           1.249-    -O -fvia-C   0.991--As the numbers above suggest, compiling without optimisation will-yield unacceptable performance.-- # Get involved!  Please report bugs via the

@@ -37,19 +37,23 @@ -- The simplest use is to generate a vector of uniformly distributed values: -- -- @---   vs <- withSystemRandom (uniformVector 100)+--   vs <- withSystemRandom . asGenST $ \gen -> uniformVector gen 100 -- @ ----- These values can be of any type which is an instance of the class 'Variate'.+-- These values can be of any type which is an instance of the class+-- 'Variate'. ----- To generate random values on demand, first 'create' a random number generator.+-- To generate random values on demand, first 'create' a random number+-- generator. -- -- @ --   gen <- create -- @ ----- Keep this generator and use it wherever random values are required. Get a random--- value using 'uniform' or 'uniformR':+-- Hold onto this generator and use it wherever random values are+-- required (creating a new generator is expensive compared to+-- generating a random number, so you don't want to throw them+-- away). Get a random value using 'uniform' or 'uniformR': -- -- @ --   v <- uniform gen@@ -62,12 +66,17 @@     (     -- * Gen: Pseudo-Random Number Generators       Gen-    , GenIO-    , GenST     , create     , initialize     , withSystemRandom +    -- ** Type helpers+    -- $typehelp+    , GenIO+    , GenST+    , asGenIO+    , asGenST+     -- * Variates: uniformly distributed values     , Variate(..)     , uniformVector@@ -294,12 +303,20 @@ -- | State of the pseudo-random number generator. newtype Gen s = Gen (M.MVector s Word32) --- | A shorter name for PRNG state in the IO monad.+-- | A shorter name for PRNG state in the 'IO' monad. type GenIO = Gen (PrimState IO) --- | A shorter name for PRNG state in the ST monad.+-- | A shorter name for PRNG state in the 'ST' monad. type GenST s = Gen (PrimState (ST s)) +-- | Constrain the type of an action to run in the 'IO' monad.+asGenIO :: (GenIO -> IO a) -> (GenIO -> IO a)+asGenIO = id++-- | Constrain the type of an action to run in the 'ST' monad.+asGenST :: (GenST s -> ST s a) -> (GenST s -> ST s a)+asGenST = id+ ioff, coff :: Int ioff = 256 coff = 257@@ -389,7 +406,7 @@   let n    = fromIntegral (numerator t) :: Word64   return [fromIntegral c, fromIntegral n, fromIntegral (n `shiftR` 32)] --- Aquire seed from /dev/urandom+-- | Acquire seed from /dev/urandom acquireSeedSystem :: IO [Word32] acquireSeedSystem = do   let nbytes = 1024@@ -400,12 +417,12 @@     peekArray (nread `div` 4) buf    -- | Seed a PRNG with data from the system's fast source of--- pseudo-random numbers (\"\/dev\/urandom\" on Unix-like systems),+-- pseudo-random numbers (\"@\/dev\/urandom@\" on Unix-like systems), -- then run the given action. ----- This is a heavyweight function, intended to be called only--- occasionally (e.g. once per thread).  You should use the `Gen` it--- creates to generate many random numbers.+-- This is a somewhat expensive function, and is intended to be called+-- only occasionally (e.g. once per thread).  You should use the `Gen`+-- it creates to generate many random numbers. -- -- /Note/: on Windows, this code does not yet use the native -- Cryptographic API as a source of random numbers (it uses the system@@ -601,3 +618,43 @@ --   (2007). Gaussian random number generators. --   /ACM Computing Surveys/ 39(4). --   <http://www.cse.cuhk.edu.hk/~phwl/mt/public/archives/papers/grng_acmcs07.pdf>++-- $typehelp+--+-- The functions in this package are deliberately written for+-- flexibility, and will run in both the 'IO' and 'ST' monads.+--+-- This can defeat the compiler's ability to infer a principal type in+-- simple (and common) cases.  For instance, we would like the+-- following to work cleanly:+--+-- > import System.Random.MWC+-- > import Data.Vector.Unboxed+-- >+-- > main = do+-- >   v <- withSystemRandom $ \gen -> uniformVector gen 20+-- >   print (v :: Vector Int)+--+-- Unfortunately, the compiler cannot tell what monad 'uniformVector'+-- should execute in.  The \"fix\" of adding explicit type annotations+-- is not pretty:+--+-- > {-# LANGUAGE ScopedTypeVariables #-}+-- >+-- > import Control.Monad.ST+-- >+-- > main = do+-- >   vs <- withSystemRandom $+-- >         \(gen::GenST s) -> uniformVector gen 20 :: ST s (Vector Int)+-- >   print vs+--+-- As a more readable alternative, this library provides 'asGenST' and+-- 'asGenIO' to constrain the types appropriately.  We can get rid of+-- the explicit type annotations as follows:+--+-- > main = do+-- >   vs <- withSystemRandom . asGenST $ \gen -> uniformVector gen 20+-- >   print (vs :: Vector Int)+--+-- This is almost as compact as the original code that the compiler+-- rejected.

@@ -0,0 +1,269 @@+{-# LANGUAGE FlexibleContexts #-}+-- |+-- Module    : System.Random.MWC.CondensedTable+-- Copyright : (c) 2012 Aleksey Khudyakov+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Table-driven generation of random variates.  This approach can+-- generate random variates in /O(1)/ time for the supported+-- distributions, at a modest cost in initialization time.+module System.Random.MWC.CondensedTable (+    -- * Condensed tables+    CondensedTable+  , CondensedTableV+  , CondensedTableU+  , genFromTable+    -- * Constructors for tables+  , tableFromProbabilities+  , tableFromWeights+  , tableFromIntWeights+    -- ** Disrete distributions+  , tablePoisson+  , tableBinomial+  ) where++import Control.Arrow           (second,(***))+import Control.Monad.Primitive (PrimMonad(..))++import Data.Word+import Data.Int+import Data.Bits+import qualified Data.Vector.Generic         as G+import           Data.Vector.Generic           ((++))+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Unboxed         as U+import qualified Data.Vector                 as V+import Data.Vector.Generic (Vector)++import Prelude hiding ((++))++import System.Random.MWC++++-- | A lookup table for arbitrary discrete distributions. It allows+-- the generation of random variates in /O(1)/. Note that probability+-- is quantized in units of @1/2^32@, and all distributions with+-- infinite support (e.g. Poisson) should be truncated.+data CondensedTable v a =+  CondensedTable+  {-# UNPACK #-} !Word64 !(v a) -- Lookup limit and first table+  {-# UNPACK #-} !Word64 !(v a) -- Second table+  {-# UNPACK #-} !Word64 !(v a) -- Third table+  !(v a)                        -- Last table++-- Implementation note. We have to store lookup limit in Word64 since+-- we need to accomodate two cases. First is when we have no values in+-- lookup table, second is when all elements are there+--+-- Both are pretty easy to realize. For first one probability of every+-- outcome should be less then 1/256, latter arise when probabilities+-- of two outcomes are [0.5,0.5]++-- | A 'CondensedTable' that uses unboxed vectors.+type CondensedTableU = CondensedTable U.Vector++-- | A 'CondensedTable' that uses boxed vectors, and is able to hold+-- any type of element.+type CondensedTableV = CondensedTable V.Vector++++-- | Generate a random value using a condensed table.+genFromTable :: (PrimMonad m, Vector v a) =>+                CondensedTable v a -> Gen (PrimState m) -> m a+{-# INLINE genFromTable #-}+genFromTable table gen = do+  w <- uniform gen+  return $ lookupTable table $ fromIntegral (w :: Word32)++lookupTable :: Vector v a => CondensedTable v a -> Word64 -> a+{-# INLINE lookupTable #-}+lookupTable (CondensedTable na aa nb bb nc cc dd) i+  | i < na    = aa `at` ( i       `shiftR` 24)+  | i < nb    = bb `at` ((i - na) `shiftR` 16)+  | i < nc    = cc `at` ((i - nb) `shiftR` 8 )+  | otherwise = dd `at` ( i - nc)+  where+    at arr j = G.unsafeIndex arr (fromIntegral j)+++----------------------------------------------------------------+-- Table generation+----------------------------------------------------------------++-- | Generate a condensed lookup table from a list of outcomes with+-- given probabilities. The vector should be non-empty and the+-- probabilites should be non-negative and sum to 1. If this is not+-- the case, this algorithm will construct a table for some+-- distribution that may bear no resemblance to what you intended.+tableFromProbabilities+    :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32)+       => v (a, Double) -> CondensedTable v a+{-# INLINE tableFromProbabilities #-}+tableFromProbabilities v+  | G.null v  = pkgError "tableFromProbabilities" "empty vector of outcomes"+  | otherwise = tableFromIntWeights $ G.map (second $ round . (* mlt)) v+  where+    mlt = 4.294967296e9 -- 2^32++-- | Same as 'tableFromProbabilities' but treats number as weights not+-- probilities. Non-positive weights are discarded, and those+-- remaining are normalized to 1.+tableFromWeights+    :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32)+       => v (a, Double) -> CondensedTable v a+{-# INLINE tableFromWeights #-}+tableFromWeights = tableFromProbabilities . normalize . G.filter ((> 0) . snd)+  where+    normalize v+      | G.null v  = pkgError "tableFromWeights" "no positive weights"+      | otherwise = G.map (second (/ s)) v+      where+        -- Explicit fold is to avoid 'Vector v Double' constraint+        s = G.foldl' (flip $ (+) . snd) 0 v+++-- | Generate a condensed lookup table from integer weights. Weights+-- should sum to @2^32@. If they don't, the algorithm will alter the+-- weights so that they do. This approach should work reasonably well+-- for rounding error.+tableFromIntWeights :: (Vector v (a,Word32), Vector v a, Vector v Word32)+                    => v (a, Word32)+                    -> CondensedTable v a+{-# INLINE tableFromIntWeights #-}+tableFromIntWeights tbl+  | n == 0    = pkgError "tableFromIntWeights" "empty table"+    -- Single element tables should be treated sepately. Otherwise+    -- they will confuse correctWeights+  | n == 1    = let m = 2^(32::Int) - 1 -- Works for both Word32 & Word64+                in CondensedTable+                   m (G.replicate 256 $ fst $ G.head tbl)+                   m  G.empty+                   m  G.empty+                      G.empty+  | otherwise = CondensedTable+                na aa+                nb bb+                nc cc+                   dd+  where+    n     = G.length tbl+    -- Corrected table+    table = uncurry G.zip $ id *** correctWeights $ G.unzip tbl+    -- Make condensed table+    mkTable  d =+      G.concatMap (\(x,w) -> G.replicate (fromIntegral $ digit d w) x) table+    len = fromIntegral . G.length+    -- Tables+    aa = mkTable 0+    bb = mkTable 1+    cc = mkTable 2+    dd = mkTable 3+    -- Offsets+    na =       len aa `shiftL` 24+    nb = na + (len bb `shiftL` 16)+    nc = nb + (len cc `shiftL` 8)+++-- Calculate N'th digit base 256+digit :: Int -> Word32 -> Word32+digit 0 x =  x `shiftR` 24+digit 1 x = (x `shiftR` 16) .&. 0xff+digit 2 x = (x `shiftR` 8 ) .&. 0xff+digit 3 x =  x .&. 0xff+digit _ _ = pkgError "digit" "the impossible happened!?"+{-# INLINE digit #-}++-- Correct integer weights so they sum up to 2^32. Array of weight+-- should contain at least 2 elements.+correctWeights :: G.Vector v Word32 => v Word32 -> v Word32+{-# INLINE correctWeights #-}+correctWeights v = G.create $ do+  let+    -- Sum of weights+    s = G.foldl' (flip $ (+) . fromIntegral) 0 v :: Int64+    -- Array size+    n = G.length v+  arr <- G.thaw v+  -- On first pass over array adjust only entries which are larger+  -- than `lim'. On second and subsequent passes `lim' is set to 1.+  --+  -- It's possibly to make this algorithm loop endlessly if all+  -- weights are 1 or 0.+  let loop lim i delta+        | delta == 0 = return ()+        | i >= n     = loop 1 0 delta+        | otherwise  = do+            w <- M.read arr i+            case () of+              _| w < lim   -> loop lim (i+1) delta+               | delta < 0 -> M.write arr i (w + 1) >> loop lim (i+1) (delta + 1)+               | otherwise -> M.write arr i (w - 1) >> loop lim (i+1) (delta - 1)+  loop 255 0 (s - 2^(32::Int))+  return arr+++-- | Create a lookup table for the Poisson distibution. Note that+-- table construction may have significant cost. For &#955; < 100 it+-- takes as much time to build table as generation of 1000-30000+-- variates.+tablePoisson :: Double -> CondensedTableU Int+tablePoisson = tableFromProbabilities . make+  where+    make lam+      | lam < 0    = pkgError "tablePoisson" "negative lambda"+      | lam < 22.8 = U.unfoldr unfoldForward (exp (-lam), 0)+      | otherwise  = U.unfoldr unfoldForward (pMax, nMax)+                  ++ U.tail (U.unfoldr unfoldBackward (pMax, nMax))+      where+        -- Number with highest probability and its probability+        nMax = floor lam :: Int+        pMax = let c = lam * exp( -lam / fromIntegral nMax )+               in  U.foldl' (\p i -> p * c / i) 1 (U.enumFromN 1 nMax)+        -- Build probability list+        unfoldForward (p,i)+          | p < minP  = Nothing+          | otherwise = Just ( (i,p)+                             , (p * lam / fromIntegral (i+1), i+1)+                             )+        -- Go down+        unfoldBackward (p,i)+          | p < minP  = Nothing+          | otherwise = Just ( (i,p)+                             , (p / lam * fromIntegral i, i-1)+                             )+    minP = 1.1641532182693481e-10 -- 2**(-33)++-- | Create a lookup table for the binomial distribution.+tableBinomial :: Int            -- ^ Number of tries+              -> Double         -- ^ Probability of success+              -> CondensedTableU Int+tableBinomial n p = tableFromProbabilities makeBinom+  where +  makeBinom+    | n <= 0         = pkgError "tableBinomial" "non-positive number of tries"+    | p == 0         = U.singleton (0,1)+    | p == 1         = U.singleton (n,1)+    | p > 0 && p < 1 = U.unfoldrN (n + 1) unfolder ((1-p)^n, 0)+    | otherwise      = pkgError "tableBinomial" "probability is out of range"+    where+      h = p / (1 - p)+      unfolder (t,i) = Just ( (i,t)+                            , (t * (fromIntegral $ n + 1 - i1) * h / fromIntegral i1, i1) )+        where i1 = i + 1++pkgError :: String -> String -> a+pkgError func err =+    error . concat $ ["System.Random.MWC.CondensedTable.", func, ": ", err]++-- $references+--+-- * Wang, J.; Tsang, W. W.; G. Marsaglia (2004), Fast Generation of+--   Discrete Random Variables, /Journal of Statistical Software,+--   American Statistical Association/, vol. 11(i03).+--   <http://ideas.repec.org/a/jss/jstsof/11i03.html>

@@ -95,7 +95,7 @@                 else return $! if neg then x - r else r - x  --- | Generate exponentially distributed random variate.+-- | Generate an exponentially distributed random variate. exponential :: PrimMonad m             => Double            -- ^ Scale parameter             -> Gen (PrimState m) -- ^ Generator@@ -139,7 +139,7 @@       a2 = 1 / sqrt(9 * a1)  --- | Random variate generator for chi square distribution.+-- | Random variate generator for the chi square distribution. chiSquare :: PrimMonad m           => Int                -- ^ Number of degrees of freedom           -> Gen (PrimState m)  -- ^ Generator

@@ -3,15 +3,23 @@ import Criterion.Main import Data.Int import Data.Word+import qualified Data.Vector.Unboxed as U import qualified System.Random as R import System.Random.MWC import System.Random.MWC.Distributions+import System.Random.MWC.CondensedTable import qualified System.Random.Mersenne as M +makeTableUniform :: Int -> CondensedTable U.Vector Int+makeTableUniform n =+  tableFromProbabilities $ U.zip (U.enumFromN 0 n) (U.replicate n (1 / fromIntegral n))+{-# INLINE makeTableUniform #-}++ main = do   mwc <- create   mtg <- M.newMTGen . Just =<< uniform mwc-  defaultMain +  defaultMain     [ bgroup "mwc"       -- One letter group names are used so they will fit on the plot.       --@@ -53,6 +61,28 @@         , bench "gamma,a<1"   (gamma 0.5 1   mwc :: IO Double)         , bench "gamma,a>1"   (gamma 2   1   mwc :: IO Double)         , bench "chiSquare"   (chiSquare 4   mwc :: IO Double)+        ]+      , bgroup "CT/gen" $ concat+        [ [ bench ("uniform "++show i)     (genFromTable (makeTableUniform i) mwc :: IO Int)+          | i <- [2..10]+          ]+        , [ bench ("poisson " ++ show l)   (genFromTable (tablePoisson l) mwc :: IO Int)+          | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]+          ]+        , [ bench ("binomial " ++ show p ++ " " ++ show n) (genFromTable (tableBinomial n p) mwc :: IO Int)+          | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]+          ]+        ]+      , bgroup "CT/table" $ concat+        [ [ bench ("uniform " ++ show i) $ whnf makeTableUniform i+          | i <- [2..30]+          ]+        , [ bench ("poisson " ++ show l) $ whnf tablePoisson l+          | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]+          ]+        , [ bench ("binomial " ++ show p ++ " " ++ show n) $ whnf (tableBinomial n) p+          | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]+          ]         ]       ]     , bgroup "random"

@@ -1,5 +1,5 @@ name:           mwc-random-version:        0.11.0.0+version:        0.12.0.0 synopsis:       Fast, high quality pseudo random number generation description:   This package contains code for generating high quality random@@ -37,6 +37,7 @@   exposed-modules:     System.Random.MWC     System.Random.MWC.Distributions+    System.Random.MWC.CondensedTable   build-depends:     base < 5,     primitive,

@@ -3,13 +3,28 @@   view.dumps <- function() {+  # Load random data from dist   load.d <- function(name) read.table(name)[,1]+  # Plots for continous distribution   plot.d <- function(name, dens, rng) {     smp <- load.d( name )     plot( density(smp), xlim=rng, main=name, col='blue', lwd=2)     hist( smp, probability=TRUE, breaks=100, add=TRUE)     plot( dens, xlim=rng, col='red', add=TRUE, lwd=2)   }+  # plots for discrete distribution+  plot.ds <- function( name, xs, prob) {+    smp <- load.d( name )+    h   <- hist( smp,+                 breaks = c( max(xs) + 0.5, xs - 0.5),+                 freq=FALSE, main = name+                )+    dh  <- sqrt( h$count ) / max( 1, sum( h$count ) )+    arrows( xs, h$density + dh,+            xs, h$density - dh,+            angle=90, code=3, length=0.2 )+    points( xs, prob(xs), pch='0', col='red', type='b')+  }   ################################################################   # Normal   plot.d ("distr/normal-0-1",@@ -59,4 +74,32 @@           function(x) dexp(x,3),           c(-0.5, 3) )   readline()+  ################################################################+  # Poisson+  plot.ds( "distr/poisson-0.1", 0:6, function(x) dpois(x, lambda=0.1) )+  readline()+  #+  plot.ds( "distr/poisson-1.0", 0:10, function(x) dpois(x, lambda=1.0) )+  readline()+  #+  plot.ds( "distr/poisson-4.5", 0:20, function(x) dpois(x, lambda=4.5) )+  readline()+  #+  plot.ds( "distr/poisson-30", 0:100, function(x) dpois(x, lambda=30) )+  readline()+  #+  ################################################################+  # Binomial+  plot.ds( "distr/binom-4-0.5", 0:4, function(x) dbinom(x, 4, 0.5) )+  readline()+  #+  plot.ds( "distr/binom-10-0.1", 0:10, function(x) dbinom(x, 10, 0.1) )+  readline()+  #+  plot.ds( "distr/binom-10-0.6", 0:10, function(x) dbinom(x, 10, 0.6) )+  readline()+  #+  plot.ds( "distr/binom-10-0.8", 0:10, function(x) dbinom(x, 10, 0.8) )+  readline()+  # }

@@ -4,18 +4,19 @@ import System.Directory  (createDirectoryIfMissing,setCurrentDirectory) import System.IO -import qualified System.Random.MWC               as MWC-import qualified System.Random.MWC.Distributions as MWC+import qualified System.Random.MWC                as MWC+import qualified System.Random.MWC.Distributions  as MWC+import qualified System.Random.MWC.CondensedTable as MWC   dumpSample :: Show a => Int -> FilePath -> IO a -> IO () dumpSample n fname gen =   withFile fname WriteMode $ \h ->      replicateM_ n (hPutStrLn h . show =<< gen)-  + main :: IO () main = MWC.withSystemRandom $ \g -> do-  let n   = 10000+  let n   = 30000       dir = "distr"   createDirectoryIfMissing True dir   setCurrentDirectory           dir@@ -31,3 +32,13 @@   -- Exponential   dumpSample n "exponential-1" $ MWC.exponential 1 g   dumpSample n "exponential-3" $ MWC.exponential 3 g+  -- Poisson+  dumpSample n "poisson-0.1"   $ MWC.genFromTable (MWC.tablePoisson 0.1) g+  dumpSample n "poisson-1.0"   $ MWC.genFromTable (MWC.tablePoisson 1.0) g+  dumpSample n "poisson-4.5"   $ MWC.genFromTable (MWC.tablePoisson 4.5) g+  dumpSample n "poisson-30"    $ MWC.genFromTable (MWC.tablePoisson 30)  g+  -- Binomial+  dumpSample n "binom-4-0.5"   $ MWC.genFromTable (MWC.tableBinomial 4  0.5) g+  dumpSample n "binom-10-0.1"  $ MWC.genFromTable (MWC.tableBinomial 10 0.1) g  +  dumpSample n "binom-10-0.6"  $ MWC.genFromTable (MWC.tableBinomial 10 0.6) g+  dumpSample n "binom-10-0.8"  $ MWC.genFromTable (MWC.tableBinomial 10 0.8) g