diff --git a/README.markdown b/README.markdown deleted file mode 100644 --- a/README.markdown +++ /dev/null @@ -1,50 +0,0 @@ -# Efficient, general purpose pseudo-random number generation - -This package provides the System.Random.MWC module, a Haskell library -for generating high-quality pseudo-random numbers in a space- and -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 -[bitbucket issue tracker](http://bitbucket.org/bos/mwc-random). - -Master [Mercurial repository](http://bitbucket.org/bos/mwc-random): - -* `hg clone http://bitbucket.org/bos/mwc-random` - -There's also a [git mirror](http://github.com/bos/mwc-random): - -* `git clone git://github.com/bos/mwc-random.git` - -(You can create and contribute changes using either Mercurial or git.) - - -# Authors - -This library is written and maintained by Bryan O'Sullivan, -. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,23 @@ +# Efficient, general purpose pseudo-random number generation +[![Build +Status](https://github.com/haskell/mwc-random/workflows/Haskell-CI/badge.svg)](https://github.com/haskell/mwc-random/actions) + +This package provides the System.Random.MWC module, a Haskell library +for generating high-quality pseudo-random numbers in a space- and +time-efficient way. + + +# Get involved! + +Please report bugs via the +[github issue tracker](https://github.com/haskell/mwc-random/issues). + +Master [git repository](https://github.com/haskell/mwc-random): + +* `git clone https://github.com/haskell/mwc-random.git` + + +# Authors + +This library is written and maintained by Bryan O'Sullivan, +. diff --git a/System/Random/MWC.hs b/System/Random/MWC.hs --- a/System/Random/MWC.hs +++ b/System/Random/MWC.hs @@ -1,41 +1,155 @@ -{-# LANGUAGE BangPatterns, CPP, DeriveDataTypeable, FlexibleContexts, - MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies, UnboxedTuples #-} +{-# LANGUAGE BangPatterns, CPP, DataKinds, DeriveDataTypeable, FlexibleContexts, + FlexibleInstances, MultiParamTypeClasses, MagicHash, Rank2Types, + ScopedTypeVariables, TypeFamilies, UnboxedTuples, TypeOperators + #-} -- | -- 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 -- Stability : experimental -- Portability : portable -- --- Pseudo-random number generation. This module contains code for --- generating high quality random numbers that follow either a uniform --- or normal distribution. +-- Pseudo-random number generation using 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. There are two representation of generator: 'Gen' which is +-- generator that uses in-place mutation and 'Seed' which is immutable +-- snapshot of generator's state. -- --- 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. +-- +-- == Initialization +-- +-- Generator could be initialized in several ways. One is to obtain +-- randomness from operating system using 'createSystemRandom', +-- 'createSystemSeed' or 'withSystemRandomST' (All examples assume +-- that @System.Random.Stateful@ is imported) +-- +-- >>> g <- createSystemRandom +-- >>> uniformM g :: IO Int +-- ... +-- +-- >>> withSystemRandomST $ \g -> uniformM g :: IO Int +-- ... +-- +-- Deterministically create generator from given seed using +-- 'initialize' function: +-- +-- >>> import Data.Int +-- >>> import qualified Data.Vector.Unboxed as U +-- >>> import System.Random.Stateful +-- >>> g <- initialize $ U.fromList [1,2,3] +-- >>> uniformRM (1,200) g :: IO Int64 +-- 101 +-- +-- Last way is to create generator with fixed seed which could be +-- useful in testing +-- +-- >>> g <- create +-- >>> uniformM g :: IO Int +-- -8765701622605876598 +-- +-- +-- == Generation of random numbers +-- +-- Recommended way of generating random numbers in simple cases like +-- generating uniformly distributed random number in range or value +-- uniformly distributed in complete type domain is to use +-- 'UniformRange' and 'Uniform' type classes. Note that while small +-- self-contained examples usually require explicit annotations +-- usually result type could be inferred. +-- +-- This example simulates 20 throws of fair 6-sided dice: +-- +-- >>> g <- create +-- >>> replicateM 20 $ uniformRM (1, 6::Integer) g +-- [3,4,3,1,4,6,1,6,1,4,2,2,3,2,4,2,5,1,3,5] +-- +-- For generating full range of possible values one could use +-- 'uniformM'. This example generates 10 random bytes, or equivalently +-- 10 throws of 256-sided dice: +-- +-- >>> g <- create +-- >>> replicateM 10 $ uniformM g :: IO [Word8] +-- [209,138,126,150,165,15,69,203,155,146] +-- +-- There are special functions for generation of @Doubles@ and @Float +-- in unit interval: 'Random.uniformDouble01M', +-- 'Random.uniformDoublePositive01M', 'Random.uniformFloat01M', +-- 'Random.uniformFloatPositive01M': +-- +-- >>> uniformDouble01M =<< create +-- 0.5248103628705498 +-- >>> uniformFloat01M =<< create +-- 0.5248104 +-- +-- For normal distribution and others see modules +-- "System.Random.MWC.Distributions" and +-- "System.Random.MWC.CondensedTable". Note that they could be used +-- with any other generator implementing 'Random.StatefulGen' API +-- +-- There're special cases for generating random vectors and +-- bytestrings. For example in order to generate random 10-byte +-- sequences as unboxed vector or bytestring: +-- +-- >>> g <- create +-- >>> uniformVector g 10 :: IO (U.Vector Word8) +-- [209,138,126,150,165,15,69,203,155,146] +-- +-- >>> import qualified Data.ByteString as BS +-- >>> g <- create +-- >>> BS.unpack <$> uniformByteStringM 10 g +-- [138,242,130,33,209,248,89,134,150,180] +-- +-- Note that 'Random.uniformByteStringM' produces different result +-- from 'uniformVector' since it uses PRNG's output more efficiently. +-- +-- +-- == State handling +-- +-- For repeatability, the state of the generator can be snapshotted +-- and replayed using the 'save' and 'restore' functions. Following +-- example shows how to save and restore generator: +-- +-- >>> g <- create +-- >>> replicateM_ 10 (uniformM g :: IO Word64) +-- >>> s <- save g +-- >>> uniformM g :: IO Word32 +-- 1771812561 +-- >>> uniformM =<< restore s :: IO Word32 +-- 1771812561 module System.Random.MWC ( - -- * Types + -- * Gen: Pseudo-Random Number Generators Gen + , create + , initialize + , createSystemSeed + , createSystemRandom + , withSystemRandomST + -- ** Type helpers + -- $typehelp , GenIO , GenST - , Seed + , asGenIO + , asGenST + + -- * Variates: uniformly distributed values + , Random.Uniform(..) + , Random.UniformRange(..) , Variate(..) - -- * Other distributions - , normal - -- * Creation - , create - , initialize - , withSystemRandom - -- * State management + , uniformVector + + -- * Seed: state management + , Seed + , fromSeed + , toSeed , save , restore - -- * Helper functions - , uniformVector + -- * Deprecated + , withSystemRandom -- * References -- $references ) where @@ -44,37 +158,33 @@ #include "MachDeps.h" #endif -import Control.Exception (IOException, catch) -import Control.Monad (ap, liftM, unless) -import Control.Monad.Primitive (PrimMonad, PrimState, unsafePrimToIO) -import Control.Monad.ST (ST) -import Data.Bits ((.&.), (.|.), xor) +import Control.Monad (unless) +import Control.Monad.Primitive (PrimMonad, PrimBase, PrimState, unsafePrimToIO, stToPrim) +import Control.Monad.ST (ST,runST) +import Data.Bits ((.&.), (.|.), shiftL, shiftR, xor) import Data.Int (Int8, Int16, Int32, Int64) -import Data.IORef (atomicModifyIORef, newIORef) -import Data.Ratio ((%), numerator) -import Data.Time.Clock.POSIX (getPOSIXTime) +import Data.IORef (IORef, atomicModifyIORef, newIORef) import Data.Typeable (Typeable) -import Data.Vector.Generic (Vector, unsafeFreeze) -import Data.Word (Word, Word8, Word16, Word32, Word64) -import Foreign.Marshal.Alloc (allocaBytes) -import Foreign.Marshal.Array (peekArray) -import GHC.Base (Int(I#)) -import GHC.Word (Word64(W64#), uncheckedShiftL64#, uncheckedShiftRL64#) -import Prelude hiding (catch) +import Data.Vector.Generic (Vector) +import Data.Word +import Data.Kind 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 (hPutStrLn, stderr) 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 +import qualified Control.Exception as E +import System.Random.MWC.SeedSource +import qualified System.Random.Stateful as Random +#if MIN_VERSION_random(1,3,0) +import Data.List.NonEmpty (NonEmpty(..), toList) +#endif --- | The class of types for which we can generate uniformly +-- | NOTE: Consider use of more principled type classes +-- 'Random.Uniform' and 'Random.UniformRange' instead. +-- +-- The class of types for which we can generate uniformly -- distributed random variates. -- -- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222) @@ -84,7 +194,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: -- @@ -96,9 +206,6 @@ -- statistical calculations that require non-zero values -- (e.g. uses of the 'log' function). -- - -- * The range of random 'Integer' variates is the same as for - -- 'Int'. - -- -- To generate a 'Float' variate with a range of [0,1), subtract -- 2**(-33). To do the same with 'Double' variates, subtract -- 2**(-53). @@ -114,49 +221,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 + uniform = uniform1 id + 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 #-} @@ -182,58 +289,54 @@ {-# INLINE uniformR #-} instance Variate Int where -#if WORD_SIZE_IN_BITS < 64 +#if WORD_SIZE_IN_BITS == 32 uniform = uniform1 fromIntegral -#else +#elif WORD_SIZE_IN_BITS == 64 uniform = uniform2 wordsTo64Bit +#else +#error "Word size is not 32 nor 64" #endif - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} instance Variate Word where -#if WORD_SIZE_IN_BITS < 64 +#if WORD_SIZE_IN_BITS == 32 uniform = uniform1 fromIntegral -#else +#elif WORD_SIZE_IN_BITS == 64 uniform = uniform2 wordsTo64Bit +#else +#error "Word size is not 32 nor 64" #endif - uniformR = uniformRange + uniformR a b = uniformRange a b {-# INLINE uniform #-} {-# INLINE uniformR #-} -{- -instance Variate Integer where - uniform g = do - u <- uniform g - return $! fromIntegral (u :: Int) - {-# INLINE uniform #-} --} - instance (Variate a, Variate b) => Variate (a,b) where - uniform g = (,) `liftM` uniform g `ap` uniform g - uniformR ((x1,y1),(x2,y2)) g = (,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g + uniform g = (,) <$> uniform g <*> uniform g + uniformR ((x1,y1),(x2,y2)) g = (,) <$> uniformR (x1,x2) g <*> uniformR (y1,y2) g {-# INLINE uniform #-} {-# INLINE uniformR #-} instance (Variate a, Variate b, Variate c) => Variate (a,b,c) where - uniform g = (,,) `liftM` uniform g `ap` uniform g `ap` uniform g + uniform g = (,,) <$> uniform g <*> uniform g <*> uniform g uniformR ((x1,y1,z1),(x2,y2,z2)) g = - (,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` uniformR (z1,z2) g + (,,) <$> uniformR (x1,x2) g <*> uniformR (y1,y2) g <*> uniformR (z1,z2) g {-# INLINE uniform #-} {-# INLINE uniformR #-} instance (Variate a, Variate b, Variate c, Variate d) => Variate (a,b,c,d) where - uniform g = (,,,) `liftM` uniform g `ap` uniform g `ap` uniform g - `ap` uniform g + uniform g = (,,,) <$> uniform g <*> uniform g <*> uniform g + <*> uniform g uniformR ((x1,y1,z1,t1),(x2,y2,z2,t2)) g = - (,,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` - uniformR (z1,z2) g `ap` uniformR (t1,t2) g + (,,,) <$> uniformR (x1,x2) g <*> uniformR (y1,y2) g <*> + uniformR (z1,z2) g <*> uniformR (t1,t2) g {-# INLINE uniform #-} {-# INLINE uniformR #-} -wordsTo64Bit :: Integral a => Word32 -> Word32 -> a +wordsTo64Bit :: (Integral a) => Word32 -> Word32 -> a wordsTo64Bit x y = - fromIntegral ((fromIntegral x `shiftL` 32) .|. fromIntegral y) + fromIntegral ((fromIntegral x `shiftL` 32) .|. fromIntegral y :: Word64) {-# INLINE wordsTo64Bit #-} wordToBool :: Word32 -> Bool @@ -249,7 +352,7 @@ wordsToDouble :: Word32 -> Word32 -> Double wordsToDouble x y = (fromIntegral u * m_inv_32 + (0.5 + m_inv_53) + - fromIntegral (v .&. 0xFFFFF) * m_inv_52) + fromIntegral (v .&. 0xFFFFF) * m_inv_52) where m_inv_52 = 2.220446049250313080847263336181640625e-16 m_inv_53 = 1.1102230246251565404236316680908203125e-16 m_inv_32 = 2.3283064365386962890625e-10 @@ -257,15 +360,25 @@ v = fromIntegral y :: Int32 {-# INLINE wordsToDouble #-} --- | State of the pseudo-random number generator. +-- | State of the pseudo-random number generator. It uses mutable +-- state so same generator shouldn't be used from the different +-- threads simultaneously. 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 @@ -282,105 +395,209 @@ -- -- Examples: -- --- > initialize (singletonU 42) +-- > initialize (singleton 42) -- --- > initialize (toU [4, 8, 15, 16, 23, 42]) +-- > initialize (fromList [4, 8, 15, 16, 23, 42]) -- -- If a seed contains fewer than 256 elements, it is first used -- verbatim, then its elements are 'xor'ed against elements of the -- default seed until 256 elements are reached. -initialize :: PrimMonad m => I.Vector Word32 -> m (Gen (PrimState m)) +-- +-- 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'@ +-- +-- In the MWC algorithm, the /carry/ value must be strictly smaller than the +-- multiplicator (see https://en.wikipedia.org/wiki/Multiply-with-carry). +-- Hence, if a seed contains exactly 258 elements, the /carry/ value, which is +-- the last of the 258 values, is moduloed by the multiplicator. +-- +-- Note that if the /first/ carry value is strictly smaller than the multiplicator, +-- all subsequent carry values are also strictly smaller than the multiplicator +-- (a proof of this is in the comments of the code of 'uniformWord32'), hence +-- when restoring a saved state, we have the guarantee that moduloing the saved +-- carry won't modify its value. +initialize :: (PrimMonad m, Vector v Word32) => + v Word32 -> m (Gen (PrimState m)) initialize seed = do q <- M.unsafeNew 258 fill q - M.unsafeWrite q ioff 255 - M.unsafeWrite q coff 362436 + if fini == 258 + then do + M.unsafeWrite q ioff $ G.unsafeIndex seed ioff .&. 255 + M.unsafeWrite q coff $ G.unsafeIndex seed coff `mod` fromIntegral aa + else do + M.unsafeWrite q ioff 255 + M.unsafeWrite q coff 362436 return (Gen q) where fill q = go 0 where go i | i == 256 = return () | otherwise = M.unsafeWrite q i s >> go (i+1) where s | i >= fini = if fini == 0 - then I.unsafeIndex defaultSeed i - else I.unsafeIndex defaultSeed i `xor` - I.unsafeIndex seed (i `mod` fini) - | otherwise = I.unsafeIndex seed i - fini = I.length seed + then G.unsafeIndex defaultSeed i + else G.unsafeIndex defaultSeed i `xor` + G.unsafeIndex seed (i `mod` fini) + | otherwise = G.unsafeIndex seed i + fini = G.length seed {-# INLINE initialize #-} - + -- | An immutable snapshot of the state of a 'Gen'. newtype Seed = Seed (I.Vector Word32) - deriving (Eq, Show, Typeable) + deriving (Eq, Show, Typeable) --- 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' +-- | Convert seed into vector. +fromSeed :: Seed -> I.Vector Word32 +fromSeed (Seed v) = v +-- | @since 0.15.0.0 +instance (s ~ PrimState m, PrimMonad m) => Random.StatefulGen (Gen s) m where + uniformWord32R u = uniformR (0, u) + {-# INLINE uniformWord32R #-} + uniformWord64R u = uniformR (0, u) + {-# INLINE uniformWord64R #-} + uniformWord8 = uniform + {-# INLINE uniformWord8 #-} + uniformWord16 = uniform + {-# INLINE uniformWord16 #-} + uniformWord32 = uniform + {-# INLINE uniformWord32 #-} + uniformWord64 = uniform + {-# INLINE uniformWord64 #-} +#if MIN_VERSION_random(1,3,0) + uniformByteArrayM isPinned n g = stToPrim (Random.fillByteArrayST isPinned n (uniform g)) + {-# INLINE uniformByteArrayM #-} +#else + uniformShortByteString n g = stToPrim (Random.genShortByteStringST n (uniform g)) + {-# INLINE uniformShortByteString #-} +#endif + +-- | @since 0.15.0.0 +instance PrimMonad m => Random.FrozenGen Seed m where + type MutableGen Seed m = Gen (PrimState m) + freezeGen = save +#if MIN_VERSION_random(1,3,0) + modifyGen gen@(Gen mv) f = do + seed <- save gen + case f seed of + (a, Seed v) -> a <$ G.copy mv v + overwriteGen (Gen mv) (Seed v) = G.copy mv v + +instance PrimMonad m => Random.ThawedGen Seed m where +#endif + thawGen = restore + +#if MIN_VERSION_random(1,3,0) +instance Random.SeedGen Seed where + type SeedSize Seed = 1032 -- == 4 * 258 + fromSeed64 seed64 = toSeed $ I.fromListN 258 + [ w32 + | !w64 <- toList seed64 + , !w32 <- [ fromIntegral (w64 `shiftR` 32) + , fromIntegral w64 ] + ] + toSeed64 vSeed = + let w32sToW64 :: Word32 -> Word32 -> Word64 + w32sToW64 w32u w32l = + (fromIntegral w32u `shiftL` 32) .|. fromIntegral w32l + v = fromSeed vSeed + evens = I.ifilter (\i _ -> even i) v + odds = I.ifilter (\i _ -> odd i) v + in case I.toList $ I.zipWith w32sToW64 evens odds of + [] -> + error $ "Impossible: Seed had an unexpected length of: " ++ show (I.length v) + x:xs -> x :| xs +#endif + +-- | Convert vector to 'Seed'. It acts similarly to 'initialize' and +-- will accept any vector. If you want to pass seed immediately to +-- restore you better call initialize directly since following law holds: +-- +-- > restore (toSeed v) = initialize v +toSeed :: (Vector v Word32) => v Word32 -> Seed +toSeed v = Seed $ I.create $ do { Gen q <- initialize v; return q } + -- | 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 <$> G.freeze q {-# INLINE save #-} -- | 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 <$> G.thaw s {-# INLINE restore #-} - --- | Using the current time as a seed, perform an action that uses a --- random variate generator. This is a horrible fallback for Windows --- systems. -withTime :: (PrimMonad m) => (Gen (PrimState m) -> m a) -> IO a -withTime act = do - c <- (numerator . (%cpuTimePrecision)) `liftM` getCPUTime - t <- toRational `liftM` getPOSIXTime - let n = fromIntegral (numerator t) :: Word64 - seed = [fromIntegral c, fromIntegral n, fromIntegral (n `shiftR` 32)] - unsafePrimToIO $ initialize (I.fromList seed) >>= act --- | Seed a PRNG with data from the system's fast source of --- pseudo-random numbers (\"\/dev\/urandom\" on Unix-like systems), --- then run the given action. + +-- $seeding -- --- /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 --- highly independent. -withSystemRandom :: PrimMonad m => (Gen (PrimState m) -> m a) -> IO a -withSystemRandom act = tryRandom `catch` \(_::IOException) -> do - seen <- atomicModifyIORef warned ((,) True) - unless seen $ do - hPutStrLn stderr ("Warning: Couldn't open " ++ show random) +-- Library provides several functions allowing to intialize generator +-- using OS-provided randomness: \"@\/dev\/urandom@\" on Unix-like +-- systems or @RtlGenRandom@ on Windows. 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. + +createSystemRandomList :: IO [Word32] +createSystemRandomList = do + acquireSeedSystem 256 `E.catch` \(_::E.IOException) -> do + seen <- atomicModifyIORef seedCreatetionWarned ((,) True) + unless seen $ E.handle (\(_::E.IOException) -> return ()) $ do + hPutStrLn stderr $ "Warning: Couldn't use randomness source " ++ randomSourceName hPutStrLn stderr ("Warning: using system clock for seed instead " ++ "(quality will be lower)") - withTime act - where tryRandom = do - let nbytes = 1024 - ws <- allocaBytes nbytes $ \buf -> do - nread <- withBinaryFile random ReadMode $ - \h -> hGetBuf h buf nbytes - peekArray (nread `div` 4) buf - unsafePrimToIO $ initialize (I.fromList ws) >>= act - random = "/dev/urandom" - warned = unsafePerformIO $ newIORef False - {-# NOINLINE warned #-} + acquireSeedTime --- | Unchecked 64-bit left shift. -shiftL :: Word64 -> Int -> Word64 -shiftL (W64# x#) (I# i#) = W64# (x# `uncheckedShiftL64#` i#) +seedCreatetionWarned :: IORef Bool +seedCreatetionWarned = unsafePerformIO $ newIORef False +{-# NOINLINE seedCreatetionWarned #-} --- | Unchecked 64-bit right shift. -shiftR :: Word64 -> Int -> Word64 -shiftR (W64# x#) (I# i#) = W64# (x# `uncheckedShiftRL64#` i#) + +-- | Generate random seed for generator using system's fast source of +-- pseudo-random numbers. +-- +-- @since 0.15.0.0 +createSystemSeed :: IO Seed +createSystemSeed = do + seed <- createSystemRandomList + return $! toSeed $ I.fromList seed + +-- | Seed a PRNG with data from the system's fast source of +-- pseudo-random numbers. +createSystemRandom :: IO GenIO +createSystemRandom = initialize . I.fromList =<< createSystemRandomList + + +-- | Seed PRNG with data from the system's fast source of +-- pseudo-random numbers and execute computation in ST monad. +-- +-- @since 0.15.0.0 +withSystemRandomST :: (forall s. Gen s -> ST s a) -> IO a +withSystemRandomST act = do + seed <- createSystemSeed + return $! runST $ act =<< restore seed + +-- | Seed a PRNG with data from the system's fast source of +-- pseudo-random numbers, then run the given action. +-- +-- This function is unsafe and for example allows STRefs or any +-- other mutable data structure to escape scope: +-- +-- >>> ref <- withSystemRandom $ \_ -> newSTRef 1 +-- >>> withSystemRandom $ \_ -> modifySTRef ref succ >> readSTRef ref +-- 2 +-- >>> withSystemRandom $ \_ -> modifySTRef ref succ >> readSTRef ref +-- 3 +withSystemRandom :: PrimBase m + => (Gen (PrimState m) -> m a) -> IO a +withSystemRandom act = do + seed <- createSystemSeed + unsafePrimToIO $ act =<< restore seed +{-# DEPRECATED withSystemRandom "Use withSystemRandomST or createSystemSeed or createSystemRandom instead" #-} + + -- | Compute the next index into the state pool. This is simply -- addition modulo 256. nextIndex :: Integral a => a -> Int @@ -388,23 +605,28 @@ where j = fromIntegral (i+1) :: Word8 {-# INLINE nextIndex #-} -a :: Word64 -a = 1540315826 -{-# INLINE a #-} +-- The multiplicator : 0x5BCF5AB2 +-- +-- Eventhough it is a 'Word64', it is important for the correctness of the proof +-- on carry value that it is /not/ greater than maxBound 'Word32'. +aa :: Word64 +aa = 1540315826 +{-# INLINE aa #-} uniformWord32 :: PrimMonad m => Gen (PrimState m) -> m Word32 +-- NOTE [Carry value] 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 + i <- nextIndex <$> M.unsafeRead q ioff + c <- fromIntegral <$> M.unsafeRead q coff + qi <- fromIntegral <$> M.unsafeRead q i + 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' #) M.unsafeWrite q i x' M.unsafeWrite q ioff (fromIntegral i) - M.unsafeWrite q coff (fromIntegral c'') + M.unsafeWrite q coff c'' return x' {-# INLINE uniformWord32 #-} @@ -416,17 +638,17 @@ uniform2 :: PrimMonad m => (Word32 -> Word32 -> a) -> Gen (PrimState m) -> m a uniform2 f (Gen q) = do - i <- nextIndex `liftM` M.unsafeRead q ioff + i <- nextIndex <$> M.unsafeRead q ioff let j = nextIndex i - 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 + c <- fromIntegral <$> M.unsafeRead q coff + qi <- fromIntegral <$> M.unsafeRead q i + qj <- fromIntegral <$> M.unsafeRead q j + 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 #) @@ -434,121 +656,87 @@ M.unsafeWrite q i x' M.unsafeWrite q j y' M.unsafeWrite q ioff (fromIntegral j) - M.unsafeWrite q coff (fromIntegral d'') + M.unsafeWrite q coff d'' return $! f x' y' {-# INLINE uniform2 #-} -- Type family for fixed size integrals. For signed data types it's --- its unsigned couterpart with same size and for unsigned data types +-- its unsigned counterpart with same size and for unsigned data types -- it's same type -type family Unsigned a :: * +type family Unsigned a :: Type type instance Unsigned Int8 = Word8 type instance Unsigned Int16 = Word16 type instance Unsigned Int32 = Word32 type instance Unsigned Int64 = Word64 -type instance Unsigned Int = Word type instance Unsigned Word8 = Word8 type instance Unsigned Word16 = Word16 type instance Unsigned Word32 = Word32 type instance Unsigned Word64 = Word64 -type instance Unsigned Word = Word +type instance Unsigned Int = Word +type instance Unsigned Word = Word + + -- Subtract two numbers under assumption that x>=y and store result in -- 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) -> Int -> m (I.Vector a) +uniformVector + :: (PrimMonad m, Random.StatefulGen g m, Random.Uniform a, Vector v a) + => g -> Int -> m (v a) +-- NOTE: We use in-place mutation in order to generate vector instead +-- of generateM because latter will go though intermediate list until +-- we're working in IO/ST monad +-- +-- See: https://github.com/haskell/vector/issues/208 for details uniformVector gen n = do - mu <- M.unsafeNew n - let go !i | i < n = uniform gen >>= M.unsafeWrite mu i >> go (i+1) - | otherwise = unsafeFreeze mu + mu <- GM.unsafeNew n + let go !i | i < n = Random.uniformM gen >>= GM.unsafeWrite mu i >> go (i+1) + | otherwise = G.unsafeFreeze mu go 0 {-# INLINE uniformVector #-} -data T = T {-# UNPACK #-} !Double {-# UNPACK #-} !Double --- | Generate a normally distributed random variate. --- --- 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. -normal :: PrimMonad m => 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 #-} - +-- This is default seed for the generator and used when no seed is +-- specified or seed is only partial. It's not known how it was +-- generated but it looks random enough defaultSeed :: I.Vector Word32 defaultSeed = I.fromList [ 0x7042e8b3, 0x06f7f4c5, 0x789ea382, 0x6fb15ad8, 0x54f7a879, 0x0474b184, @@ -598,20 +786,89 @@ -- $references -- --- * Doornik, J.A. (2005) An improved ziggurat method to generate --- normal random samples. Mimeo, Nuffield College, University of --- Oxford. +-- * Marsaglia, G. (2003) Seeds for random number generators. +-- /Communications of the ACM/ 46(5):90–93. +-- -- -- * Doornik, J.A. (2007) Conversion of high-period random numbers to -- floating point. -- /ACM Transactions on Modeling and Computer Simulation/ 17(1). -- + + +-- $typehelp -- --- * Marsaglia, G. (2003) Seeds for random number generators. --- /Communications of the ACM/ 46(5):90–93. --- +-- The functions in this package are deliberately written for +-- flexibility, and will run in both the 'IO' and 'ST' monads. -- --- * Thomas, D.B.; Leong, P.G.W.; Luk, W.; Villasenor, J.D. --- (2007). Gaussian random number generators. --- /ACM Computing Surveys/ 39(4). --- +-- 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. + + + +-- $setup +-- +-- >>> import Control.Monad +-- >>> import Data.Word +-- >>> import Data.STRef +-- >>> :set -Wno-deprecations + + +-- NOTE [Carry value] +-- ------------------ +-- This is proof of statement: +-- +-- > if the carry value is strictly smaller than the multiplicator, +-- > the next carry value is also strictly smaller than the multiplicator. +-- +-- Even though the proof is written in terms of the actual value of the +-- multiplicator, it holds for any multiplicator value /not/ greater +-- than maxBound 'Word32' +-- +-- (In the code, the multiplicator is aa, the carry value is c, +-- the next carry value is c''.) +-- +-- So we'll assume that c < aa, and show that c'' < aa : +-- +-- by definition, aa = 0x5BCF5AB2, qi <= 0xFFFFFFFF (because it is a 'Word32') +-- +-- Then we get following: +-- +-- aa*qi <= 0x5BCF5AB200000000 - 0x5BCF5AB2. +-- t < 0x5BCF5AB200000000 (because t = aa * qi + c and c < 0x5BCF5AB2) +-- t <= 0x5BCF5AB1FFFFFFFF +-- c' < 0x5BCF5AB1 +-- c'' < 0x5BCF5AB2, +-- c'' < aa, which is what we wanted to prove. 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,285 @@ +{-# 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 + -- * References + -- $references + ) where + +import Control.Arrow (second,(***)) + +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 Numeric.SpecFunctions (logFactorial) +import System.Random.Stateful + +import Prelude hiding ((++)) + + + +-- | 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 :: (StatefulGen g m, Vector v a) => CondensedTable v a -> g -> m a +{-# INLINE genFromTable #-} +genFromTable table gen = do + w <- uniformM 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 +-- probabilities 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 tbl = pkgError "tableFromProbabilities" "empty vector of outcomes" + | otherwise = tableFromIntWeights $ G.map (second $ toWeight . (* mlt)) tbl + where + -- 2^32. N.B. This number is exatly representable. + mlt = 4.294967296e9 + -- Drop non-positive probabilities + tbl = G.filter ((> 0) . snd) v + -- Convert Double weight to Word32 and avoid overflow at the same + -- time. It's especially dangerous if one probability is + -- approximately 1 and others are 0. + toWeight w | w > mlt - 1 = 2^(32::Int) - 1 + | otherwise = round w + + +-- | 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@ at least approximately. This function will +-- correct small deviations from @2^32@ such as arising from rounding +-- errors. But for large deviations it's likely to product incorrect +-- result with terrible performance. +tableFromIntWeights :: (Vector v (a,Word32), Vector v a, Vector v Word32) + => v (a, Word32) + -> CondensedTable v a +{-# INLINE tableFromIntWeights #-} +tableFromIntWeights v + | n == 0 = pkgError "tableFromIntWeights" "empty table" + -- Single element tables should be treated separately. 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 + -- We must filter out zero-probability outcomes because they may + -- confuse weight correction algorithm + tbl = G.filter ((/=0) . snd) v + 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 distribution. 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 + -- + -- FIXME: this is not ideal precision-wise. Check if code + -- from statistics gives better precision. + nMax = floor lam :: Int + pMax = exp $ fromIntegral nMax * log lam - lam - logFactorial 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) + ) + -- Minimal representable probability for condensed tables + 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 new file mode 100644 --- /dev/null +++ b/System/Random/MWC/Distributions.hs @@ -0,0 +1,589 @@ +{-# LANGUAGE MultiWayIf #-} +{-# LANGUAGE BangPatterns, GADTs, FlexibleContexts, ScopedTypeVariables #-} +-- | +-- 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 + -- ** Continuous distributions + normal + , standard + , exponential + , truncatedExp + , gamma + , chiSquare + , beta + -- ** Discrete distribution + , categorical + , logCategorical + , geometric0 + , geometric1 + , bernoulli + , binomial + , poisson + -- ** Multivariate + , dirichlet + -- * Permutations + , uniformPermutation + , uniformShuffle + , uniformShuffleM + -- * References + -- $references + ) where + +import Prelude hiding (mapM) +import Control.Monad.Primitive (PrimMonad, PrimState) +import Data.Bits ((.&.)) +import Data.Foldable (foldl') +import Data.Traversable (mapM) +import Data.Word (Word32) +import System.Random.Stateful (StatefulGen(..),Uniform(..),UniformRange(..),uniformDoublePositive01M, uniformDouble01M) +import qualified Data.Vector.Unboxed as I +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Generic.Mutable as M +import Numeric.SpecFunctions (logFactorial) + +-- Unboxed 2-tuple +data T = T {-# UNPACK #-} !Double {-# UNPACK #-} !Double + + +-- | Generate a normally distributed random variate with given mean +-- and standard deviation. +normal :: StatefulGen g m + => Double -- ^ Mean + -> Double -- ^ Standard deviation + -> g + -> 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 :: StatefulGen g m => g -> m Double +{-# INLINE standard #-} +standard gen = loop + where + loop = do + u <- subtract 1 . (*2) <$> uniformDoublePositive01M gen + ri <- uniformM 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 <- uniformDoublePositive01M gen + if e + c * (d - e) < 1 + then return x + else loop + normalTail neg = tailing + where tailing = do + x <- (/ rNorm) . log <$> uniformDoublePositive01M gen + y <- log <$> uniformDoublePositive01M gen + if y * (-2) < x * x + then tailing + else return $! if neg then x - rNorm else rNorm - x + +-- Constants used by standard/normal. They are floated to the top +-- level to avoid performance regression (Bug #16) when blocks/ratios +-- are recalculated on each call to standard/normal. It's also +-- somewhat difficult to trigger reliably. +blocks :: I.Vector Double +blocks = (`I.snoc` 0) . I.cons (v/f) . I.cons rNorm . I.unfoldrN 126 go $! T rNorm 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 * rNorm * rNorm) +{-# NOINLINE blocks #-} + +rNorm :: Double +rNorm = 3.442619855899 + +ratios :: I.Vector Double +ratios = I.zipWith (/) (I.tail blocks) blocks +{-# NOINLINE ratios #-} + + + +-- | Generate an exponentially distributed random variate. +exponential :: StatefulGen g m + => Double -- ^ Scale parameter + -> g -- ^ Generator + -> m Double +{-# INLINE exponential #-} +exponential b gen = do + x <- uniformDoublePositive01M gen + return $! - log x / b + + +-- | Generate truncated exponentially distributed random variate. +truncatedExp :: StatefulGen g m + => Double -- ^ Scale parameter + -> (Double,Double) -- ^ Range to which distribution is + -- truncated. Values may be negative. + -> g -- ^ Generator. + -> m Double +{-# INLINE truncatedExp #-} +truncatedExp scale (a,b) gen = do + -- We shift a to 0 and then generate distribution truncated to [0,b-a] + -- It's easier + let delta = b - a + p <- uniformDoublePositive01M gen + return $! a - log ( (1 - p) + p*exp(-scale*delta)) / scale + +-- | Random variate generator for gamma distribution. +gamma :: (StatefulGen g m) + => Double -- ^ Shape parameter + -> Double -- ^ Scale parameter + -> g -- ^ 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 <- uniformDoublePositive01M 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 <- uniformDoublePositive01M 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 the chi square distribution. +chiSquare :: StatefulGen g m + => Int -- ^ Number of degrees of freedom + -> g -- ^ 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 + +-- | Random variate generator for the geometric distribution, +-- computing the number of failures before success. Distribution's +-- support is [0..]. +geometric0 :: StatefulGen g m + => Double -- ^ /p/ success probability lies in (0,1] + -> g -- ^ Generator + -> m Int +{-# INLINE geometric0 #-} +geometric0 p gen + | p == 1 = return 0 + | p > 0 && p < 1 = do q <- uniformDoublePositive01M gen + -- FIXME: We want to use log1p here but it will + -- introduce dependency on math-functions. + return $! floor $ log q / log (1 - p) + | otherwise = pkgError "geometric0" "probability out of [0,1] range" + +-- | Random variate generator for geometric distribution for number of +-- trials. Distribution's support is [1..] (i.e. just 'geometric0' +-- shifted by 1). +geometric1 :: StatefulGen g m + => Double -- ^ /p/ success probability lies in (0,1] + -> g -- ^ Generator + -> m Int +{-# INLINE geometric1 #-} +geometric1 p gen = do n <- geometric0 p gen + return $! n + 1 + +-- | Random variate generator for Beta distribution +beta :: StatefulGen g m + => Double -- ^ alpha (>0) + -> Double -- ^ beta (>0) + -> g -- ^ Generator + -> m Double +{-# INLINE beta #-} +beta a b gen = do + x <- gamma a 1 gen + y <- gamma b 1 gen + return $! x / (x+y) + +-- | Random variate generator for Dirichlet distribution +dirichlet :: (StatefulGen g m, Traversable t) + => t Double -- ^ container of parameters + -> g -- ^ Generator + -> m (t Double) +{-# INLINE dirichlet #-} +dirichlet t gen = do + t' <- mapM (\x -> gamma x 1 gen) t + let total = foldl' (+) 0 t' + return $ fmap (/total) t' + +-- | Random variate generator for Bernoulli distribution +bernoulli :: StatefulGen g m + => Double -- ^ Probability of success (returning True) + -> g -- ^ Generator + -> m Bool +{-# INLINE bernoulli #-} +bernoulli p gen = (< p) <$> uniformDoublePositive01M gen + +-- | Random variate generator for categorical distribution. +-- +-- Note that if you need to generate a lot of variates functions +-- "System.Random.MWC.CondensedTable" will offer better +-- performance. If only few is needed this function will faster +-- since it avoids costs of setting up table. +categorical :: (StatefulGen g m, G.Vector v Double) + => v Double -- ^ List of weights [>0] + -> g -- ^ Generator + -> m Int +{-# INLINE categorical #-} +categorical v gen + | G.null v = pkgError "categorical" "empty weights!" + | otherwise = do + let cv = G.scanl1' (+) v + p <- (G.last cv *) <$> uniformDoublePositive01M gen + return $! case G.findIndex (>=p) cv of + Just i -> i + Nothing -> pkgError "categorical" "bad weights!" + +-- | Random variate generator for categorical distribution where the +-- weights are in the log domain. It's implemented in terms of +-- 'categorical'. +logCategorical :: (StatefulGen g m, G.Vector v Double) + => v Double -- ^ List of logarithms of weights + -> g -- ^ Generator + -> m Int +{-# INLINE logCategorical #-} +logCategorical v gen + | G.null v = pkgError "logCategorical" "empty weights!" + | otherwise = categorical (G.map (exp . subtract m) v) gen + where + m = G.maximum v + +-- | Random variate generator for uniformly distributed permutations. +-- It returns random permutation of vector /[0 .. n-1]/. +-- +-- This is the Fisher-Yates shuffle +uniformPermutation :: forall g m v. (StatefulGen g m, PrimMonad m, G.Vector v Int) + => Int + -> g + -> m (v Int) +{-# INLINE uniformPermutation #-} +uniformPermutation n gen + | n < 0 = pkgError "uniformPermutation" "size must be >=0" + | otherwise = uniformShuffle (G.generate n id :: v Int) gen + +-- | Random variate generator for a uniformly distributed shuffle (all +-- shuffles are equiprobable) of a vector. It uses Fisher-Yates +-- shuffle algorithm. +uniformShuffle :: (StatefulGen g m, PrimMonad m, G.Vector v a) + => v a + -> g + -> m (v a) +{-# INLINE uniformShuffle #-} +uniformShuffle vec gen + | G.length vec <= 1 = return vec + | otherwise = do + mvec <- G.thaw vec + uniformShuffleM mvec gen + G.unsafeFreeze mvec + +-- | In-place uniformly distributed shuffle (all shuffles are +-- equiprobable)of a vector. +uniformShuffleM :: (StatefulGen g m, PrimMonad m, M.MVector v a) + => v (PrimState m) a + -> g + -> m () +{-# INLINE uniformShuffleM #-} +uniformShuffleM vec gen + | M.length vec <= 1 = return () + | otherwise = loop 0 + where + n = M.length vec + lst = n-1 + loop i | i == lst = return () + | otherwise = do j <- uniformRM (i,lst) gen + M.unsafeSwap vec i j + loop (i+1) + + +sqr :: Double -> Double +sqr x = x * x +{-# INLINE sqr #-} + +pkgError :: String -> String -> a +pkgError func msg = error $ "System.Random.MWC.Distributions." ++ func ++ + ": " ++ msg + +-- | Random variate generator for Binomial distribution. Will throw +-- exception when parameters are out range. +-- +-- The probability of getting exactly k successes in n trials is +-- given by the probability mass function: +-- +-- \[ +-- f(k;n,p) = \Pr(X = k) = \binom n k p^k(1-p)^{n-k} +-- \] +-- +-- @since 0.15.1.0 +binomial :: forall g m . StatefulGen g m + => Int -- ^ Number of trials, must be positive. + -> Double -- ^ Probability of success \(p \in [0,1]\) + -> g -- ^ Generator + -> m Int +{-# INLINE binomial #-} +binomial n p gen + | n <= 0 = pkgError "binomial" "number of trials must be positive" + | p < 0.0 || p > 1.0 = pkgError "binomial" "probability must be >= 0 and <= 1" + | p == 0.0 = return 0 + | p == 1.0 = return n + | p <= 0.5 = if + | fromIntegral n * p < inv_thr -> binomialInv n p gen + | otherwise -> binomialTPE n p gen + | p > 0.5 = do + ix <- case 1 - p of + p' | fromIntegral n * p' < inv_thr -> binomialInv n p' gen + | otherwise -> binomialTPE n p' gen + pure $! n - ix + -- Reachable when p is NaN + | otherwise = pkgError "binomial" "probability must be >= 0 and <= 1" + where + -- Threshold for preferring the BINV algorithm / inverse cdf + -- logic. The paper suggests 10, Ranlib uses 30, R uses 30, Rust uses + -- 10 and GSL uses 14. + inv_thr = 10 + +-- Binomial-Triangle-Parallelogram-Exponential algorithm (BTPE) +-- described in Kachitvichyanukul1988 +binomialTPE :: forall g m . StatefulGen g m => Int -> Double -> g -> m Int +{-# INLINE binomialTPE #-} +binomialTPE n p g = loop + where + -- Main accept/reject loop + loop = do + u <- uniformRM (0.0, p4) g + v <- uniformDoublePositive01M g + selectArea u v + -- Acceptance / rejection of sample [step 5] + acceptReject :: Int -> Double -> m Int + acceptReject !ix !v + | var <= accept = return ix + | otherwise = loop + where + var = log v + accept = logFactorial bigM + logFactorial (n - bigM) - + logFactorial ix - logFactorial (n - ix) + + fromIntegral (ix - bigM) * log (p / q) + -- Select area to be used [Steps 1-4] + selectArea :: Double -> Double -> m Int + selectArea !u !v + -- Triangular region + | u <= p1 = return $! floor $ xm - p1 * v + u + -- Parallelogram region + | u <= p2 = do let x = xl + (u - p1) / c + w = v * c + 1.0 - abs (x - xm) / p1 + if w > 1 || w <= 0 + then loop + else do let ix = floor x + acceptReject ix w + -- Left tail + | u <= p3 = case floor $ xl + log v / lambdaL of + ix | ix < 0 -> loop + | otherwise -> do let w = v * (u - p2) * lambdaL + acceptReject ix w + -- Right tail + | otherwise = case floor $ xr - log v / lambdaR of + ix | ix > n -> loop + | otherwise -> do let w = v * (u - p3) * lambdaR + acceptReject ix w + ---------------------------------------- + -- Constants used in algorithm. See [Step 0] + q = 1 - p + np = fromIntegral n * p + ffm = np + p + bigM = floor ffm + -- Half integer mean (tip of triangle) + xm = fromIntegral bigM + 0.5 + -- p1: the distance to the left and right edges of the triangle + -- region below the target distribution; since height=1, also: + -- area of region (half base * height) + !p1 = let npq = np * q + in fromIntegral (floor (2.195 * sqrt npq - 4.6 * q) :: Int) + 0.5 + xl = xm - p1 -- Left edge of triangle + xr = xm + p1 -- Right edge of triangle + c = 0.134 + 20.5 / (15.3 + fromIntegral bigM) + -- p1 + area of parallelogram region + !p2 = p1 * (1.0 + c + c) + -- + lambdaL = let al = (ffm - xl) / (ffm - xl * p) + in al * (1.0 + 0.5 * al) + lambdaR = let ar = (xr - ffm) / (xr * q) + in ar * (1.0 + 0.5 * ar) + -- p2 + area of left tail + !p3 = p2 + c / lambdaL + -- p3 + area of right tail + !p4 = p3 + c / lambdaR + + +-- Compute binomial variate using inversion method (BINV in +-- Kachitvichyanukul1988) +binomialInv :: StatefulGen g m => Int -> Double -> g -> m Int +{-# INLINE binomialInv #-} +binomialInv n p g = do + u <- uniformDoublePositive01M g + return $! invertBinomial n p u + +-- This function is defined on top level to avoid inlining it since it's rather +-- large and we don't need specializations since it's monomorphic anyway +invertBinomial + :: Int -- N of trials + -> Double -- probability of success + -> Double -- Output of PRNG + -> Int +invertBinomial !n !p !u0 = invert (q^n) u0 0 + where + -- We forcing s&a in order to avoid allocating thunks. Those are + -- more expensive than computing them unconditionally + q = 1 - p + !s = p / q + !a = fromIntegral (n + 1) * s + -- + invert !r !u !x + | u <= r = x + | otherwise = invert r' u' x' + where + u' = u - r + x' = x + 1 + r' = r * ((a / fromIntegral x') - s) + + +-- | Random variate generate for Poisson distribution. +-- +-- If parameter λ is within 10 σ or greater than @maxBound :: Int@ +-- error is raised since result may not be representable as @Int@ +-- +-- @since 0.15.3.0 +poisson + :: StatefulGen g m + => Double -- ^ Rate parameter, also known as \( \lambda \) + -> g -- ^ Generator + -> m Int +{-# INLINE poisson #-} +poisson lambda gen + | lambda > maxPoissonLam + = pkgError "poisson" + $ "lambda is too large: "++show lambda++" Result won't fit into Int" + | lambda > 10 = poissonAtkinson lambda gen + | lambda >= 0 = poissonInterArrival lambda gen + | otherwise + = pkgError "poisson" "Lambda parameter must be greater than zero" + +maxPoissonLam :: Double +maxPoissonLam = m - 10 * sqrt m where + m = fromIntegral (maxBound :: Int) + + +-- This uses the fact that if N(t) is a Poisson process +-- with rate lambda, then the counting process N(t) can +-- be represented as interarrival times X[1], X[2],... with +-- X[i] ~ Exp(lambda). +poissonInterArrival :: StatefulGen g m + => Double -- ^ Rate parameter, also known as lambda + -> g -- ^ Generator + -> m Int +{-# INLINE poissonInterArrival #-} +poissonInterArrival lambda gen = do + loop 0 1.0 + where + loop !k !p = do + p' <- (*p) <$> uniformDouble01M gen + if p' <= bigL then return $! k else loop (k+1) p' + bigL = exp (negate lambda) + +-- Attributed to Atkinson, via Casella. Uses a rejection +-- method that uses logistic distribution as the envelope +-- distribution. +poissonAtkinson :: forall g m . StatefulGen g m + => Double -- ^ Rate parameter, also known as lambda + -> g -- ^ Generator + -> m Int +{-# INLINE poissonAtkinson #-} +poissonAtkinson lambda gen = loop + where loop :: m Int + loop = do + bigU <- uniformDouble01M gen + let x = (alpha - log ((1.0 - bigU) / bigU)) / bbeta + if x < (-0.5) + then loop + else do + bigV <- uniformDouble01M gen + let n = floor (x + 0.5) :: Int + y = alpha - bbeta * x + logFacN = logFactorial n + lhs = y + log (bigV / (1.0 + exp y)**2) + rhs = bigK + fromIntegral n * logLambda - logFacN + if lhs <= rhs + then return n + else loop + bigC,alpha,bbeta,bigK,logLambda :: Double + bigC = 0.767 - 3.36 / lambda + bbeta = pi / sqrt (3.0 * lambda) + alpha = bbeta * lambda + bigK = log bigC - lambda - log bbeta + logLambda = log lambda + +-- $references +-- +-- * Doornik, J.A. (2005) An improved ziggurat method to generate +-- normal random samples. Mimeo, Nuffield College, University of +-- Oxford. +-- +-- * Thomas, D.B.; Leong, P.G.W.; Luk, W.; Villasenor, J.D. +-- (2007). Gaussian random number generators. +-- /ACM Computing Surveys/ 39(4). +-- +-- +-- * Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random +-- Variate Generation. Communications of the ACM, 31, 2 (February, +-- 1988) 216. +-- Here's an example of how the algorithm's sampling regions look +-- ![Something](docs/RecreateFigure.svg) +-- +-- * Devroye, L. (1986) Non-uniform Random Variate Generation. +-- Springer Verlag. Chapter 10: Discrete Univariate Distributions. +-- +-- +-- * Robert, C.P. & Casella, G. Monte Carlo Statistical Methods. +-- Springer Texts in Statistics. Algorithm A6: Atkinson's Method +-- for Generating Poisson Random Variables. +-- diff --git a/System/Random/MWC/SeedSource.hs b/System/Random/MWC/SeedSource.hs new file mode 100644 --- /dev/null +++ b/System/Random/MWC/SeedSource.hs @@ -0,0 +1,97 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE ScopedTypeVariables #-} +-- | +-- Low level source of random values for seeds. It should work on both +-- unices and windows +module System.Random.MWC.SeedSource ( + acquireSeedSystem + , acquireSeedTime + , randomSourceName + ) where + +import Data.Word (Word32,Word64) +import Data.Bits (shiftR) +import Data.Ratio ((%), numerator) +import Data.Time.Clock.POSIX (getPOSIXTime) + +import Foreign.Storable +import Foreign.Marshal.Alloc (allocaBytes) +import Foreign.Marshal.Array (peekArray) +#if defined(mingw32_HOST_OS) +import Foreign.Ptr +import Foreign.C.Types +#else +import System.IO (IOMode(..), hGetBuf, withBinaryFile) +#endif +import System.CPUTime (cpuTimePrecision, getCPUTime) + +-- Acquire seed from current time. This is horrible fallback for +-- Windows system. +acquireSeedTime :: IO [Word32] +acquireSeedTime = do + c <- numerator . (% cpuTimePrecision) <$> getCPUTime + t <- toRational <$> getPOSIXTime + let n = fromIntegral (numerator t) :: Word64 + return [fromIntegral c, fromIntegral n, fromIntegral (n `shiftR` 32)] + +-- | Acquire seed from the system entropy source. On Unix machines, +-- this will attempt to use @/dev/urandom@. On Windows, it will internally +-- use @RtlGenRandom@. +acquireSeedSystem :: forall a. Storable a => Int -> IO [a] +acquireSeedSystem nElts = do + let eltSize = sizeOf (undefined :: a) + nbytes = nElts * eltSize +#if !defined(mingw32_HOST_OS) + allocaBytes nbytes $ \buf -> do + nread <- withBinaryFile "/dev/urandom" ReadMode $ \h -> hGetBuf h buf nbytes + peekArray (nread `div` eltSize) buf +#else + -- Generate 256 random Word32s from RtlGenRandom + allocaBytes nbytes $ \buf -> do + ok <- c_RtlGenRandom buf (fromIntegral nbytes) + if ok then return () else fail "Couldn't use RtlGenRandom" + peekArray nElts buf + +-- Note: on 64-bit Windows, the 'stdcall' calling convention +-- isn't supported, so we use 'ccall' instead. +#if defined(i386_HOST_ARCH) +# define WINDOWS_CCONV stdcall +#elif defined(x86_64_HOST_ARCH) +# define WINDOWS_CCONV ccall +#else +# error Unknown mingw32 architecture! +#endif + +-- Note: On Windows, the typical convention would be to use +-- the CryptoGenRandom API in order to generate random data. +-- However, here we use 'SystemFunction036', AKA RtlGenRandom. +-- +-- This is a commonly used API for this purpose; one bonus is +-- that it avoids having to bring in the CryptoAPI library, +-- and completely sidesteps the initialization cost of CryptoAPI. +-- +-- While this function is technically "subject to change" that is +-- extremely unlikely in practice: rand_s in the Microsoft CRT uses +-- this, and they can't change it easily without also breaking +-- backwards compatibility with e.g. statically linked applications. +-- +-- The name 'SystemFunction036' is the actual link-time name; the +-- display name is just for giggles, I guess. +-- +-- See also: +-- - http://blogs.msdn.com/b/michael_howard/archive/2005/01/14/353379.aspx +-- - https://bugzilla.mozilla.org/show_bug.cgi?id=504270 +-- +foreign import WINDOWS_CCONV unsafe "SystemFunction036" + c_RtlGenRandom :: Ptr a -> CULong -> IO Bool +#endif + + +-- | Name of source of randomness. It should be used in error messages +randomSourceName :: String +#if !defined(mingw32_HOST_OS) +randomSourceName = "/dev/urandom" +#else +randomSourceName = "RtlGenRandom" +#endif diff --git a/bench-papi/Bench.hs b/bench-papi/Bench.hs new file mode 100644 --- /dev/null +++ b/bench-papi/Bench.hs @@ -0,0 +1,14 @@ +-- | +-- Here we reexport definitions of tasty-bench +module Bench + ( whnf + , nf + , nfIO + , whnfIO + , bench + , bgroup + , defaultMain + , benchIngredients + ) where + +import Test.Tasty.PAPI diff --git a/bench-time/Bench.hs b/bench-time/Bench.hs new file mode 100644 --- /dev/null +++ b/bench-time/Bench.hs @@ -0,0 +1,14 @@ +-- | +-- Here we reexport definitions of tasty-bench +module Bench + ( whnf + , nf + , nfIO + , whnfIO + , bench + , bgroup + , defaultMain + , benchIngredients + ) where + +import Test.Tasty.Bench diff --git a/bench/Benchmark.hs b/bench/Benchmark.hs new file mode 100644 --- /dev/null +++ b/bench/Benchmark.hs @@ -0,0 +1,161 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +module Main(main) where + +import Control.Exception +import Data.Int +import Data.Word +import Data.Proxy +import qualified Data.Vector.Unboxed as U +import qualified System.Random as R +import System.Random.Stateful (StatefulGen) +import System.Random.MWC +import System.Random.MWC.Distributions +import System.Random.MWC.CondensedTable +import qualified System.Random.Mersenne as M + +import Test.Tasty.Options +import Test.Tasty.Runners +import Test.Tasty (includingOptions) +import Bench + +-- | Size of vector used in benchmarks +newtype Iterations = Iterations Int + +instance IsOption Iterations where + defaultValue = Iterations 10000 + parseValue = fmap Iterations . safeRead + optionName = pure "iter" + optionHelp = pure "Number of iteration in sampling benchmarks" + + +loop :: Iterations -> IO a -> IO () +{-# INLINE loop #-} +loop (Iterations n) act = go n where + go i | i <= 0 = pure () + | otherwise = do _ <- evaluate =<< act + go (i - 1) + +makeTableUniform :: Int -> CondensedTable U.Vector Int +makeTableUniform n = + tableFromProbabilities $ U.zip (U.enumFromN 0 n) (U.replicate n (1 / fromIntegral n)) +{-# INLINE makeTableUniform #-} + +main :: IO () +main = do + -- Set up tasty + let tasty_opts = [Option (Proxy :: Proxy Iterations)] + ingredients = includingOptions tasty_opts : benchIngredients + opts <- parseOptions ingredients (bench "Fake" (nf id ())) + let iter = lookupOption opts + -- Set up RNG + mwc <- create + seed <- save mwc + mtg <- M.newMTGen . Just =<< uniform mwc + defaultMainWithIngredients ingredients $ bgroup "All" + [ bgroup "mwc" + -- One letter group names are used so they will fit on the plot. + -- + -- U - uniform + -- R - uniformR + -- D - distribution + [ bgroup "U" + [ bench "Double" $ whnfIO $ loop iter (uniform mwc :: IO Double) + , bench "Int" $ whnfIO $ loop iter (uniform mwc :: IO Int) + , bench "Int8" $ whnfIO $ loop iter (uniform mwc :: IO Int8) + , bench "Int16" $ whnfIO $ loop iter (uniform mwc :: IO Int16) + , bench "Int32" $ whnfIO $ loop iter (uniform mwc :: IO Int32) + , bench "Int64" $ whnfIO $ loop iter (uniform mwc :: IO Int64) + , bench "Word" $ whnfIO $ loop iter (uniform mwc :: IO Word) + , bench "Word8" $ whnfIO $ loop iter (uniform mwc :: IO Word8) + , bench "Word16" $ whnfIO $ loop iter (uniform mwc :: IO Word16) + , bench "Word32" $ whnfIO $ loop iter (uniform mwc :: IO Word32) + , bench "Word64" $ whnfIO $ loop iter (uniform mwc :: IO Word64) + ] + , bgroup "R" + -- I'm not entirely convinced that this is right way to test + -- uniformR. /A.Khudyakov/ + [ bench "Double" $ whnfIO $ loop iter (uniformR (-3.21,26) mwc :: IO Double) + , bench "Int" $ whnfIO $ loop iter (uniformR (-12,679) mwc :: IO Int) + , bench "Int8" $ whnfIO $ loop iter (uniformR (-12,4) mwc :: IO Int8) + , bench "Int16" $ whnfIO $ loop iter (uniformR (-12,679) mwc :: IO Int16) + , bench "Int32" $ whnfIO $ loop iter (uniformR (-12,679) mwc :: IO Int32) + , bench "Int64" $ whnfIO $ loop iter (uniformR (-12,679) mwc :: IO Int64) + , bench "Word" $ whnfIO $ loop iter (uniformR (34,633) mwc :: IO Word) + , bench "Word8" $ whnfIO $ loop iter (uniformR (34,63) mwc :: IO Word8) + , bench "Word16" $ whnfIO $ loop iter (uniformR (34,633) mwc :: IO Word16) + , bench "Word32" $ whnfIO $ loop iter (uniformR (34,633) mwc :: IO Word32) + , bench "Word64" $ whnfIO $ loop iter (uniformR (34,633) mwc :: IO Word64) + ] + , bgroup "D" + [ bench "standard" $ whnfIO $ loop iter (standard mwc :: IO Double) + , bench "normal" $ whnfIO $ loop iter (normal 1 3 mwc :: IO Double) + , bench "exponential" $ whnfIO $ loop iter (exponential 3 mwc :: IO Double) + , bench "gamma,a<1" $ whnfIO $ loop iter (gamma 0.5 1 mwc :: IO Double) + , bench "gamma,a>1" $ whnfIO $ loop iter (gamma 2 1 mwc :: IO Double) + , bench "beta" $ whnfIO $ loop iter (beta 2 3 mwc :: IO Double) + , bench "chiSquare" $ whnfIO $ loop iter (chiSquare 4 mwc :: IO Double) + -- NOTE: We switch between algorithms when Np=10 + , bgroup "binomial" + [ bench (show p ++ " " ++ show n) $ whnfIO $ loop iter (binomial n p mwc) + | (n,p) <- [ (2, 0.2), (2, 0.5), (2, 0.8) + , (10, 0.1), (10, 0.9) + , (20, 0.2), (20, 0.8) + -- + , (60, 0.2), (60, 0.8) + , (600, 0.2), (600, 0.8) + , (6000, 0.2), (6000, 0.8) + ] + ] + ] + , bgroup "poisson" + [ bench (show lam) $ whnfIO $ loop iter (poisson lam mwc) + | lam <- [0.1, 1, 8, 12, 100, 1000, 10000] + ] + -- Test sampling performance. Table creation must be floated out! + , bgroup "CT/gen" $ concat + [ [ bench ("uniform "++show i) $ whnfIO $ loop iter (genFromTable tbl mwc) + | i <- [2..10] + , let tbl = makeTableUniform i + ] + , [ bench ("poisson " ++ show l) $ whnfIO $ loop iter (genFromTable tbl mwc) + | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000] + , let tbl = tablePoisson l + ] + , [ bench ("binomial " ++ show p ++ " " ++ show n) $ whnfIO $ loop iter (genFromTable tbl mwc) + | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)] + , let tbl = tableBinomial n p + ] + ] + -- Benchmarking of setting up table (no need to use iterations + -- here!). Setting up is rather expensive + , 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" + [ + bench "Double" $ whnfIO $ loop iter (R.randomIO :: IO Double) + , bench "Int" $ whnfIO $ loop iter (R.randomIO :: IO Int) + ] + , bgroup "mersenne" + [ + bench "Double" $ whnfIO $ loop iter (M.random mtg :: IO Double) + , bench "Int" $ whnfIO $ loop iter (M.random mtg :: IO Int) + ] +#if MIN_VERSION_random(1,3,0) + , bgroup "seed" + [ bench "SeedGen.fromSeed" $ let rseed = R.toSeed seed :: R.Seed Seed + in whnf R.fromSeed rseed + , bench "SeedGen.toSeed" $ whnf R.toSeed seed + ] +#endif + ] diff --git a/benchmarks/Benchmark.hs b/benchmarks/Benchmark.hs deleted file mode 100644 --- a/benchmarks/Benchmark.hs +++ /dev/null @@ -1,39 +0,0 @@ -import Control.Exception -import Control.Monad.ST -import Criterion.Main -import Data.Int -import Data.Word -import qualified System.Random as R -import System.Random.MWC -import qualified System.Random.Mersenne as M - -main = do - mwc <- create - 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) - ] - , bgroup "random" - [ - bench "Double" (R.randomIO >>= evaluate :: IO Double) - , bench "Int" (R.randomIO >>= evaluate :: IO Int) - ] - , bgroup "mersenne" - [ - bench "Double" (M.random mtg :: IO Double) - , bench "Int" (M.random mtg :: IO Int) - ] - ] diff --git a/benchmarks/Quickie.hs b/benchmarks/Quickie.hs deleted file mode 100644 --- a/benchmarks/Quickie.hs +++ /dev/null @@ -1,13 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -import System.Random.MWC (create, uniform) -import Control.Monad.ST (ST, runST) - -u :: ST s Double -u = do - let last = 1000000 :: Int - gen <- create - let loop !n !i | n == last = return i - | otherwise = uniform gen >>= loop (n+1) - loop 0 0 - -main = print (runST u) diff --git a/benchmarks/mwc-random-benchmarks.cabal b/benchmarks/mwc-random-benchmarks.cabal deleted file mode 100644 --- a/benchmarks/mwc-random-benchmarks.cabal +++ /dev/null @@ -1,18 +0,0 @@ -name: mwc-random-benchmarks -version: 0 -synopsis: Benchmarks for the mwc-random package -description: Benchmarks for the mwc-random package -license: BSD3 -license-file: ../LICENSE -build-type: Simple -cabal-version: >= 1.6 - -executable bm - main-is: Benchmark.hs - - build-depends: - base < 5, - criterion, - mersenne-random, - mwc-random, - random diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,128 @@ +## Changes in 0.15.3.0 [2025-12-30] + + * `poisson` samples for Poisson distribution is added + +## Changes in 0.15.2.0 [2025-01-18] + + * Support for `random-1.3`. + + +## Changes in 0.15.1.0 [2024-07-09] + + * Additon of binomial sampler using the rejection sampling method in + Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random + Variate Generation. Communications of the ACM, 31, 2 (February, + 1988) 216. . A more + efficient basis for e.g. the beta binomial distribution: + `beta a b g >>= \p -> binomial n p g`. + +## Changes in 0.15.0.2 [2021-08-15] + + * Doctests on 32-bit platforms are fixed. (#79) + + +## Changes in 0.15.0.1 [2020-08-08] + + * Bug in generation of Int/Word in both uniform and uniformR is fixed. (#75) + + +## Changes in 0.15.0.0 [2020-07-31] + + * `withSystemRandomST` and `createSystemSeed` are added. + + * `withSystemRandom` is deprecated. + + * `random>=1.2` is dependency of `mwc-random`. + + * Instances for type classes `StatefulGen` & `FrozenGen` defined in random-1.2 + are added for `Gen`. + + * Functions in `System.Random.MWC.Distributions` and + `System.Random.MWC.CondensedTable` now work with arbitrary `StatefulGen` + + * `System.Random.MWC.uniformVector` now works with arbitrary `StatefulGen` as + well and uses in-place initialization instead of `generateM`. It should be + faster for anything but IO and ST (those shoud remain same). + + +## Changes in 0.14.0.0 + + * Low level functions for acquiring random data for initialization + of PRGN state is moved to `System.Random.MWC.SeedSource` module + + * Ensure that carry is always correct when restoring PRNG state from + seed. Only affects users who create 258 element seed manually. + (#63, #65) + + +## Changes in 0.13.6.0 + + * `tablePoisson` now can handle λ>1923, see #59 for details. + That required intoduction of dependency on math-functions. + + +## Changes in 0.13.5.0 + + * `logCategorical` added + +## Changes in 0.13.4.0 + + * `withSystemRandom` uses RtlGenRandom for seeding generator on windows + + +## Changes in 0.13.3.1 + + * primitive-0.6 compatibility + + +## Changes in 0.13.3.0 + + * Monadic variant of vector shuffle added: `uniformShuffleM` + + * Context on `uniformShuffle` loosened + + +## Changes in 0.13.2.2 + + * Fixed crash during gen. initialization on Windows when stderr + is not available (#36). + +## Changes in 0.13.2.0 + + * Generators for beta, Bernoully, Dirichlet and categorical distributions + added. + + * Functions for generating random shuffles added. + + +## Changes in 0.13.1.2 + + * GHC 7.9 support + + +## Changes in 0.13.1.1 + + * Long standing performance problem in normal distribution fixed (#16) + + +## Changes in 0.13.1.0 + + * `createSystemRandom` added + + +## Changes in 0.13.0.0 + + * Workaround for GHC bug 8072 (bug 25). GHC 7.6 on 32-bit platrofms is + affected. + + * Generators for truncated exponential and geometric distributions + added. + + +## Changes in 0.12.0.0 + + * Fucntion `asGenIO` and `asGenST` added. + + * Generation of discrete random variates using condensed tables + methed. Tables for Poisson and binomial distributions are + provided. diff --git a/docs/RecreateFigure.svg b/docs/RecreateFigure.svg new file mode 100644 --- /dev/null +++ b/docs/RecreateFigure.svg @@ -0,0 +1,3 @@ + +0.51.01.568.6231.8859.540.5FourThreeTwoOnep = 0.25n = 200Target PDFMinorizing FunctionMajorizing FunctionProbability Density (Unnormalised)Number of Successes diff --git a/mwc-random.cabal b/mwc-random.cabal --- a/mwc-random.cabal +++ b/mwc-random.cabal @@ -1,5 +1,17 @@ +cabal-version: 3.0 +build-type: Simple name: mwc-random -version: 0.8.0.5 +version: 0.15.3.0 +license: BSD-2-Clause +license-file: LICENSE +copyright: 2009, 2010, 2011 Bryan O'Sullivan + +author: Bryan O'Sullivan +maintainer: Alexey Khudyakov +homepage: https://github.com/haskell/mwc-random +bug-reports: https://github.com/haskell/mwc-random/issues + +category: Math, Statistics synopsis: Fast, high quality pseudo random number generation description: This package contains code for generating high quality random @@ -14,45 +26,125 @@ Compared to the mersenne-random package, this package has a more convenient API, is faster, and supports more statistical distributions. -license: BSD3 -license-file: LICENSE -homepage: https://bitbucket.org/bos/mwc-random -bug-reports: https://bitbucket.org/bos/mwc-random/issues -author: Bryan O'Sullivan -maintainer: Bryan O'Sullivan -copyright: 2009, 2010, 2011 Bryan O'Sullivan -category: Math, Statistics -build-type: Simple -cabal-version: >= 1.6 -extra-source-files: - README.markdown - benchmarks/*.hs - benchmarks/Quickie.hs - benchmarks/mwc-random-benchmarks.cabal -library - exposed-modules: - System.Random.MWC - build-depends: - base < 5, - primitive, - time, - vector >= 0.5 - if impl(ghc >= 6.10) - build-depends: - base >= 4 - -- gather extensive profiling data for now - ghc-prof-options: -auto-all +extra-source-files: + README.md - ghc-options: -Wall -funbox-strict-fields - if impl(ghc >= 6.8) - ghc-options: -fwarn-tabs +extra-doc-files: + docs/*.svg + changelog.md -source-repository head - type: mercurial - location: https://bitbucket.org/bos/mwc-random +tested-with: + GHC ==8.0.2 + || ==8.2.2 + || ==8.4.4 + || ==8.6.5 + || ==8.8.4 + || ==8.10.7 + || ==9.0.2 + || ==9.2.8 + || ==9.4.8 + || ==9.6.7 + || ==9.8.4 + || ==9.10.2 + || ==9.12.2 + source-repository head type: git - location: git://github.com/bos/mwc-random + location: https://github.com/haskell/mwc-random.git + +flag BenchPAPI + Description: Enable building of benchmarks which use instruction counters. + It requires libpapi and only works on Linux so it's protected by flag + Default: False + Manual: True + +library + default-language: Haskell2010 + exposed-modules: System.Random.MWC + System.Random.MWC.Distributions + System.Random.MWC.CondensedTable + System.Random.MWC.SeedSource + build-depends: base >= 4.9 && < 5 + , primitive >= 0.6.2 + , random >= 1.2 + , time + , vector >= 0.7 + , math-functions >= 0.2.1.0 + + ghc-options: -O2 -Wall -funbox-strict-fields -fwarn-tabs + + +-- We want to be able to build benchmarks using both tasty-bench and tasty-papi. +-- They have similar API so we just create two shim modules which reexport +-- definitions from corresponding library and pick one in cabal file. +common bench-stanza + ghc-options: -Wall + default-language: Haskell2010 + build-depends: base < 5 + , vector >= 0.11 + , mersenne-random + , mwc-random + , random + , tasty >=1.3.1 + +benchmark mwc-bench + import: bench-stanza + type: exitcode-stdio-1.0 + hs-source-dirs: bench bench-time + main-is: Benchmark.hs + Other-modules: Bench + build-depends: tasty-bench >= 0.3 + +benchmark mwc-bench-papi + import: bench-stanza + type: exitcode-stdio-1.0 + if impl(ghcjs) || !flag(BenchPAPI) || impl(ghc < 8.2) + buildable: False + hs-source-dirs: bench bench-papi + main-is: Benchmark.hs + Other-modules: Bench + build-depends: tasty-papi >= 0.1.2 + +test-suite mwc-prop-tests + type: exitcode-stdio-1.0 + hs-source-dirs: tests + main-is: props.hs + default-language: Haskell2010 + ghc-options: + -Wall -threaded -rtsopts + + build-depends: base + , mwc-random + , QuickCheck >=2.2 + , vector >=0.12.1 + , tasty >=1.3.1 + , tasty-quickcheck >=0.10.2 + , tasty-hunit + , random >=1.2 + , mtl + , math-functions >=0.3.4 + +test-suite mwc-doctests + type: exitcode-stdio-1.0 + main-is: doctests.hs + hs-source-dirs: tests + default-language: Haskell2010 + if impl(ghcjs) || impl(ghc < 8.0) + Buildable: False + -- Linker on macos prints warnings to console which confuses doctests. + -- We simply disable doctests on ma for older GHC + -- > warning: -single_module is obsolete + if os(darwin) && impl(ghc < 9.6) + buildable: False + build-depends: + base -any + , mwc-random -any + , doctest >=0.15 && <0.24 + -- + , bytestring + , primitive + , vector >=0.11 + , random >=1.2 diff --git a/tests/doctests.hs b/tests/doctests.hs new file mode 100644 --- /dev/null +++ b/tests/doctests.hs @@ -0,0 +1,4 @@ +import Test.DocTest (doctest) + +main :: IO () +main = doctest ["-fobject-code", "System/Random/MWC.hs"] diff --git a/tests/props.hs b/tests/props.hs new file mode 100644 --- /dev/null +++ b/tests/props.hs @@ -0,0 +1,316 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ViewPatterns #-} +import Control.Monad +import Data.Word +import Data.Proxy +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Mutable as MVU +import Numeric.SpecFunctions (logChoose,incompleteGamma,log1p,logFactorial) + +import Test.Tasty +import Test.Tasty.QuickCheck +import Test.Tasty.Runners +import Test.Tasty.Options +import Test.Tasty.HUnit +import Test.QuickCheck.Monadic + +import System.Random.MWC +import System.Random.MWC.Distributions +import System.Random.Stateful (StatefulGen) +#if MIN_VERSION_random(1,3,0) +import qualified System.Random.Stateful as Random (SeedGen(..)) +#endif + +---------------------------------------------------------------- +-- +---------------------------------------------------------------- + +-- | Average number of events per bin for binned statistical tests +newtype NPerBin = NPerBin Int + +instance IsOption NPerBin where + defaultValue = NPerBin 100 + parseValue = fmap NPerBin . safeRead + optionName = pure "n-per-bin" + optionHelp = pure "Average number of events per bin" + + +-- | P-value for statistical test. +newtype PValue = PValue Double + +instance IsOption PValue where + defaultValue = PValue 1e-9 + parseValue = fmap PValue . safeRead + optionName = pure "pvalue" + optionHelp = pure "P-value for statistical test" + +---------------------------------------------------------------- +-- +---------------------------------------------------------------- + +main :: IO () +main = do + -- Set up tasty + let tasty_opts = [ Option (Proxy :: Proxy NPerBin) + , Option (Proxy :: Proxy PValue) + , Option (Proxy :: Proxy QuickCheckTests) + , Option (Proxy :: Proxy QuickCheckReplay) + , Option (Proxy :: Proxy QuickCheckShowReplay) + , Option (Proxy :: Proxy QuickCheckMaxSize) + , Option (Proxy :: Proxy QuickCheckMaxRatio) + , Option (Proxy :: Proxy QuickCheckVerbose) + , Option (Proxy :: Proxy QuickCheckMaxShrinks) + ] + ingredients = includingOptions tasty_opts : defaultIngredients + opts <- parseOptions ingredients (testCase "Fake" (pure ())) + let n_per_bin = lookupOption opts :: NPerBin + p_val = lookupOption opts + -- + g0 <- createSystemRandom + defaultMainWithIngredients ingredients $ testGroup "mwc" + [ testProperty "save/restore" $ prop_SeedSaveRestore g0 +#if MIN_VERSION_random(1,3,0) + , testProperty "SeedGen" $ prop_SeedGen g0 +#endif + , testCase "user save/restore" $ saveRestoreUserSeed + , testCase "empty seed data" $ emptySeed + , testCase "output correct" $ do + g <- create + xs <- replicateM 513 (uniform g) + assertEqual "[Word32]" xs golden + , testCase "beta binomial mean" $ prop_betaBinomialMean + , testProperty "binomial is binomial" $ prop_binomial_PMF n_per_bin p_val g0 + , testProperty "Poisson is Poisson" $ prop_Poisson_PMF n_per_bin p_val g0 + , testCase "poisson mean and variance of 1" $ prop_poissonMeanAndVar 1 + , testCase "poisson mean and variance of 5" $ prop_poissonMeanAndVar 5 + , testCase "poisson mean and variance of 40" $ prop_poissonMeanAndVar 40 + , testCase "poisson mean and variance of 150" $ prop_poissonMeanAndVar 150 + , testCase "poisson 0" $ do + n <- poisson 0 g0 + 0 @=? n + ] + +updateGenState :: GenIO -> IO () +updateGenState g = replicateM_ 250 (uniform g :: IO Word32) + +prop_SeedSaveRestore :: GenIO -> Property +prop_SeedSaveRestore g = monadicIO $ do + run $ updateGenState g + seed <- run $ save g + seed' <- run $ save =<< restore seed + return $ seed == seed' + +#if MIN_VERSION_random(1,3,0) +prop_SeedGen :: GenIO -> Property +prop_SeedGen g = monadicIO $ do + run $ updateGenState g + seed <- run $ save g + return $ seed == (Random.fromSeed . Random.toSeed) seed +#endif + +saveRestoreUserSeed :: IO () +saveRestoreUserSeed = do + let seed = toSeed $ U.replicate 258 0 + seed' <- save =<< restore seed + assertEqual "Seeds must be equal" seed' seed + +emptySeed :: IO () +emptySeed = do + let seed = toSeed U.empty + seed' <- save =<< create + assertEqual "Seeds must be equal" seed' seed + +-- First 513 values generated from seed made using create +golden :: [Word32] +golden = + [ 2254043345, 562229898, 1034503294, 2470032534, 2831944869, 3042560015, 838672965, 715056843 + , 3122641307, 2300516242, 4079538318, 3722020688, 98524204, 1450170923, 2669500465, 2890402829 + , 114212910, 1914313000, 2389251496, 116282477, 1771812561, 1606473512, 1086420617, 3652430775 + , 1165083752, 3599954795, 3006722175, 341614641, 3000394300, 1378097585, 1551512487, 81211762 + , 604209599, 3949866361, 77745071, 3170410267, 752447516, 1213023833, 1624321744, 3251868348 + , 1584957570, 2296897736, 3305840056, 1158966242, 2458014362, 1919777052, 3203159823, 3230279656 + , 755741068, 3005087942, 2478156967, 410224731, 1196248614, 3302310440, 3295868805, 108051054 + , 1010042411, 2725695484, 2201528637, 667561409, 79601486, 50029770, 566202616, 3217300833 + , 2162817014, 925506837, 1527015413, 3079491438, 927252446, 118306579, 499811870, 2973454232 + , 2979271640, 4078978924, 1864075883, 197741457, 296365782, 1784247291, 236572186, 464208268 + , 1769568958, 827682258, 4247376295, 2959098022, 1183860331, 2475064236, 3952901213, 1953014945 + , 393081236, 1616500498, 2201176136, 1663813362, 2167124739, 630903810, 113470040, 924745892 + , 1081531735, 4039388931, 4118728223, 107819176, 2212875141, 1941653033, 3660517172, 192973521 + , 3653156164, 1878601439, 3028195526, 2545631291, 3882334975, 456082861, 2775938704, 3813508885 + , 1758481462, 3332769695, 3595846251, 3745876981, 152488869, 2555728588, 3058747945, 39382408 + , 520595021, 2185388418, 3502636573, 2650173199, 1077668433, 3548643646, 71562049, 2726649517 + , 494210825, 1208915815, 620990806, 2877290965, 3253243521, 804166732, 2481889113, 623399529 + , 44880343, 183645859, 3283683418, 2214754452, 419328482, 4224066437, 1102669380, 1997964721 + , 2437245376, 985749802, 858381069, 116806511, 1771295365, 97352549, 341972923, 2971905841 + , 110707773, 950500868, 1237119233, 691764764, 896381812, 1528998276, 1269357470, 2567094423 + , 52141189, 2722993417, 80628658, 3919817965, 3615946076, 899371181, 46940285, 4010779728 + , 318101834, 30736609, 3577200709, 971882724, 1478800972, 3769640027, 3706909300, 3300631811 + , 4057825972, 4285058790, 2329759553, 2967563409, 4080096760, 2762613004, 2518395275, 295718526 + , 598435593, 2385852565, 2608425408, 604857293, 2246982455, 919156819, 1721573814, 2502545603 + , 643962859, 587823425, 3508582012, 1777595823, 4119929334, 2833342174, 414044876, 2469473258 + , 289159600, 3715175415, 966867024, 788102818, 3197534326, 3571396978, 3508903890, 570753009 + , 4273926277, 3301521986, 1411959102, 2766249515, 4071012597, 959442028, 1962463990, 1098904190 + , 714719899, 562204808, 1658783410, 1471669042, 2565780129, 1616648894, 4236521717, 1788863789 + , 3068674883, 191936470, 253084644, 1915647866, 276372665, 2117183118, 3704675319, 218791054 + , 3680045802, 406662689, 3844864229, 91140313, 3834015630, 25116147, 904830493, 3152559113 + , 820358622, 1301896358, 296152699, 2202014455, 4256659428, 1175171414, 3287520873, 2028006499 + , 327448717, 2095642873, 3798661296, 58567008, 3907537112, 3691259011, 1730142328, 2373011713 + , 3387040741, 3189417655, 2949233059, 1238379614, 1813238023, 1064726446, 1339055235, 1744523609 + , 279811576, 2934103599, 283542302, 994488448, 418691747, 1062780152, 102211875, 4071713296 + , 1790834038, 1035092527, 2374272359, 3558280982, 1927663822, 3645417844, 3481790745, 3566282546 + , 2000290859, 505518126, 363501589, 4075468679, 3247300709, 3705242654, 2731103609, 2836871038 + , 589640144, 2546495106, 84767518, 1376911639, 2400770705, 527489676, 3804134352, 150084021 + , 240070593, 3807594859, 3518576690, 659503830, 2239678479, 1273668921, 4271050554, 3090482972 + , 401956859, 1772128561, 4438455, 1989666158, 2521484677, 3960178700, 4220196277, 1033999035 + , 2214785840, 3428469341, 428564336, 2517446784, 3935757188, 3294001677, 1037971963, 3590324170 + , 1220969729, 1719719817, 807688972, 77076422, 4251553858, 3963852375, 326128795, 3277818295 + , 3671513069, 549617771, 1683950556, 3352913781, 409318429, 2456264774, 4036950639, 1162718475 + , 83888874, 5578966, 172866494, 1542278848, 455546979, 1296511553, 4263636440, 2450589064 + , 372411483, 211216338, 2632256495, 2393754408, 1336054289, 4087203071, 3159642437, 1933346856 + , 2914152714, 3805541979, 2769740793, 1161287028, 2289749561, 4124509890, 2128452935, 210531695 + , 4250709834, 390950534, 1421430300, 3030519715, 3228987297, 3086837053, 2866915453, 2335948692 + , 1684378991, 2575634059, 4153427304, 2426048796, 4197556954, 2605152326, 2909410733, 2424889219 + , 654577921, 811955499, 118126602, 504071559, 1278756230, 3896458168, 4105558075, 750276169 + , 1120805572, 1762689330, 993728154, 1104363215, 774344996, 4077568952, 2183487324, 994724370 + , 3323036885, 3880704963, 746305447, 961608310, 2030117337, 453935768, 800490463, 1034636 + , 2323633564, 602565693, 806061242, 1899269713, 162686347, 467541008, 1529175313, 282891502 + , 2529616339, 2930657178, 464272784, 2878535316, 807165854, 3209080518, 4080120278, 347748171 + , 3972126063, 284174728, 2498328933, 1723872460, 143845955, 4223866687, 1761495357, 1544646770 + , 4206103283, 3771574626, 642165282, 1119501013, 3514063332, 1443320304, 4056369796, 3602131475 + , 1422908288, 804093687, 431176780, 40108717, 2998264213, 3705835674, 169805085, 454593842 + , 2781536994, 2385225212, 4137367775, 2631435125, 2347082354, 629238010, 3283635219, 3815791831 + , 1340400558, 4061846985, 3803921868, 3196119096, 718610843, 3694290834, 2169960411, 2407155570 + , 2557480499, 16164105, 480957288, 2155919829, 2490067282, 2356287132, 511737296, 1602800634 + , 1802275249, 3316832299, 50286484, 2106622541, 2352302834, 2538374315, 344766394, 2777260569 + , 1215135803, 2229011963, 114632277, 1645499402, 1111617833, 3833259754, 928611385, 686744723 + , 1898396834, 2461932251, 2665457318, 3797019621, 868313114, 2366635205, 481934875, 1170532970 + , 642610859, 3150733309, 3508548582, 666714469, 711663449, 2436617656, 2681476315, 1637296693 + , 2487349478, 4174144946, 2793869557, 559398604, 1898140528, 991962870, 864792875, 3861665129 + , 4024051364, 3383200293, 773730975, 33517291, 2660126073, 689133464, 2248134097, 3874737781 + , 3358012678] + +-- We can test two for the price of one +betaBinomial :: StatefulGen g m => + Double -> Double -> Int -> g -> m Int +betaBinomial a b n g = do + p <- beta a b g + binomial n p g + +nSamples :: Int +nSamples = 10000 + +alpha, delta :: Double +alpha = 600.0 +delta = 400.0 + +nTrials :: Int +nTrials = 10 + +prop_betaBinomialMean :: IO () +prop_betaBinomialMean = do + g <- create + ss <- replicateM nSamples $ betaBinomial alpha delta nTrials g + let m = fromIntegral (sum ss) / fromIntegral nSamples + let x1 = fromIntegral nTrials * alpha / (alpha + delta) + assertBool ("Mean is " ++ show x1 ++ " but estimated as " ++ show m) (abs (m - x1) < 0.001) + +prop_poissonMeanAndVar :: Double -> IO () +prop_poissonMeanAndVar lambda = do + gen <- create + ss <- replicateM nSamples $ poisson lambda gen + let m = fromIntegral (sum ss) / fromIntegral nSamples :: Double + let v = (fromIntegral (sum (map (^ 2) ss)) / fromIntegral nSamples) - (m ** 2) + assertBool + ("True mean and var: " ++ show lambda ++ " but estimated as " ++ show (m, v) ++ " respectively") + ((abs (lambda - m) / lambda < 0.025) && (abs (lambda - v) / lambda < 0.025)) + + +-- Test that `binomial` really samples from binomial distribution. +-- +-- If we have binomial random variate with number of trials N and +-- sample it M times. Then number of events with K successes is +-- described by multinomial distribution and we can test whether +-- experimental distribution is described using likelihood ratio test +prop_binomial_PMF :: NPerBin -> PValue -> GenIO -> Property +prop_binomial_PMF (NPerBin n_per_bin) (PValue p_val) g = property $ do + p <- choose (0, 1.0) -- Success probability + n_trial <- choose (2, 100) -- Number of trials in binomial distribution + -- Number of binomial samples to generate + let n_samples = n_trial * n_per_bin + n_samples' = fromIntegral n_samples + -- Compute number of outcomes + pure $ ioProperty $ do + hist <- do + buf <- MVU.new (n_trial + 1) + replicateM_ n_samples $ + MVU.modify buf (+(1::Int)) =<< binomial n_trial p g + U.unsafeFreeze buf + -- Here we compute twice log of likelihood ratio. Alternative + -- hypothesis is some distribution which fits data perfectly + -- + -- Asymtotically it's ditributed as χ² with n_trial-1 degrees of + -- freedom + let likelihood _ 0 + = 0 + likelihood k (fromIntegral -> n_obs) + = n_obs * (log (n_obs / n_samples') - logProbBinomial n_trial p k) + let logL = 2 * U.sum (U.imap likelihood hist) + let significance = 1 - cumulativeChi2 (n_trial - 1) logL + pure $ counterexample ("p = " ++ show p) + $ counterexample ("N = " ++ show n_trial) + $ counterexample ("p-val = " ++ show significance) + $ counterexample ("chi2 = " ++ show logL) + $ significance > p_val + +-- Similar test for Poisson distribution. We build histogram by moving +-- all values >λ+2σ into overflow bin. +prop_Poisson_PMF :: NPerBin -> PValue -> GenIO -> Property +prop_Poisson_PMF (NPerBin n_per_bin) (PValue p_val) g = property $ do + lam <- choose (0, 100.0) -- Poisson rate + let max_n = ceiling $ lam + 2 * sqrt lam -- Maximum histogram bin + n_samples = max_n * n_per_bin -- Number of samples to generate + n_samples' = fromIntegral n_samples + pure $ ioProperty $ do + hist <- do + buf <- MVU.new (max_n + 1) + replicateM_ n_samples $ do + k <- poisson lam g + MVU.modify buf (+(1::Int)) (min max_n k) + U.unsafeFreeze buf + let likelihood _ 0 + = 0 + likelihood k (fromIntegral -> n_obs) + = n_obs * (log (n_obs / n_samples') - p) + where + p | k == max_n = logComlCdfPoisson lam k + | otherwise = logProbPoisson lam k + let logL = 2 * U.sum (U.imap likelihood hist) + let significance = 1 - cumulativeChi2 max_n logL + pure $ counterexample ("lambda = " ++ show lam) + $ counterexample ("p-val = " ++ show significance) + $ counterexample ("chi2 = " ++ show logL) + $ significance > p_val + +---------------------------------------------------------------- +-- Statistical helpers +---------------------------------------------------------------- + +-- Logarithm of probability for binomial distribution +logProbBinomial :: Int -> Double -> Int -> Double +logProbBinomial n p k + = logChoose n k + log p * k' + log1p (-p) * nk' + where + k' = fromIntegral k + nk' = fromIntegral $ n - k + +logProbPoisson :: Double -> Int -> Double +logProbPoisson lam i = log lam * fromIntegral i - logFactorial i - lam + +logComlCdfPoisson :: Double -> Int -> Double +logComlCdfPoisson lam i = log $ incompleteGamma (fromIntegral (i :: Int)) lam + + +cumulativeChi2 :: Int -> Double -> Double +cumulativeChi2 (fromIntegral -> ndf) x + | x <= 0 = 0 + | otherwise = incompleteGamma (ndf/2) (x/2)