diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -3,6 +3,8 @@ {-# LANGUAGE FlexibleInstances #-} module Main ( main ) where +import Data.Word + import Control.Eff import Control.Eff.State.Strict import Control.Applicative @@ -41,9 +43,28 @@ main = defaultMainWith config (return ()) [ bcompare [ - bench "fastIntDist" (fastIntDist 1 1000 :: Eff (State Random :> ()) Int) - , bench "uniformIntDist" (uniformIntDist 1 1000 :: Eff (State Random :> ()) Integer) + bench "fastIntDist" (fastIntDist 1 1000 :: Eff (State Random :> ()) Int) + , bench "uniformIntDist" (uniformIntDist 1 1000 :: Eff (State Random :> ()) Integer) + , bench "uniformIntegralDist" (uniformIntegralDist 1 1000 :: Eff (State Random :> ()) Word64) ] + , bcompare [ + bench "normalDist" (normalDist 0 10 :: Eff (State Random :> ()) Double) + , bench "lognormalDist" (lognormalDist 0 10 :: Eff (State Random :> ()) Double) + ] + , bcompare [ + bench "randomWord32" (randomWord :: Eff (State Random :> ()) Word) + , bench "randomWord64" (randomWord64 :: Eff (State Random :> ()) Word64) + , bench "randomBits32" (randomBits :: Eff (State Random :> ()) Word32) + , bench "randomBits64" (randomBits :: Eff (State Random :> ()) Word64) + , bench "randomBitList64" (randomBitList 64 :: Eff (State Random :> ()) [Bool]) + , bench "randomBitList128" (randomBitList 128 :: Eff (State Random :> ()) [Bool]) + ] + , bcompare [ + bench "randomDouble" (randomDouble :: Eff (State Random :> ()) Double) + , bench "uniformRealDist" (uniformRealDist 0 10 :: Eff (State Random :> ()) Double) + , bench "linearRealDist" (linearRealDist 0 10 :: Eff (State Random :> ()) Double) + ] + , bench "binomialDist" (binomialDist 1000 0.5 :: Eff (State Random :> ()) Int) ] where config = defaultConfig { diff --git a/src/System/Random/Effect.hs b/src/System/Random/Effect.hs --- a/src/System/Random/Effect.hs +++ b/src/System/Random/Effect.hs @@ -3,6 +3,8 @@ {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} -- | A random number effect, using a pure mersenne twister under -- the hood. This should be plug-and-play with any application -- making use of extensible effects. @@ -13,12 +15,16 @@ -- * Seeding , mkRandom , mkRandomIO + , mkSecureRandomIO -- * Running + , forRandEff , runRandomState -- * Uniform Distributions , uniformIntDist , uniformIntegralDist , uniformRealDist + -- * Linear Distributions + , linearRealDist -- * Bernoulli Distributions , bernoulliDist , binomialDist @@ -42,6 +48,10 @@ , buildDDH , discreteDist , piecewiseConstantDist + , piecewiseLinearDist + -- * Shuffling + , knuthShuffle + , knuthShuffleM -- * Raw Generators , randomInt , randomInt64 @@ -54,15 +64,17 @@ import Control.Applicative import Control.Monad +import Control.Monad.Primitive import Control.Monad.ST import Data.Bits -import Data.Int import Data.List import Data.Ratio import Data.Typeable import Data.Vector.Algorithms.Search -import Data.Vector ( Vector ) +import Data.Vector ( Vector, (//), (!) ) import qualified Data.Vector as V +import Data.Vector.Mutable (MVector) +import qualified Data.Vector.Mutable as M import Data.Word import Statistics.Distribution @@ -75,79 +87,23 @@ --import qualified Statistics.Distribution.Poisson as DP import qualified Statistics.Distribution.StudentT as DS -import qualified System.Random.Mersenne.Pure64 as SR - import Control.Eff import Control.Eff.Lift import Control.Eff.State.Strict --- | A pure mersenne twister pseudo-random number generator. -data Random = Random {-# UNPACK #-} !SR.PureMT - deriving Typeable - --- | Create a random number generator from a 'Word64' seed. -mkRandom :: Word64 -> Random -mkRandom = Random . SR.pureMT -{-# INLINE mkRandom #-} - --- | Create a new random number generator, using the clocktime as the base for --- the seed. This must be called from a computation with a lifted base effect --- of 'IO'. -mkRandomIO :: SetMember Lift (Lift IO) r - => Eff r Random -mkRandomIO = lift (fmap Random SR.newPureMT) -{-# INLINE mkRandomIO #-} - --- | Runs an effectful random computation, returning the computation's result. -runRandomState :: Random - -> Eff (State Random :> r) w - -> Eff r w -runRandomState seed computation = - snd <$> runState seed computation -{-# INLINE runRandomState #-} - --- | A generalized form of generating a random number of the correct type --- from System.Random.Mersenne.Pure64. -randomF :: Member (State Random) r - => (SR.PureMT -> (a, SR.PureMT)) - -> Eff r a -randomF f = do - (Random old) <- get - let (val, new) = f old - put (Random new) - return val -{-# INLINE randomF #-} - --- | Yield a new 'Int' value from the generator. The full 64 bits will be used --- on a 64 bit machine. -randomInt :: Member (State Random) r => Eff r Int -randomInt = randomF SR.randomInt -{-# INLINE randomInt #-} - --- | Yield a new 'Word' value from the generator. -randomWord :: Member (State Random) r => Eff r Word -randomWord = randomF SR.randomWord -{-# INLINE randomWord #-} - --- | Yield a new 'Int64' value from the generator. -randomInt64 :: Member (State Random) r => Eff r Int64 -randomInt64 = randomF SR.randomInt64 -{-# INLINE randomInt64 #-} - --- | Yield a new 'Word64' value from the generator. -randomWord64 :: Member (State Random) r => Eff r Word64 -randomWord64 = randomF SR.randomWord64 -{-# INLINE randomWord64 #-} +import System.Random.Effect.Raw --- | Yield a new 53-bit precise 'Double' value from the generator. --- The returned number will be in the range [0, 1). -randomDouble :: Member (State Random) r => Eff r Double -randomDouble = randomF SR.randomDouble +randomDouble :: Member (State Random) r + => Eff r Double +randomDouble = do + r <- randomWord64 + return (fromIntegral (r `div` 2048) / 9007199254740992) {-# INLINE randomDouble #-} -- | Yields a set of random from the internal generator, -- using 'randomWord64' internally. -randomBits :: (Bits x, Member (State Random) r) +randomBits :: ( Member (State Random) r + , Bits x) => Eff r x randomBits = do let z = clearBit (bit 0) 0 -- zero, so we can get the number of bits @@ -172,14 +128,10 @@ return $ take k (concatMap breakBits word64s) -- | Returns the maximum set bit in an integer. -maxBit :: Integer -> Int -maxBit x - | x == 0 = 0 - | otherwise = - let loop y !k - | y == 1 = k - | otherwise = loop (y `unsafeShiftR` 1) (k+1) - in loop x 0 +maxBit :: Num a => Integer -> a +maxBit = let loop !k 0 = k + loop !k x = loop (k+1) (x `unsafeShiftR` 1) + in loop 0 -- | Repeat a computation until it succeeds a test. loopUntil :: Monad m => (a -> Bool) -> m a -> m a @@ -238,8 +190,8 @@ -- since it relaxes type constraints, but passing in -- constant bounds such as @uniformIntegralDist 0 10@ -- will warn with -Wall. -uniformIntegralDist :: (Member (State Random) r - , Integral a) +uniformIntegralDist :: ( Member (State Random) r + , Integral a ) => a -- ^ a -> a -- ^ b -> Eff r a @@ -260,7 +212,7 @@ return (d * range + a) -- | Generates a uniformly distributed random number in --- the inclusive range [a, b]. +-- the range [a, b). -- -- NOTE: This code might not be correct, in that the -- returned value may not be perfectly uniformly @@ -280,6 +232,36 @@ in uniformRealDist' a range {-# INLINE uniformRealDist #-} +-- | Generates a linearly-distributed random number in the range @[a, b)@; +-- @a@ with a probability of 0. +-- This code is not guaranteed to be correct. +linearRealDist :: Member (State Random) r + => Double + -> Double + -> Eff r Double +linearRealDist a' b' = do + let (a, b) = (min a' b', max a' b') + range = b - a + -- P(x) = + -- { 0 x < a + -- { k*(x-a) a <= x < b + -- { 0 b <= x + -- CDF(x) = + -- { 0 x < a + -- { k*(1/2)*(x-a)^2 a <= x < b + -- { 1 b <= x + -- Choose k such that CDF is continuous: + -- lim + -- x->b- CDF(x) = k*(1/2)*(b-a)^2 = 1 + -- k = 2/(b-a)^2 = 2/range^2 + -- Thus CDF(x, a<=x (x, a<=x. -sampleContDist :: (Member (State Random) r - , ContDistr d) +sampleContDist :: ( Member (State Random) r + , ContDistr d ) => d -- ^ The distribution to sample. -> Eff r Double sampleContDist d = @@ -335,7 +317,7 @@ -- | The value represents the number of failures in a series of -- independent yes/no trials (each succeeds with probability p), --- before exactly k successes occur. +-- before exactly k successes occur. -- -- p must be in the range (0, 1] -- k must be >= 0 @@ -406,15 +388,16 @@ gammaDist alpha beta = sampleContDist (DG.gammaDistr alpha beta) --- | ??? --- --- Warning: NOT IMPLEMENTED! +-- | Generates random numbers as sampled from a Weibull +-- distribution. It was originally identified to describe +-- particle size distribution. weibullDist :: Member (State Random) r => Double -- ^ α. The shape parameter. -> Double -- ^ β. The scale parameter. -> Eff r Double -weibullDist = - error "system-random-effect: weibullDist: TODO: Patches welcome!" +weibullDist alpha beta = do -- lambda = beta, k = alpha + r <- randomDouble + return $ beta * (-log (1-r)) ** (1/alpha) -- | ??? -- @@ -426,6 +409,7 @@ extremeValueDist = error "system-random-effect: extremeValueDist: TODO: Patches welcome!" + -- | Generates random numbers as sampled from the -- normal distribution. normalDist :: Member (State Random) r @@ -485,7 +469,7 @@ -- measurements, each with additive errors of unknown standard -- deviation, as in physical measurements. Or, alternatively, -- when estimating the unknown mean of a normal distribution --- with unknown standard deviation, given n+1 samples. +-- with unknown standard deviation, given n+1 samples. studentTDist :: Member (State Random) r => Double -- ^ The number of degrees of freedom -> Eff r Double @@ -526,30 +510,87 @@ mv <- V.unsafeThaw xs binarySearch mv y +-- Select an interval at random from a collection of intervals and +-- their weights, then run a provided function between its endpoints. +piecewiseDist :: Member (State Random) r + => (Double -> Double -> Eff r Double) + -> [Double] + -> DiscreteDistributionHelper + -> Eff r Double +piecewiseDist f intervals weights@(DDH rs) + | V.length rs + 1 /= length intervals = + error $ "system-random-effect: piecewiseDist:" + ++ " Incongruent parameter lengths." + ++ " intervals=" ++ show intervals + ++ " weights=" ++ show rs + | otherwise = do + idx <- discreteDist weights + + let vints = V.fromList intervals + (l, r) = (vints V.! idx, vints V.! (idx+1)) + f l r +{-# INLINE piecewiseDist #-} + -- | This function produces random floating-point numbers, which -- are uniformly distributed within each of the several subintervals -- [b_i, b_(i+1)), each with its own weight w_i. The set of interval -- boundaries and the set of weights are the parameters of this -- distribution. -- --- For example, `piecewiseConstantDist [ 0, 1, 10, 15 ] --- (buildDDH [ 1, 0, 1 ])` +-- For example, @piecewiseConstantDist [ 0, 1, 10, 15 ] +-- (buildDDH [ 1, 0, 1 ])@ -- will produce values between 0 and 1 half the time, and values -- between 10 and 15 the other half of the time. piecewiseConstantDist :: Member (State Random) r => [Double] -- ^ Intervals -> DiscreteDistributionHelper -- ^ Weights -> Eff r Double -piecewiseConstantDist intervals weights@(DDH rs) - | V.length rs + 1 /= length intervals = - error $ "system-random-effect: piecewiseConstantDist:" - ++ " Incongruent parameter lengths." - ++ " intervals=" ++ show intervals - ++ " weights=" ++ show rs - | otherwise = do - idx <- discreteDist weights +piecewiseConstantDist = piecewiseDist uniformRealDist - let vints = V.fromList intervals - (l, r) = (vints V.! idx, vints V.! (idx+1)) +-- | This function produces random floating-point numbers, which +-- are distributed with linearly-increasing probability within +-- each of the several subintervals [b_i, b_(i+1)), each with its own +-- weight w_i. The set of interval boundaries and the set of weights +-- are the parameters of this distribution. +-- +-- For example, `piecewiseLinearDist [ 0, 1, 10, 15 ] +-- (buildDDH [ 1, 0, 1 ])` +-- will produce values between 0 and 1 half the time, and values +-- between 10 and 15 the other half of the time. +piecewiseLinearDist :: Member (State Random) r + => [Double] -- ^ Intervals + -> DiscreteDistributionHelper -- ^ Weights + -> Eff r Double +piecewiseLinearDist = piecewiseDist linearRealDist - uniformRealDist l r +-- | Shuffle a mutable vector. +knuthShuffleM :: ( PrimMonad m + , Applicative m + , Typeable1 m + , Member (State Random) r + , SetMember Lift (Lift m) r + ) + => MVector (PrimState m) a + -> Eff r () +knuthShuffleM v = forM_ [0..iLast] + $ \i -> uniformIntegralDist 0 iLast + >>= lift . swap i + where + iLast = M.length v - 1 + + swap i j = do + [vi, vj] <- mapM (M.read v) [i, j] + mapM_ (uncurry $ M.write v) [(i, vj), (j, vi)] + +-- | Shuffle an immutable vector. +knuthShuffle :: Member (State Random) r + => Vector a + -> Eff r (Vector a) +knuthShuffle v0 = foldM swap v0 [0..iLast] + where + iLast = V.length v0 - 1 + + swap :: Member (State Random) r => Vector a -> Int -> Eff r (Vector a) + swap v i = do + j <- uniformIntegralDist 0 iLast + return $ v // zip [j, i] ((v !) <$> [i, j]) diff --git a/src/System/Random/Effect/Raw.hs b/src/System/Random/Effect/Raw.hs new file mode 100644 --- /dev/null +++ b/src/System/Random/Effect/Raw.hs @@ -0,0 +1,146 @@ +-- | Contains a random number generator for the deterministic, +-- fast, but NOT SECURE mt19937 random number generator. This +-- is a very versatile (and pure!) generator that works for +-- pretty much everything EXCEPT security. +-- +-- Seriously. Don't use this for crypto. +-- +-- This is the minimal definition for a random number generator. +-- new ones can be implemented by adding to the sum type 'Random'. +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE BangPatterns #-} +module System.Random.Effect.Raw ( Random + -- * Seeding + , mkRandom + , mkRandomIO + , mkSecureRandomIO + -- * Running + , forRandEff + , runRandomState + -- * Raw Generators + , randomInt + , randomWord + , randomInt64 + , randomWord64 + ) where + +import Control.Applicative +import Control.Arrow ( second ) +import Data.Bits +import qualified Data.ByteString as B +import Data.Typeable +import Data.Int +import Data.Word + +import qualified Crypto.Random as CR +import qualified System.Random.Mersenne.Pure64 as SR + +import Control.Eff +import Control.Eff.Lift +import Control.Eff.State.Strict + +-- | A random number generator. Either a fast, insecure mersenne +-- twister or a secure one, depending on which smart constructor +-- is used to construct this type. +data Random = FastRandom {-# UNPACK #-} !SR.PureMT + | SecureRandom {-# UNPACK #-} !CR.SystemRandom + deriving Typeable + +randomInt :: Member (State Random) r + => Eff r Int +randomInt = randomF SR.randomInt srandomBits +{-# INLINE randomInt #-} + +randomWord :: Member (State Random) r + => Eff r Word +randomWord = randomF SR.randomWord srandomBits +{-# INLINE randomWord #-} + +randomInt64 :: Member (State Random) r + => Eff r Int64 +randomInt64 = randomF SR.randomInt64 srandomBits +{-# INLINE randomInt64 #-} + +randomWord64 :: Member (State Random) r + => Eff r Word64 +randomWord64 = randomF SR.randomWord64 srandomBits +{-# INLINE randomWord64 #-} + +-- | Create a random number generator from a 'Word64' seed. +-- This uses the insecure (but fast) mersenne twister. +mkRandom :: Word64 -> Random +mkRandom = FastRandom . SR.pureMT +{-# INLINE mkRandom #-} + +-- | Create a new random number generator, using the clocktime as the base for +-- the seed. This must be called from a computation with a lifted base effect +-- of 'IO'. +-- +-- This is just a conveniently seeded mersenne twister. +mkRandomIO :: SetMember Lift (Lift IO) r + => Eff r Random +mkRandomIO = + FastRandom <$> lift SR.newPureMT +{-# INLINE mkRandomIO #-} + +-- | Creates a new random number generator, using the system entropy source +-- as a seed. The random number generator returned from this function is +-- cryptographically secure, but not nearly as fast as the one returned +-- by 'mkRandom' and 'mkRandomIO'. +mkSecureRandomIO :: SetMember Lift (Lift IO) r + => Eff r Random +mkSecureRandomIO = do + SecureRandom <$> lift CR.newGenIO +{-# INLINE mkSecureRandomIO #-} + +-- | Use a non-random effect as the Random source for running a random effect. +forRandEff :: Eff r Random -> Eff (State Random :> r) w -> Eff r w +forRandEff rndgen e = rndgen >>= (`runRandomState` e) +{-# INLINE forRandEff #-} + +-- | Runs an effectful random computation, returning the computation's result. +runRandomState :: Random + -> Eff (State Random :> r) w + -> Eff r w +runRandomState seed computation = + snd <$> runState seed computation +{-# INLINE runRandomState #-} + +foldBits :: (Bits a, Num a) + => B.ByteString + -> a +foldBits bs = + let addByte byte bits = + (bits `unsafeShiftL` 8) .|. fromIntegral byte + in B.foldr' addByte 0 bs +{-# INLINE foldBits #-} + +-- | Securely generate some random bits. +srandomBits :: ( Bits a + , Num a ) + => CR.SystemRandom + -> (a, CR.SystemRandom) +srandomBits sr = + let z = clearBit (bit 0) 0 + nBytes = bitSize z `div` 8 + in case CR.genBytes nBytes sr of + Left err -> error $ "system-random-effect: System.Random.Effect.Secure: genBytes: " ++ show err + Right (bs, sr') -> (z .|. foldBits bs, sr') +{-# INLINE srandomBits #-} + +-- | A generalized form of generating a random number of the correct type +-- from System.Random.Mersenne.Pure64. +randomF :: Member (State Random) r + => (SR.PureMT -> (a, SR.PureMT)) + -> (CR.SystemRandom -> (a, CR.SystemRandom)) + -> Eff r a +randomF f s = do + old <- get + let (val, new) = case old of + (FastRandom r) -> second FastRandom (f r) + (SecureRandom r) -> second SecureRandom (s r) + put new + return val +{-# INLINE randomF #-} diff --git a/system-random-effect.cabal b/system-random-effect.cabal --- a/system-random-effect.cabal +++ b/system-random-effect.cabal @@ -1,5 +1,5 @@ name: system-random-effect -version: 0.3.1 +version: 0.4.0 synopsis: Random number generation for extensible effects. homepage: https://github.com/wowus/system-random-effect license: BSD3 @@ -14,6 +14,7 @@ hs-source-dirs: src/ ghc-options: -Wall exposed-modules: System.Random.Effect + other-modules: System.Random.Effect.Raw build-depends: base == 4.6.* , extensible-effects == 1.2.* @@ -21,6 +22,9 @@ , statistics == 0.10.* , vector == 0.10.* , vector-algorithms == 0.5.* + , primitive == 0.5.* + , crypto-api == 0.12.* + , bytestring == 0.10.* default-language: Haskell2010 @@ -29,9 +33,10 @@ main-is: Test.hs hs-source-dirs: test/ - ghc-options: -rtsopts=all -threaded + ghc-options: -rtsopts=all -threaded -O2 build-depends: base == 4.6.* + , vector == 0.10.* , QuickCheck == 2.* , HUnit == 1.2.* , test-framework == 0.8.* @@ -47,11 +52,13 @@ hs-source-dirs: bench/ ghc-options: -Wall -rtsopts=all -threaded -fno-warn-orphans -O2 main-is: Bench.hs + x-uses-tf: true build-depends: base == 4.6.* , deepseq == 1.3.* , criterion == 0.8.* , extensible-effects == 1.2.* + , vector == 0.10.* , system-random-effect default-language: Haskell2010 diff --git a/test/Test.hs b/test/Test.hs --- a/test/Test.hs +++ b/test/Test.hs @@ -1,22 +1,28 @@ {-# LANGUAGE TypeOperators #-} {-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE BangPatterns #-} module Main ( main ) where +import Prelude hiding (all) + import Control.Eff +import Control.Eff.Lift import Control.Eff.State.Strict import System.Random.Effect -import Control.Monad (void) +import Control.Applicative +import Control.Monad.ST +import Data.Vector ( Vector ) +import qualified Data.Vector as V +import qualified Data.Vector.Mutable as MV import Data.Word -import Data.Typeable -import Test.Framework (defaultMain, testGroup) -import Test.Framework.Providers.HUnit +import Test.Framework ( defaultMain, Test ) import Test.Framework.Providers.QuickCheck2 -import Test.HUnit hiding (State) import Test.QuickCheck +import Test.QuickCheck.Property ( morallyDubiousIOProperty ) main :: IO () main = defaultMain tests @@ -24,6 +30,9 @@ runWithSeed :: Word64 -> Eff (State Random :> ()) a -> a runWithSeed seed = run . runRandomState (mkRandom seed) +runIOWithSeed :: Word64 -> Eff (State Random :> Lift IO :> ()) a -> IO a +runIOWithSeed seed = runLift . runRandomState (mkRandom seed) + checkRange :: (Integer, Integer) -> Integer -> Bool checkRange (low, high) x = x >= low && x <= high @@ -36,33 +45,37 @@ in checkRange (low, high) . runWithSeed seed $ do uniformIntDist a b -testDiscreteDistributionInRange :: [Word64] -> Word64 -> Bool -testDiscreteDistributionInRange xs seed = +newtype DiscreteWeights = DW [Word64] + deriving Show + +instance Arbitrary DiscreteWeights where + arbitrary = DW <$> suchThat (listOf arbitrary) ((> 0) . sum) + shrink (DW xs) = map DW (shrink xs) + +testDiscreteDistributionInRange :: DiscreteWeights -> Word64 -> Bool +testDiscreteDistributionInRange (DW xs) seed = let ddh = buildDDH xs minVal = 0 maxVal = length xs - 1 - in sum xs == 0 || - ((\x -> x >= minVal && x <= maxVal) . runWithSeed seed $ - discreteDist ddh) + in (\x -> x >= minVal && x <= maxVal) . runWithSeed seed $ + discreteDist ddh -testNoZeroDiscreteDistributionPick :: [Word64] -> Word64 -> Bool -testNoZeroDiscreteDistributionPick xs seed = +testNoZeroDiscreteDistributionPick :: DiscreteWeights -> Word64 -> Bool +testNoZeroDiscreteDistributionPick (DW xs) seed = let ddh = buildDDH xs - in sum xs == 0 || - ((\x -> (xs !! x) /= 0) . runWithSeed seed $ - discreteDist ddh) + in (\x -> (xs !! x) /= 0) . runWithSeed seed $ + discreteDist ddh -testUnsafeThaw :: [Word64] -> Word64 -> Bool -testUnsafeThaw xs seed = +testUnsafeThaw :: DiscreteWeights -> Word64 -> Bool +testUnsafeThaw (DW xs) seed = let ddh = buildDDH xs - in sum xs == 0 || - (runWithSeed seed $ do + in runWithSeed seed $ do _ <- discreteDist ddh _ <- discreteDist ddh _ <- discreteDist ddh _ <- discreteDist ddh _ <- discreteDist ddh - return True) + return True testUniformIntegralDist :: Integer -> Integer -> Word64 -> Bool testUniformIntegralDist a b seed = @@ -70,10 +83,88 @@ r2 = runWithSeed seed $ uniformIntegralDist a b in r1 == r2 +testKnuthShuffle :: [Int] -> Word64 -> Bool +testKnuthShuffle xs' seed = + let xs = V.fromList xs' + countIf f = V.length . V.filter f + shuffled = runWithSeed seed (knuthShuffle xs) + sameCount v1 v2 = V.all id + $ V.map (\x -> countIf (== x) v1 + == countIf (== x) v2) v1 + in sameCount xs shuffled + +testKnuthShuffleM :: [Int] -> Word64 -> IO Bool +testKnuthShuffleM xs' seed = do + let xs = V.fromList xs' + countIf f = V.length . V.filter f + shuffled = do + vs <- V.thaw xs + runIOWithSeed seed (knuthShuffleM vs) + V.freeze vs + sameCount v1 v2 = V.all id + $ V.map (\x -> countIf (== x) v1 + == countIf (== x) v2) v1 + shuf <- shuffled + return (sameCount xs shuf) + +testKnuthShuffleEquivalence :: [Int] -> Word64 -> IO Bool +testKnuthShuffleEquivalence xs seed = do + let vs = V.fromList xs + ks1 = runWithSeed seed (knuthShuffle vs) + + xs' <- V.thaw vs + runIOWithSeed seed (knuthShuffleM xs') + ks2 <- V.freeze xs' + + return (ks1 == ks2) + +testSecureRandom :: Integer -> Integer -> IO Bool +testSecureRandom a b = do + let low = min a b + high = max a b + + runLift $ do + rng <- mkSecureRandomIO + return $ run $ runRandomState rng $ + checkRange (low, high) <$> uniformIntDist a b + +(|>) :: b -> (b -> c) -> c +(|>) = flip ($) + +histogram :: Vector Integer -> Vector Int +histogram v = runST $ do + mv <- MV.replicate (fromIntegral (V.maximum v + 1)) 0 + V.forM_ v $ \i' -> do + let i = (fromIntegral i') :: Int + !x <- MV.read mv i + MV.write mv i (x+1) + V.unsafeFreeze mv + +-- check if all uniformly distributed numbers are within 10% of the mean. +-- 10% was a number chosen arbitrarily. +simpleUniformIntDistTest :: Word64 -> Bool +simpleUniformIntDistTest seed = + let nBuckets = 5 :: Int + samplesPerBucket = 4000 :: Int + nSamples = nBuckets * samplesPerBucket + maxDelta = (fromIntegral samplesPerBucket) `div` 10 + nums = uniformIntDist 0 (fromIntegral (nBuckets - 1)) + |> V.replicateM nSamples + |> runWithSeed seed + hist = histogram nums + in V.all (\x -> x >= samplesPerBucket - maxDelta + && x <= samplesPerBucket + maxDelta) hist + +tests :: [Test] tests = [ testProperty "random range" testUniformRandom , testProperty "discrete dist range" testDiscreteDistributionInRange , testProperty "no non-zero discrete dist pick" testNoZeroDiscreteDistributionPick , testProperty "unsafeThaw is okay to use" testUnsafeThaw , testProperty "testUniformIntegralDist == testUniformIntDist" testUniformIntegralDist + , testProperty "knuth shuffle" testKnuthShuffle + , testProperty "knuth shuffle (monadic)" (\xs seed -> morallyDubiousIOProperty $ testKnuthShuffleM xs seed) + , testProperty "knuth shuffle equivalence" (\xs seed -> morallyDubiousIOProperty $ testKnuthShuffleEquivalence xs seed) + , testProperty "secure random" (\a b -> morallyDubiousIOProperty $ testSecureRandom a b) + , testProperty "uniformIntDist is uniform-ish" simpleUniformIntDistTest ]