diff --git a/README.markdown b/README.markdown --- a/README.markdown +++ b/README.markdown @@ -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 diff --git a/System/Random/MWC.hs b/System/Random/MWC.hs --- a/System/Random/MWC.hs +++ b/System/Random/MWC.hs @@ -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). -- + +-- $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. diff --git a/System/Random/MWC/CondensedTable.hs b/System/Random/MWC/CondensedTable.hs new file mode 100644 --- /dev/null +++ b/System/Random/MWC/CondensedTable.hs @@ -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 λ < 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). +-- diff --git a/System/Random/MWC/Distributions.hs b/System/Random/MWC/Distributions.hs --- a/System/Random/MWC/Distributions.hs +++ b/System/Random/MWC/Distributions.hs @@ -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 diff --git a/benchmarks/Benchmark.hs b/benchmarks/Benchmark.hs --- a/benchmarks/Benchmark.hs +++ b/benchmarks/Benchmark.hs @@ -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" diff --git a/mwc-random.cabal b/mwc-random.cabal --- a/mwc-random.cabal +++ b/mwc-random.cabal @@ -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, diff --git a/test/visual.R b/test/visual.R --- a/test/visual.R +++ b/test/visual.R @@ -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() + # } diff --git a/test/visual.hs b/test/visual.hs --- a/test/visual.hs +++ b/test/visual.hs @@ -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