random-variates 0.1.1.0 → 0.1.3.0

raw patch · 7 files changed

+224/−74 lines, 7 filesdep +binarydep +bytestringdep ~basenew-component:exe:GenPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: binary, bytestring

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• random-variates.cabal +19/−9
• src/Stochastic/Distributions.hs +30/−6
• src/Stochastic/Distributions/Continuous.hs +0/−11
• src/Stochastic/Distributions/Discrete.hs +10/−12
• src/Stochastic/Generator.hs +52/−36
• src/Stochastic/Uniform.hs +69/−0
• tests/gen.hs +44/−0

@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                random-variates-version:             0.1.1.0+version:             0.1.3.0 synopsis:            "Uniform RNG => Non-Uniform RNGs" description:         "Collection of transforms uniform random number generators (RNGs) into any of a dozen common RNGs. Each presenting several common interfaces. Additionally Empirical distributions can be sampled from and tested (chi-squared) against theoretical distributions."    license:             MIT@@ -23,6 +23,7 @@ library   exposed-modules:                   Stochastic.Generator+                  Stochastic.Uniform                   Stochastic.Distribution.Continuous                   Stochastic.Distribution.Discrete                   Stochastic.Distributions@@ -33,16 +34,25 @@    other-modules:    build-depends:       -                         base >=4.6 && <5.0-                       , containers >= 0.5.7.0-                       , lens >=4.13-                       , random >=1.1-                       , reinterpret-cast >= 0.1.0-                       , mtl >= 2.2-                       , erf >= 2.0+                       base >=4.6 && <5.0+                     , containers >= 0.5.7.0+                     , lens >=4.13+                     , random >=1.1+                     , reinterpret-cast >= 0.1.0+                     , mtl >= 2.2+                     , erf >= 2.0+                     , bytestring >= 0.10+                     ,  binary    hs-source-dirs:      src   default-language:    Haskell2010+  default-extensions: DeriveGeneric,  DeriveDataTypeable +Executable Gen+  main-is: tests/gen.hs+  build-depends:       base+                     , random-variates >=0.1+  default-language:    Haskell2010+ Test-Suite vis   type:                exitcode-stdio-1.0   main-is:             tests/vis.hs@@ -60,4 +70,4 @@                      , random-variates >=0.1   default-language:    Haskell2010 -  +   

