mwc-random 0.12.0.0 → 0.12.0.1

raw patch · 6 files changed

+148/−17 lines, 6 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• System/Random/MWC.hs +7/−8
• System/Random/MWC/CondensedTable.hs +19/−7
• mwc-random.cabal +2/−2
• test/KS.hs +48/−0
• test/QC.hs +56/−0
• test/tests.hs +16/−0

@@ -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")

@@ -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

@@ -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,

@@ -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 ()

@@ -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 (<=))

@@ -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+    ]