diff --git a/System/Random/MWC.hs b/System/Random/MWC.hs --- a/System/Random/MWC.hs +++ b/System/Random/MWC.hs @@ -96,7 +96,6 @@ #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) @@ -110,13 +109,13 @@ import Data.Word (Word, Word8, Word16, Word32, Word64) import Foreign.Marshal.Alloc (allocaBytes) import Foreign.Marshal.Array (peekArray) -import Prelude hiding (catch) import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as I import qualified Data.Vector.Unboxed.Mutable as M import System.CPUTime (cpuTimePrecision, getCPUTime) import System.IO (IOMode(..), hGetBuf, hPutStrLn, stderr, withBinaryFile) import System.IO.Unsafe (unsafePerformIO) +import qualified Control.Exception as E @@ -292,7 +291,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 @@ -370,11 +369,11 @@ | otherwise = G.unsafeIndex seed i fini = G.length seed {-# INLINE initialize #-} - + -- | An immutable snapshot of the state of a 'Gen'. -newtype Seed = Seed { +newtype Seed = Seed { -- | Convert seed into vector. - fromSeed :: I.Vector Word32 + fromSeed :: I.Vector Word32 } deriving (Eq, Show, Typeable) @@ -415,7 +414,7 @@ nread <- withBinaryFile random ReadMode $ \h -> hGetBuf h buf nbytes peekArray (nread `div` 4) buf - + -- | Seed a PRNG with data from the system's fast source of -- pseudo-random numbers (\"@\/dev\/urandom@\" on Unix-like systems), -- then run the given action. @@ -430,7 +429,7 @@ -- highly independent. withSystemRandom :: PrimMonad m => (Gen (PrimState m) -> m a) -> IO a withSystemRandom act = do - seed <- acquireSeedSystem `catch` \(_::IOException) -> do + seed <- acquireSeedSystem `E.catch` \(_::E.IOException) -> do seen <- atomicModifyIORef warned ((,) True) unless seen $ do hPutStrLn stderr ("Warning: Couldn't open /dev/urandom") diff --git a/System/Random/MWC/CondensedTable.hs b/System/Random/MWC/CondensedTable.hs --- a/System/Random/MWC/CondensedTable.hs +++ b/System/Random/MWC/CondensedTable.hs @@ -102,20 +102,29 @@ -- 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) + :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32, Show a) => v (a, Double) -> CondensedTable v a {-# INLINE tableFromProbabilities #-} tableFromProbabilities v - | G.null v = pkgError "tableFromProbabilities" "empty vector of outcomes" - | otherwise = tableFromIntWeights $ G.map (second $ round . (* mlt)) v + | G.null tbl = pkgError "tableFromProbabilities" "empty vector of outcomes" + | otherwise = tableFromIntWeights $ G.map (second $ toWeight . (* mlt)) tbl where - mlt = 4.294967296e9 -- 2^32 + -- 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) + :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32, Show a) => v (a, Double) -> CondensedTable v a {-# INLINE tableFromWeights #-} tableFromWeights = tableFromProbabilities . normalize . G.filter ((> 0) . snd) @@ -131,12 +140,12 @@ -- | Generate a condensed lookup table from integer weights. Weights -- should sum to @2^32@. If they don't, the algorithm will alter the -- weights so that they do. This approach should work reasonably well --- for rounding error. +-- for rounding errors. tableFromIntWeights :: (Vector v (a,Word32), Vector v a, Vector v Word32) => v (a, Word32) -> CondensedTable v a {-# INLINE tableFromIntWeights #-} -tableFromIntWeights tbl +tableFromIntWeights v | n == 0 = pkgError "tableFromIntWeights" "empty table" -- Single element tables should be treated sepately. Otherwise -- they will confuse correctWeights @@ -152,6 +161,9 @@ 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 diff --git a/mwc-random.cabal b/mwc-random.cabal --- a/mwc-random.cabal +++ b/mwc-random.cabal @@ -1,5 +1,5 @@ name: mwc-random -version: 0.12.0.0 +version: 0.12.0.1 synopsis: Fast, high quality pseudo random number generation description: This package contains code for generating high quality random @@ -61,12 +61,12 @@ main-is: tests.hs other-modules: KS QC - Uniform ghc-options: -Wall -threaded -rtsopts build-depends: + vector >= 0.7, HUnit, QuickCheck, base, diff --git a/test/KS.hs b/test/KS.hs new file mode 100644 --- /dev/null +++ b/test/KS.hs @@ -0,0 +1,48 @@ +-- Kolmogorov-Smirnov tests for distribution +-- +-- Note that it's not most powerful test for normality. +module KS ( + tests + ) where + +import qualified Data.Vector.Unboxed as U + +import Statistics.Test.KolmogorovSmirnov + +import Statistics.Distribution +import Statistics.Distribution.Exponential +import Statistics.Distribution.Gamma +import Statistics.Distribution.Normal +import Statistics.Distribution.Uniform + +import qualified System.Random.MWC as MWC +import qualified System.Random.MWC.Distributions as MWC + +import Test.HUnit hiding (Test) +import Test.Framework +import Test.Framework.Providers.HUnit + + +tests :: MWC.GenIO -> Test +tests g = testGroup "Kolmogorov-Smirnov" + [ testCase "standard" $ testKS standard MWC.standard g + , testCase "normal m=1 s=2" $ testKS (normalDistr 1 2) (MWC.normal 1 2) g + -- Gamma distribution + , testCase "gamma k=1 θ=1" $ testKS (gammaDistr 1 1 ) (MWC.gamma 1 1 ) g + , testCase "gamma k=0.3 θ=0.4" $ testKS (gammaDistr 0.3 0.4) (MWC.gamma 0.3 0.4) g + , testCase "gamma k=0.3 θ=3" $ testKS (gammaDistr 0.3 3 ) (MWC.gamma 0.3 3 ) g + , testCase "gamma k=3 θ=0.4" $ testKS (gammaDistr 3 0.4) (MWC.gamma 3 0.4) g + , testCase "gamma k=3 θ=3" $ testKS (gammaDistr 3 3 ) (MWC.gamma 3 3 ) g + -- Uniform + , testCase "uniform -2 .. 3" $ testKS (uniformDistr (-2) 3) (MWC.uniformR (-2,3)) g + -- Exponential + , testCase "exponential l=1" $ testKS (exponential 1) (MWC.exponential 1) g + , testCase "exponential l=3" $ testKS (exponential 3) (MWC.exponential 3) g + ] + +testKS :: (Distribution d) => d -> (MWC.GenIO -> IO Double) -> MWC.GenIO -> IO () +testKS distr generator g = do + sample <- U.replicateM 1000 (generator g) + case kolmogorovSmirnovTest distr 0.01 sample of + Significant -> assertFailure "KS test failed" + NotSignificant -> return () diff --git a/test/QC.hs b/test/QC.hs new file mode 100644 --- /dev/null +++ b/test/QC.hs @@ -0,0 +1,56 @@ +-- QC tests for random number generators +-- +-- Require QuickCheck >= 2.2 +module QC ( + tests + ) where + +import Control.Applicative + +import Data.Word (Word8,Word16,Word32,Word64,Word) +import Data.Int (Int8, Int16, Int32, Int64 ) + +import Test.QuickCheck +import Test.QuickCheck.Monadic +import Test.Framework +import Test.Framework.Providers.QuickCheck2 + +import System.Random.MWC + + + +---------------------------------------------------------------- + +tests :: GenIO -> Test +tests g = testGroup "Range" + [ testProperty "Int8" $ (prop_InRange g :: InRange Int8) + , testProperty "Int16" $ (prop_InRange g :: InRange Int16) + , testProperty "Int32" $ (prop_InRange g :: InRange Int32) + , testProperty "Int64" $ (prop_InRange g :: InRange Int64) + , testProperty "Word8" $ (prop_InRange g :: InRange Word8) + , testProperty "Word16" $ (prop_InRange g :: InRange Word16) + , testProperty "Word32" $ (prop_InRange g :: InRange Word32) + , testProperty "Word64" $ (prop_InRange g :: InRange Word64) + , testProperty "Int" $ (prop_InRange g :: InRange Int) + , testProperty "Word64" $ (prop_InRange g :: InRange Word) + , testProperty "Float" $ (prop_InRange g :: InRange Float) + , testProperty "Double" $ (prop_InRange g :: InRange Double) + ] + + + +---------------------------------------------------------------- + +-- Test that values generated with uniformR never lie outside range. +prop_InRange :: (Variate a, Ord a,Num a) => GenIO -> OrderedPair a -> Property +prop_InRange g (OrderedPair (x1,x2)) = monadicIO $ do + r <- run $ uniformR (x1,x2) g + assert (x1 <= r && r <= x2) + +type InRange a = OrderedPair a -> Property + +-- Ordered pair (x,y) for which x <= y +newtype OrderedPair a = OrderedPair (a,a) + deriving Show +instance (Ord a, Arbitrary a) => Arbitrary (OrderedPair a) where + arbitrary = OrderedPair <$> suchThat arbitrary (uncurry (<=)) diff --git a/test/tests.hs b/test/tests.hs new file mode 100644 --- /dev/null +++ b/test/tests.hs @@ -0,0 +1,16 @@ +import Test.Framework (defaultMain) +import System.Random.MWC (withSystemRandom) + +import qualified QC +import qualified ChiSquare +import qualified KS + + +main :: IO () +main = + withSystemRandom $ \g -> + defaultMain + [ QC.tests g + , ChiSquare.tests g + , KS.tests g + ]