@@ -2,28 +2,52 @@ module Stochastic.Distributions(   UniformBase(rDouble)   ,stdBase+  ,seededBase   ,Empirical(..)   ,mkEmpirical   ) where +import Stochastic.Uniform+import qualified Data.ByteString.Char8 as B+import qualified Data.ByteString.Lazy as LBS+ import System.Random import Control.Monad.State.Lazy import Stochastic.Tools+import System.IO+import Data.Word+import Data.Binary.Get  data UniformBase = UniformBase {   rDouble :: (Double, UniformBase) }  -stdGen2Uni gen = UniformBase {-  rDouble = mapTuple-            (id)-            (stdGen2Uni)-            (randomR (0,1) gen)+readWord64 :: Handle -> IO Word64+readWord64 h = do+  w1 <- hGetChar h+  w2 <- hGetChar h+  w3 <- hGetChar h+  w4 <- hGetChar h+  w5 <- hGetChar h+  w6 <- hGetChar h+  w7 <- hGetChar h+  w8 <- hGetChar h+  let words = [w1,w2,w3,w4,w5,w6,w7,w8] :: String+  return $ runGet getWord64host $ LBS.fromStrict $ B.pack words++mk g = UniformBase {+  rDouble = mapTuple (id) (mk) (random g)   } -stdBase s = stdGen2Uni (mkStdGen s)+stdBase :: Integer -> UniformBase+stdBase s = mk $ xorshift128plus s +seededBase :: IO UniformBase+seededBase = do+  word <- withBinaryFile "/dev/random/" ReadMode (readWord64)+  let seed = toInteger word+  return $ stdBase seed  data Empirical = Empirical {   degreesOfFreedom :: Int,

@@ -1,5 +1,3 @@-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE TypeFamilies #-} module Stochastic.Distributions.Continuous(   mkUniform   ,mkExp@@ -25,15 +23,6 @@   cdf  _ x = x   cdf' _ p = p   degreesOfFreedom _ = 0---instance Generator Dist where-  type (From Dist) = Double-  nextG = state $ \ g0 -> rand g0--instance Generator UniformBase where-  type (From UniformBase) = Double-  nextG = state $ \ g0 -> rDouble g0   data Dist =

@@ -1,4 +1,3 @@-{-# LANGUAGE TypeFamilies #-} module Stochastic.Distributions.Discrete(   mkBinomial   ,mkBernoulli@@ -12,15 +11,11 @@ import Data.Maybe import Control.Monad.State.Lazy import Stochastic.Tools-import Stochastic.Generator(Generator(..), foldWhile)+import Stochastic.Generator(foldGenWhile) import Stochastic.Distributions import Stochastic.Distribution.Discrete import qualified Stochastic.Distributions.Continuous as C -instance Generator Dist where-  type (From Dist) = Int-  nextG = state $ \ g0 -> rand g0- data Dist =   Uniform Int Int UniformBase   | Poisson C.Dist@@ -96,7 +91,7 @@     mapTuple     (\x -> length x)     (Poisson)-    (runState (foldWhile (+) 0.0 (<1.0)) g0)+    ((foldGenWhile (C.rand) (+) 0.0 (<1.0)) g0)   rand (Bernoulli p g0) =     mapTuple     (\x -> if (x >= p) then 1 else 0)@@ -132,12 +127,12 @@   cdf (Uniform a b _) x = toDbl (x-a) / toDbl (b-a)    cdf' g@(Poisson (C.Exponential y _)) x =-    sum . fst $ runState (fold) [1..]+    (sum . fst) (fold [1..])     where-      reduce :: Int -> Double -> Double-      reduce = (\y p -> (pmf g y) + p)-      fold :: State [Int] [From [Int]]-      fold = foldWhile (reduce) 0 (<x)+      reduce :: Double -> Int -> Double+      reduce = (\p y -> (pmf g y) + p)+      fold :: [Int] -> ([Int], [Int])+      fold = foldGenWhile (myUncons) (reduce) 0 (<x)   cdf' (Geometric p _) x =     ceiling $ (log (1-x)) / (log (1-p))   cdf' (Bernoulli p _) x@@ -166,6 +161,9 @@       (_, r, _) = head $ filter (\(w,_,_) -> (w==k)) cache   pmf (Uniform a b _) x = toDbl x/toDbl (b-a) ++myUncons :: [a] -> (a, [a])+myUncons (x:xs) = (x, xs)  toDbl = fromInteger . toInteger 

@@ -1,45 +1,61 @@-{-# LANGUAGE TypeFamilies #-} module Stochastic.Generator where-+import Control.Concurrent.MVar import Control.Monad.State.Lazy -class Generator g where-  type From g-  nextG   :: State g (From g)-  nextN  :: Int -> State g [(From g)]-  nextN 0 = state $ \g0 -> ([], g0)-  nextN n =-    do-      x  <- nextG-      xs <- nextN (n-1)-      return (x:xs)+type Gen g a = (g -> (a,g))+data IOGen g a = IOGen (g -> (a, g)) (MVar g) -instance Generator [a] where-  type From [a] = a-  nextG = state $ \ g0 -> (head g0, tail g0)+liftGen :: (g -> (a, g)) -> g -> IO (IOGen g a)+liftGen f g = do+  var <- newMVar g+  return $ IOGen f var -foldWhile :: Generator g-         => (From g -> a -> a)-         -> a-         -> (a -> Bool)-         -> State g [From g]-foldWhile f z p =-  do-   x  <- nextG-   let y = f x z-   xs <- if (p y)-         then foldWhile f y p -         else return []-   return (x:xs)+nextIO :: IOGen g a -> IO a+nextIO (IOGen f var) = do+  val <- takeMVar var+  let (x, g') = f val+  putMVar var g'+  return x -while :: Generator g => ((From g) -> Bool) -> State g [From g]-while p =-  do-   x  <- nextG-   xs <- if (p x)-         then while p-         else return []-   return (x:xs)+foldGenWhile :: (g -> (a,g))+          -> (b -> a -> b)+          -> b+          -> (b -> Bool)+          -> (g -> ([a], g))+foldGenWhile nxt f zz p = h zz+  where+    h z g0 +      | not (p z) = ([], g0)+      | otherwise = (x:xs, g2)+      where+        (xs, g2) = h (f z x) g1+        (x, g1)  = nxt g0 +genWhile :: (g -> (a, g)) -> (a -> Bool) -> (g -> ([a], g))+genWhile nxt p = h+  where+    h g0 = let (x, g1) = nxt g0 in+            let (xs, g2) = h g1 in+            if (p x) then (x:xs, g2) else ([], g1) +genTake :: (g -> (a,g)) -> Integer -> (g -> ([a], g))+genTake f 0 g0 = ([], g0)+genTake f n g0 = ((x:xs), g2)+  where+    (x, g1) = f g0+    (xs, g2) = genTake f (n-1) g1++dropGen :: (g -> (a,g)) -> Integer -> g -> g+dropGen f = d+  where+    d 0 g0 = g0+    d n g0 = d (n-1) $! (snd $ f g0)++++dropIO :: IOGen g a -> Integer -> IO ()+dropIO _ 0 = return ()+dropIO ioG n = do+  nextIO ioG+  dropIO ioG (n-1) 

@@ -0,0 +1,69 @@+{-# LANGUAGE BangPatterns #-}+module Stochastic.Uniform(xorshift128plus,+                          UniformRandom,+                          nWayAllocate,+                          splitAllocate,+                          RandomGen(..)) where++import Stochastic.Generator+import Data.Word+import Data.Bits+import Data.Typeable+import Control.Exception(throw, Exception)+import System.Random(RandomGen(..))++data UniformRandom = XorShift128Plus Word64 Word64 Integer++data EntropyExhausted = EntropyExhausted+                      deriving(Eq, Typeable)++instance Exception EntropyExhausted where+instance Show EntropyExhausted where+  show e = "EntropyExhausted"++{- |+For information on the performance of the xorshift-128-plus PRNG, please see: <http://vigna.di.unimi.it/ftp/papers/xorshiftplus.pdf Vigna et al.>+-}+xorshift128plus :: Integer -> UniformRandom+xorshift128plus seed = XorShift128Plus high low entropy+  where+    high = fromInteger seed+    low = fromInteger seed+    entropy = (2^127)++nWayAllocate :: Integer -> Integer -> UniformRandom -> ([UniformRandom], UniformRandom)+nWayAllocate _ 0 g0 = ([], g0)+nWayAllocate size n g0 = ((g1:gs), g3)+  where+    !(gs,g3) = nWayAllocate size (n-1) g2+    !(g1,g2) = splitAllocate size g0++splitAllocate :: Integer -> UniformRandom -> (UniformRandom, UniformRandom)+splitAllocate count g@(XorShift128Plus high low entropy) =+    ((XorShift128Plus high low count), g')+    where+      !g' = step (next) (count) g+++instance RandomGen UniformRandom where+  next (XorShift128Plus high low entropy) +    | entropy == 0 = throw EntropyExhausted+    | otherwise    = final +    where+      -- eagerly evaluate this function, retain no intermediaries or we might blow the stack+      final    = (ret, XorShift128Plus high' low' entropy')+      !ret     = fromInteger $ toInteger $ high' + high+      x        = low `xor` (low `shift` 23)+      high'    = x `xor` high `xor` (x `shift` (-17)) `xor` (high `shift` (-26))+      low'     = high+      entropy' = entropy - 1+  split g@(XorShift128Plus high low entropy) =+    ((XorShift128Plus high low entropy'), g')+    where+      !entropy' = (2^32)+      !g' = step (next) (entropy') g++step :: (g -> (a,g)) -> Integer -> g -> g+step f 0 g'' = g''+step f n g'' = step (f) (n-1) $! (snd $ f g'')+

@@ -0,0 +1,44 @@+{-# LANGUAGE BangPatterns #-}+module Main where++import Stochastic.Generator+import Stochastic.Uniform+import Stochastic.Distributions+import System.Environment+import System.IO+import Data.IORef++main :: IO ()+main = do+  [count] <- getArgs+  let n  = read count :: Integer+  let g  = xorshift128plus 42+  let g' = dropGen (next) n g+  iog <- liftGen (next) g'+  x <- nextIO iog+  putStrLn (show x)+--  ref <- newIORef g+--  invert $ step n (f ref)+  return ()+++myUncons (x:xs) = (x, xs)++step :: Integer -> IO Int -> [IO Int]+step 0 _  = []+step n io = (io:(step (n-1) io))++invert :: [IO a] -> IO [a]+invert [] = return []+invert (io:ios) = do+  x <- io+  xs <- invert ios+  return (x:xs)++f :: IORef UniformRandom -> IO Int+f ref = do+  g <- readIORef ref+  let (x, g') = next g+  writeIORef ref g'+  return x+