diff --git a/System/Random/MWC.hs b/System/Random/MWC.hs --- a/System/Random/MWC.hs +++ b/System/Random/MWC.hs @@ -2,7 +2,7 @@ MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies, UnboxedTuples #-} -- | -- Module : System.Random.MWC --- Copyright : (c) 2009, 2010, 2011 Bryan O'Sullivan +-- Copyright : (c) 2009-2012 Bryan O'Sullivan -- License : BSD3 -- -- Maintainer : bos@serpentine.com @@ -10,34 +10,75 @@ -- Portability : portable -- -- Pseudo-random number generation. This module contains code for --- generating high quality random numbers that follow either a uniform --- or normal distribution. +-- generating high quality random numbers that follow a uniform +-- distribution. -- +-- For non-uniform distributions, see the +-- 'System.Random.MWC.Distributions' module. +-- -- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222) -- multiply-with-carry generator, which has a period of 2^8222 and -- fares well in tests of randomness. It is also extremely fast, -- between 2 and 3 times faster than the Mersenne Twister. +-- +-- The generator state is stored in the 'Gen' data type. It can be +-- created in several ways: +-- +-- 1. Using the 'withSystemRandom' call, which creates a random state. +-- +-- 2. Supply your own seed to 'initialize' function. +-- +-- 3. Finally, 'create' makes a generator from a fixed seed. +-- Generators created in this way aren't really random. +-- +-- For repeatability, the state of the generator can be snapshotted +-- and replayed using the 'save' and 'restore' functions. +-- +-- The simplest use is to generate a vector of uniformly distributed values: +-- +-- @ +-- vs <- withSystemRandom (uniformVector 100) +-- @ +-- +-- 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. +-- +-- @ +-- gen <- create +-- @ +-- +-- Keep this generator and use it wherever random values are required. Get a random +-- value using 'uniform' or 'uniformR': +-- +-- @ +-- v <- uniform gen +-- @ +-- +-- @ +-- v <- uniformR (1, 52) gen +-- @ module System.Random.MWC ( - -- * Types + -- * Gen: Pseudo-Random Number Generators Gen , GenIO , GenST - , Seed - , fromSeed - , toSeed - , Variate(..) - -- * Other distributions - , normal - -- * Creation , create , initialize , withSystemRandom - -- * State management + + -- * Variates: uniformly distributed values + , Variate(..) + , uniformVector + + -- * Seed: state management + , Seed + , fromSeed + , toSeed , save , restore - -- * Helper functions - , uniformVector + -- * References -- $references ) where @@ -56,24 +97,20 @@ import Data.Ratio ((%), numerator) import Data.Time.Clock.POSIX (getPOSIXTime) import Data.Typeable (Typeable) -import Data.Vector.Generic (Vector, unsafeFreeze) +import Data.Vector.Generic (Vector) import Data.Word (Word, Word8, Word16, Word32, Word64) import Foreign.Marshal.Alloc (allocaBytes) import Foreign.Marshal.Array (peekArray) import Prelude hiding (catch) import qualified Data.Vector.Generic as G -import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Unboxed as I import qualified Data.Vector.Unboxed.Mutable as M import System.CPUTime (cpuTimePrecision, getCPUTime) import System.IO (IOMode(..), hGetBuf, hPutStrLn, stderr, withBinaryFile) import System.IO.Unsafe (unsafePerformIO) --- FIXME: removal of Unbox constraint leads to severe (~10x) --- performance drop with GHC 6.12. For details see bug #33 in the --- vector bug tracker[1] --- [1] http://trac.haskell.org/vector/ticket/33 + -- | The class of types for which we can generate uniformly -- distributed random variates. -- @@ -84,7 +121,7 @@ -- -- /Note/: Marsaglia's PRNG is not known to be cryptographically -- secure, so you should not use it for cryptographic operations. -class M.Unbox a => Variate a where +class Variate a where -- | Generate a single uniformly distributed random variate. The -- range of values produced varies by type: -- @@ -111,49 +148,49 @@ instance Variate Int8 where uniform = uniform1 fromIntegral - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Int16 where uniform = uniform1 fromIntegral - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Int32 where uniform = uniform1 fromIntegral - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Int64 where uniform = uniform2 wordsTo64Bit - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Word8 where uniform = uniform1 fromIntegral - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Word16 where uniform = uniform1 fromIntegral - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Word32 where uniform = uniform1 fromIntegral - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Word64 where uniform = uniform2 wordsTo64Bit - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} @@ -184,7 +221,7 @@ #else uniform = uniform2 wordsTo64Bit #endif - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} @@ -194,7 +231,7 @@ #else uniform = uniform2 wordsTo64Bit #endif - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} @@ -287,10 +324,12 @@ -- verbatim, then its elements are 'xor'ed against elements of the -- default seed until 256 elements are reached. -- --- If a seed contains exactly 258 elements last two elements are used --- to set generator state. It's to ensure that @gen' == gen@ +-- If a seed contains exactly 258 elements, then the last two elements +-- are used to set the generator's initial state. This allows for +-- complete generator reproducibility, so that e.g. @gen' == gen@ in +-- the following example: -- --- > gen' <- initialize . fromSeed =<< save +-- @gen' <- 'initialize' . 'fromSeed' =<< 'save'@ initialize :: (PrimMonad m, Vector v Word32) => v Word32 -> m (Gen (PrimState m)) initialize seed = do @@ -330,29 +369,14 @@ toSeed :: (Vector v Word32) => v Word32 -> Seed toSeed v = Seed $ I.create $ do { Gen q <- initialize v; return q } --- Safe version of unsafeFreeze. --- NOTE: vector-0.7 will provide function `freeze' with same --- functionality. This function shall be removed when support for --- vector<=0.6 is dropped -safeFreeze :: (PrimMonad m, Vector v a) => G.Mutable v (PrimState m) a -> m (v a) -safeFreeze v = do - v' <- GM.unsafeNew (GM.length v) - GM.unsafeCopy v' v - unsafeFreeze v' - -- | Save the state of a 'Gen', for later use by 'restore'. save :: PrimMonad m => Gen (PrimState m) -> m Seed -save (Gen q) = Seed `liftM` safeFreeze q +save (Gen q) = Seed `liftM` G.freeze q {-# INLINE save #-} --- NOTE: with vector-0.7 all code could be replaced with `clone' -- | Create a new 'Gen' that mirrors the state of a saved 'Seed'. restore :: PrimMonad m => Seed -> m (Gen (PrimState m)) -restore (Seed s) = M.unsafeNew n >>= fill - where fill q = go 0 where - go !i | i >= n = return $! Gen q - | otherwise = M.unsafeWrite q i (I.unsafeIndex s i) >> go (i+1) - n = I.length s +restore (Seed s) = Gen `liftM` G.thaw s {-# INLINE restore #-} @@ -379,6 +403,10 @@ -- 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. +-- -- /Note/: on Windows, this code does not yet use the native -- Cryptographic API as a source of random numbers (it uses the system -- clock instead). As a result, the sequences it generates may not be @@ -404,16 +432,16 @@ where j = fromIntegral (i+1) :: Word8 {-# INLINE nextIndex #-} -a :: Word64 -a = 1540315826 -{-# INLINE a #-} +aa :: Word64 +aa = 1540315826 +{-# INLINE aa #-} uniformWord32 :: PrimMonad m => Gen (PrimState m) -> m Word32 uniformWord32 (Gen q) = do i <- nextIndex `liftM` M.unsafeRead q ioff c <- fromIntegral `liftM` M.unsafeRead q coff qi <- fromIntegral `liftM` M.unsafeRead q i - let t = a * qi + c + let t = aa * qi + c c' = fromIntegral (t `shiftR` 32) x = fromIntegral t + c' (# x', c'' #) | x < c' = (# x + 1, c' + 1 #) @@ -437,12 +465,12 @@ c <- fromIntegral `liftM` M.unsafeRead q coff qi <- fromIntegral `liftM` M.unsafeRead q i qj <- fromIntegral `liftM` M.unsafeRead q j - let t = a * qi + c + let t = aa * qi + c c' = fromIntegral (t `shiftR` 32) x = fromIntegral t + c' (# x', c'' #) | x < c' = (# x + 1, c' + 1 #) | otherwise = (# x, c' #) - u = a * qj + fromIntegral c'' + u = aa * qj + fromIntegral c'' d' = fromIntegral (u `shiftR` 32) y = fromIntegral u + d' (# y', d'' #) | y < d' = (# y + 1, d' + 1 #) @@ -475,31 +503,36 @@ -- unsigned data type of same size sub :: (Integral a, Integral (Unsigned a)) => a -> a -> Unsigned a sub x y = fromIntegral x - fromIntegral y +{-# INLINE sub #-} add :: (Integral a, Integral (Unsigned a)) => a -> Unsigned a -> a add m x = m + fromIntegral x - --- Generate uniform value in the range [0,n). Values must be --- unsigned. Second parameter is random number generator -unsignedRange :: (PrimMonad m, Integral a, Bounded a) => a -> m a -> m a -unsignedRange n rnd = go - where - buckets = maxBound `div` n - maxN = buckets * n - go = do x <- rnd - if x < maxN then return (x `div` buckets) - else go -{-# INLINE unsignedRange #-} +{-# INLINE add #-} --- Generate unformly distributed value in inclusive range. +-- Generate uniformly distributed value in inclusive range. +-- +-- NOTE: This function must be fully applied. Otherwise it won't be +-- inlined, which will cause a severe performance loss. +-- +-- > uniformR = uniformRange -- won't be inlined +-- > uniformR a b = uniformRange a b -- will be inlined uniformRange :: ( PrimMonad m , Integral a, Bounded a, Variate a , Integral (Unsigned a), Bounded (Unsigned a), Variate (Unsigned a)) => (a,a) -> Gen (PrimState m) -> m a uniformRange (x1,x2) g - | x1 == minBound && x2 == maxBound = uniform g - | otherwise = do x <- unsignedRange (sub x2 x1 + 1) (uniform g) - return $! add x1 x + | n == 0 = uniform g -- Abuse overflow in unsigned types + | otherwise = loop + where + -- Allow ranges where x2 Gen (PrimState m) -> m Double -normal gen = loop - where - loop = do - u <- (subtract 1 . (*2)) `liftM` uniform gen - ri <- uniform gen - let i = fromIntegral ((ri :: Word32) .&. 127) - bi = I.unsafeIndex blocks i - bj = I.unsafeIndex blocks (i+1) - if abs u < I.unsafeIndex ratios i - then return $! u * bi - else if i == 0 - then normalTail (u < 0) - else do - let x = u * bi - xx = x * x - d = exp (-0.5 * (bi * bi - xx)) - e = exp (-0.5 * (bj * bj - xx)) - c <- uniform gen - if e + c * (d - e) < 1 - then return x - else loop - blocks = let f = exp (-0.5 * r * r) - in (`I.snoc` 0) . I.cons (v/f) . I.cons r . - I.unfoldrN 126 go $! T r f - where - go (T b g) = let !u = T h (exp (-0.5 * h * h)) - h = sqrt (-2 * log (v / b + g)) - in Just (h, u) - v = 9.91256303526217e-3 - {-# NOINLINE blocks #-} - r = 3.442619855899 - ratios = I.zipWith (/) (I.tail blocks) blocks - {-# NOINLINE ratios #-} - normalTail neg = tailing - where tailing = do - x <- ((/r) . log) `liftM` uniform gen - y <- log `liftM` uniform gen - if y * (-2) < x * x - then tailing - else return $! if neg then x - r else r - x -{-# INLINE normal #-} - defaultSeed :: I.Vector Word32 defaultSeed = I.fromList [ 0x7042e8b3, 0x06f7f4c5, 0x789ea382, 0x6fb15ad8, 0x54f7a879, 0x0474b184, @@ -609,15 +592,6 @@ {-# NOINLINE defaultSeed #-} -- $references --- --- * Doornik, J.A. (2005) An improved ziggurat method to generate --- normal random samples. Mimeo, Nuffield College, University of --- Oxford. --- --- * Doornik, J.A. (2007) Conversion of high-period random numbers to --- floating point. --- /ACM Transactions on Modeling and Computer Simulation/ 17(1). --- -- -- * Marsaglia, G. (2003) Seeds for random number generators. -- /Communications of the ACM/ 46(5):90–93. diff --git a/System/Random/MWC/Distributions.hs b/System/Random/MWC/Distributions.hs new file mode 100644 --- /dev/null +++ b/System/Random/MWC/Distributions.hs @@ -0,0 +1,171 @@ +{-# LANGUAGE BangPatterns #-} +-- | +-- Module : System.Random.MWC.Distributions +-- Copyright : (c) 2012 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Pseudo-random number generation for non-uniform distributions. + +module System.Random.MWC.Distributions + ( + -- * Variates: non-uniformly distributed values + normal + , standard + , exponential + , gamma + , chiSquare + + -- * References + -- $references + ) where + +import Control.Monad (liftM) +import Control.Monad.Primitive (PrimMonad, PrimState) +import Data.Bits ((.&.)) +import Data.Word (Word32) +import System.Random.MWC (Gen, uniform) +import qualified Data.Vector.Unboxed as I + +-- Unboxed 2-tuple +data T = T {-# UNPACK #-} !Double {-# UNPACK #-} !Double + + +-- | Generate a normally distributed random variate with given mean +-- and standard deviation. +normal :: PrimMonad m + => Double -- ^ Mean + -> Double -- ^ Standard deviation + -> Gen (PrimState m) + -> m Double +{-# INLINE normal #-} +normal m s gen = do + x <- standard gen + return $! m + s * x + +-- | Generate a normally distributed random variate with zero mean and +-- unit variance. +-- +-- The implementation uses Doornik's modified ziggurat algorithm. +-- Compared to the ziggurat algorithm usually used, this is slower, +-- but generates more independent variates that pass stringent tests +-- of randomness. +standard :: PrimMonad m => Gen (PrimState m) -> m Double +{-# INLINE standard #-} +standard gen = loop + where + loop = do + u <- (subtract 1 . (*2)) `liftM` uniform gen + ri <- uniform gen + let i = fromIntegral ((ri :: Word32) .&. 127) + bi = I.unsafeIndex blocks i + bj = I.unsafeIndex blocks (i+1) + case () of + _| abs u < I.unsafeIndex ratios i -> return $! u * bi + | i == 0 -> normalTail (u < 0) + | otherwise -> do + let x = u * bi + xx = x * x + d = exp (-0.5 * (bi * bi - xx)) + e = exp (-0.5 * (bj * bj - xx)) + c <- uniform gen + if e + c * (d - e) < 1 + then return x + else loop + blocks = (`I.snoc` 0) . I.cons (v/f) . I.cons r . I.unfoldrN 126 go $! T r f + where + go (T b g) = let !u = T h (exp (-0.5 * h * h)) + h = sqrt (-2 * log (v / b + g)) + in Just (h, u) + v = 9.91256303526217e-3 + f = exp (-0.5 * r * r) + {-# NOINLINE blocks #-} + r = 3.442619855899 + ratios = I.zipWith (/) (I.tail blocks) blocks + {-# NOINLINE ratios #-} + normalTail neg = tailing + where tailing = do + x <- ((/r) . log) `liftM` uniform gen + y <- log `liftM` uniform gen + if y * (-2) < x * x + then tailing + else return $! if neg then x - r else r - x + + +-- | Generate exponentially distributed random variate. +exponential :: PrimMonad m + => Double -- ^ Scale parameter + -> Gen (PrimState m) -- ^ Generator + -> m Double +{-# INLINE exponential #-} +exponential beta gen = do + x <- uniform gen + return $! - log x / beta + + +-- | Random variate generator for gamma distribution. +gamma :: PrimMonad m + => Double -- ^ Shape parameter + -> Double -- ^ Scale parameter + -> Gen (PrimState m) -- ^ Generator + -> m Double +{-# INLINE gamma #-} +gamma a b gen + | a <= 0 = pkgError "gamma" "negative alpha parameter" + | otherwise = mainloop + where + mainloop = do + T x v <- innerloop + u <- uniform gen + let cont = u > 1 - 0.331 * sqr (sqr x) + && log u > 0.5 * sqr x + a1 * (1 - v + log v) -- Rarely evaluated + case () of + _| cont -> mainloop + | a >= 1 -> return $! a1 * v * b + | otherwise -> do y <- uniform gen + return $! y ** (1 / a) * a1 * v * b + -- inner loop + innerloop = do + x <- standard gen + case 1 + a2*x of + v | v <= 0 -> innerloop + | otherwise -> return $! T x (v*v*v) + -- constants + a' = if a < 1 then a + 1 else a + a1 = a' - 1/3 + a2 = 1 / sqrt(9 * a1) + + +-- | Random variate generator for chi square distribution. +chiSquare :: PrimMonad m + => Int -- ^ Number of degrees of freedom + -> Gen (PrimState m) -- ^ Generator + -> m Double +{-# INLINE chiSquare #-} +chiSquare n gen + | n <= 0 = pkgError "chiSquare" "number of degrees of freedom must be positive" + | otherwise = do x <- gamma (0.5 * fromIntegral n) 1 gen + return $! 2 * x + + +sqr :: Double -> Double +sqr x = x * x +{-# INLINE sqr #-} + +pkgError :: String -> String -> a +pkgError func msg = error $ "System.Random.MWC.Distributions." ++ func ++ + ": " ++ msg + +-- $references +-- +-- * Doornik, J.A. (2005) An improved ziggurat method to generate +-- normal random samples. Mimeo, Nuffield College, University of +-- Oxford. +-- +-- * Doornik, J.A. (2007) Conversion of high-period random numbers to +-- floating point. +-- /ACM Transactions on Modeling and Computer Simulation/ 17(1). +-- diff --git a/benchmarks/Benchmark.hs b/benchmarks/Benchmark.hs --- a/benchmarks/Benchmark.hs +++ b/benchmarks/Benchmark.hs @@ -5,6 +5,7 @@ import Data.Word import qualified System.Random as R import System.Random.MWC +import System.Random.MWC.Distributions import qualified System.Random.Mersenne as M main = do @@ -12,19 +13,47 @@ mtg <- M.newMTGen . Just =<< uniform mwc defaultMain [ bgroup "mwc" - [ bench "Double" (uniform mwc :: IO Double) - , bench "Int" (uniform mwc :: IO Int) - , bench "Int8" (uniform mwc :: IO Int8) - , bench "Int16" (uniform mwc :: IO Int16) - , bench "Int32" (uniform mwc :: IO Int32) - , bench "Int64" (uniform mwc :: IO Int64) - , bench "Word" (uniform mwc :: IO Word) - , bench "Word8" (uniform mwc :: IO Word8) - , bench "Word16" (uniform mwc :: IO Word16) - , bench "Word32" (uniform mwc :: IO Word32) - , bench "Word64" (uniform mwc :: IO Word64) - , bench "Integer" (uniform mwc :: IO Word64) - , bench "normal" (normal mwc :: IO Double) + -- One letter group names are used so they will fit on the plot. + -- + -- U - uniform + -- R - uniformR + -- D - distribution + [ bgroup "U" + [ bench "Double" (uniform mwc :: IO Double) + , bench "Int" (uniform mwc :: IO Int) + , bench "Int8" (uniform mwc :: IO Int8) + , bench "Int16" (uniform mwc :: IO Int16) + , bench "Int32" (uniform mwc :: IO Int32) + , bench "Int64" (uniform mwc :: IO Int64) + , bench "Word" (uniform mwc :: IO Word) + , bench "Word8" (uniform mwc :: IO Word8) + , bench "Word16" (uniform mwc :: IO Word16) + , bench "Word32" (uniform mwc :: IO Word32) + , bench "Word64" (uniform mwc :: IO Word64) + ] + , bgroup "R" + -- I'm not entirely convinced that this is right way to test + -- uniformR. /A.Khudyakov/ + [ bench "Double" (uniformR (-3.21,26) mwc :: IO Double) + , bench "Int" (uniformR (-12,679) mwc :: IO Int) + , bench "Int8" (uniformR (-12,4) mwc :: IO Int8) + , bench "Int16" (uniformR (-12,679) mwc :: IO Int16) + , bench "Int32" (uniformR (-12,679) mwc :: IO Int32) + , bench "Int64" (uniformR (-12,679) mwc :: IO Int64) + , bench "Word" (uniformR (34,633) mwc :: IO Word) + , bench "Word8" (uniformR (34,63) mwc :: IO Word8) + , bench "Word16" (uniformR (34,633) mwc :: IO Word16) + , bench "Word32" (uniformR (34,633) mwc :: IO Word32) + , bench "Word64" (uniformR (34,633) mwc :: IO Word64) + ] + , bgroup "D" + [ bench "standard" (standard mwc :: IO Double) + , bench "normal" (normal 1 3 mwc :: IO Double) + , bench "exponential" (exponential 3 mwc :: IO Double) + , 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 "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.10.0.1 +version: 0.11.0.0 synopsis: Fast, high quality pseudo random number generation description: This package contains code for generating high quality random @@ -23,21 +23,25 @@ copyright: 2009, 2010, 2011 Bryan O'Sullivan category: Math, Statistics build-type: Simple -cabal-version: >= 1.6 +cabal-version: >= 1.8 extra-source-files: README.markdown benchmarks/*.hs benchmarks/Quickie.hs benchmarks/mwc-random-benchmarks.cabal + test/*.R + test/*.sh + test/visual.hs library exposed-modules: System.Random.MWC + System.Random.MWC.Distributions build-depends: base < 5, primitive, time, - vector >= 0.6.0.2 + vector >= 0.7 if impl(ghc >= 6.10) build-depends: base >= 4 @@ -48,6 +52,28 @@ ghc-options: -Wall -funbox-strict-fields if impl(ghc >= 6.8) ghc-options: -fwarn-tabs + +test-suite tests + buildable: False + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: tests.hs + other-modules: KS + QC + Uniform + + ghc-options: + -Wall -threaded -rtsopts + + build-depends: + HUnit, + QuickCheck, + base, + mwc-random, + statistics >= 0.10.1.0, + test-framework, + test-framework-hunit, + test-framework-quickcheck2 source-repository head type: git diff --git a/test/run-dieharder-test.sh b/test/run-dieharder-test.sh new file mode 100644 --- /dev/null +++ b/test/run-dieharder-test.sh @@ -0,0 +1,26 @@ +#!/bin/sh +# +# Run dieharder set of tests for PRNG. All command line parameters are +# passed directly to the dieharder. If no parameters are given -a flag +# is passed which runs all available tests. Full list of dieharder +# options is available at dieharder manpage +# +# NOTE: +# Full set of test require a lot of time to complete. From several +# hours to a few days depending on CPU speed and thoroughness +# settings. +# +# dieharder-source.hs is enthropy source for this test. +# +# This test require dieharder to be installed. It is available at: +# http://www.phy.duke.edu/~rgb/General/dieharder.php + +which dieharder > /dev/null || { echo "dieharder is not found. Aborting"; exit 1; } + +ghc -fforce-recomp -O2 diehard-source +( + date + ./diehard-source | \ + if [ $# = 0 ]; then dieharder -a -g 200; else dieharder "$@" -g 200; fi + date +) | tee diehard.log diff --git a/test/visual.R b/test/visual.R new file mode 100644 --- /dev/null +++ b/test/visual.R @@ -0,0 +1,62 @@ +# Ugly script for displaying distributions alogside with theoretical +# distribution. + + +view.dumps <- function() { + load.d <- function(name) read.table(name)[,1] + 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) + } + ################################################################ + # Normal + plot.d ("distr/normal-0-1", + function(x) dnorm( x, 0, 1 ), + c(-4,4) ) + readline() + # + plot.d ("distr/normal-1-2", + function(x) dnorm( x, 1, 2 ), + c(-6,8) ) + readline(); + + ################################################################ + # Gamma + plot.d ("distr/gamma-1.0-1.0", + function(x) dgamma( x, 1, 1 ), + c(-1,8) ) + readline(); + # + plot.d ("distr/gamma-0.3-0.4", + function(x) dgamma( x, 0.3, scale=0.4 ), + c(-0.25,2) ) + readline(); + # + plot.d ("distr/gamma-0.3-3.0", + function(x) dgamma( x, 0.3, scale=3.0 ), + c(-1,5) ) + readline(); + # + plot.d ("distr/gamma-3.0-0.4", + function(x) dgamma( x, 3.0, scale=0.4 ), + c(-1,6) ) + readline(); + # + plot.d ("distr/gamma-3.0-3.0", + function(x) dgamma( x, 3.0, scale=3.0 ), + c(-1,32) ) + readline(); + ################################################################ + # Exponential + plot.d ("distr/exponential-1", + function(x) dexp(x,1), + c(-0.5, 9) ) + readline() + # + plot.d ("distr/exponential-3", + function(x) dexp(x,3), + c(-0.5, 3) ) + readline() +} diff --git a/test/visual.hs b/test/visual.hs new file mode 100644 --- /dev/null +++ b/test/visual.hs @@ -0,0 +1,33 @@ +-- Generates samples of value for display with visual.R +import Control.Monad + +import System.Directory (createDirectoryIfMissing,setCurrentDirectory) +import System.IO + +import qualified System.Random.MWC as MWC +import qualified System.Random.MWC.Distributions 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 + dir = "distr" + createDirectoryIfMissing True dir + setCurrentDirectory dir + -- Normal + dumpSample n "normal-0-1" $ MWC.normal 0 1 g + dumpSample n "normal-1-2" $ MWC.normal 1 2 g + -- Gamma + dumpSample n "gamma-1.0-1.0" $ MWC.gamma 1.0 1.0 g + dumpSample n "gamma-0.3-0.4" $ MWC.gamma 0.3 0.4 g + dumpSample n "gamma-0.3-3.0" $ MWC.gamma 0.3 3.0 g + dumpSample n "gamma-3.0-0.4" $ MWC.gamma 3.0 0.4 g + dumpSample n "gamma-3.0-3.0" $ MWC.gamma 3.0 3.0 g + -- Exponential + dumpSample n "exponential-1" $ MWC.exponential 1 g + dumpSample n "exponential-3" $ MWC.exponential 3 g