random-fu 0.0.2.1 → 0.0.3
raw patch · 7 files changed
+175/−88 lines, 7 filesdep ~stateref
Dependency ranges changed: stateref
Files
- random-fu.cabal +7/−6
- src/Data/Random/Distribution/Discrete.hs +72/−41
- src/Data/Random/Distribution/Normal.hs +5/−7
- src/Data/Random/Distribution/Uniform.hs +43/−17
- src/Data/Random/Internal/Fixed.hs +22/−0
- src/Data/Random/Source/PureMT.hs +13/−9
- src/Data/Random/Source/StdGen.hs +13/−8
random-fu.cabal view
@@ -1,5 +1,5 @@ name: random-fu-version: 0.0.2.1+version: 0.0.3 stability: experimental cabal-version: >= 1.2@@ -27,27 +27,28 @@ Data.Random.Distribution.Bernoulli Data.Random.Distribution.Beta Data.Random.Distribution.Binomial- Data.Random.Distribution.Exponential Data.Random.Distribution.Discrete+ Data.Random.Distribution.Exponential Data.Random.Distribution.Gamma Data.Random.Distribution.Normal- Data.Random.Distribution.Ziggurat Data.Random.Distribution.Poisson Data.Random.Distribution.Rayleigh Data.Random.Distribution.Triangular Data.Random.Distribution.Uniform+ Data.Random.Distribution.Ziggurat Data.Random.Internal.Find+ Data.Random.Internal.Fixed Data.Random.Internal.TH Data.Random.Internal.Words- Data.Random.List Data.Random.Lift+ Data.Random.List Data.Random.RVar Data.Random.Sample Data.Random.Source Data.Random.Source.DevRandom- Data.Random.Source.StdGen Data.Random.Source.PureMT Data.Random.Source.Std+ Data.Random.Source.StdGen if flag(base4) build-depends: base >= 4 && <5, syb@@ -62,6 +63,6 @@ mtl, random, random-shuffle,- stateref,+ stateref >= 0.3 && < 0.4, storablevector, template-haskell
src/Data/Random/Distribution/Discrete.hs view
@@ -13,8 +13,11 @@ import Data.Random.Distribution.Uniform import Data.Random.List (randomElement) +import Control.Arrow import Control.Monad import Control.Applicative+import Data.Foldable (Foldable(foldMap))+import Data.Traversable (Traversable(traverse, sequenceA)) import Data.List import Data.Function@@ -22,6 +25,12 @@ discrete :: Distribution (Discrete p) a => [(p,a)] -> RVar a discrete ps = rvar (Discrete ps) +empirical :: (Num p, Ord a) => [a] -> Discrete p a+empirical xs = Discrete bins+ where bins = [ (genericLength bin, x)+ | bin@(x:_) <- group (sort xs)+ ]+ newtype Discrete p a = Discrete [(p, a)] deriving (Eq, Show) @@ -33,11 +42,21 @@ when (any (<0) ps) $ fail "negative probability in discrete distribution" - let totalProb = last cs- if totalProb <= 0- then randomElement xs -- this probably makes the monad instance incorrect for discarding zero-probability events...+ let totalWeight = last cs+ if totalWeight <= 0+ -- if all events are "equally impossible", just pick an arbitrary one.+ then randomElement xs else do- u <- uniform 0 totalProb+ let getU = do+ u <- uniform 0 totalWeight+ -- reject 0; this causes integral weights to behave sensibly (although it + -- potentially wastes up to 50% of all sampled integral values in the+ -- case where totalWeight = 1) and is still valid for fractional weights+ -- (it only prevents zero-probability events from ever occurring, which + -- is reasonable).+ if u == 0 then getU+ else return u+ u <- getU return $ head [ x | (c,x) <- zip cs xs@@ -47,59 +66,71 @@ instance Functor (Discrete p) where fmap f (Discrete ds) = Discrete [(p, f x) | (p, x) <- ds] +instance Foldable (Discrete p) where+ foldMap f (Discrete ds) = foldMap (f . snd) ds++instance Traversable (Discrete p) where+ traverse f (Discrete ds) = Discrete <$> traverse (\(p,e) -> (\e -> (p,e)) <$> f e) ds+ sequenceA (Discrete ds) = Discrete <$> traverse (\(p,e) -> (\e -> (p,e)) <$> e) ds+ -- We want each subset of cases in fx derived from a given case --- in x to have the same total weight as the set in x from whence they came.------ thus, w(f x) == w (x) is sufficient (although not necessary), where w() is--- the weight.------ TODO: consider establishing normalization invariant?--- TODO: Consider what should happen when f x returns multiple events with zero weight-instance (Fractional p, Ord p) => Monad (Discrete p) where+-- in x to have the same relative weight as the set in x from whence they came.+instance Num p => Monad (Discrete p) where return x = Discrete [(1, x)] (Discrete x) >>= f = Discrete $ do (p, x) <- x let Discrete fx = f x- let qx = fx- -- TODO - check out whether this is valid when not requiring normalization...- -- qx = [ (q, x)- -- | (q, x) <- fx- -- , q > 0 -- should this be done? Consider case where all results have 0 weight...- -- ]- qs = sum (map fst qx) -- either (qx == []) or (sum qs > 0)- scale - | qs > 0- = recip qs- | otherwise- = recip (fromIntegral (length qx))- - (q, x) <- qx+ (q, x) <- fx - return (p * q * scale, x)+ return (p * q, x) -instance (Fractional p, Ord p) => Applicative (Discrete p) where+instance Num p => Applicative (Discrete p) where pure = return (<*>) = ap +-- |Like 'fmap', but for the weights of a discrete distribution.+mapDiscreteWeights :: (p -> q) -> Discrete p e -> Discrete q e+mapDiscreteWeights f (Discrete ds) = Discrete [(f p, x) | (p, x) <- ds]++-- |Adjust all the weights of a discrete distribution so that they +-- sum to unity. If not possible, returns the original distribution +-- unchanged.+normalizeDiscreteWeights :: (Fractional p) => Discrete p e -> Discrete p e+normalizeDiscreteWeights orig@(Discrete ds) = + -- For practical purposes the scale factor is strict anyway,+ -- so check if it's 0 or 1 and, if so, skip the actual scaling part.+ if ws `elem` [0,1]+ then orig+ else Discrete+ [ (w * scale, e)+ | (w, e) <- ds+ ] + where+ ws = sum (map fst ds)+ scale = recip ws++-- |Simplify a discrete distribution by combining equivalent events (the new+-- event will have a weight equal to the sum of all the originals). collectDiscreteEvents :: (Ord e, Num p, Ord p) => Discrete p e -> Discrete p e-collectDiscreteEvents (Discrete ds) = - Discrete . concatMap (uncurry combine . unzip) . groupEvents . sortEvents $ ds+collectDiscreteEvents = collectDiscreteEventsBy compare sum head+ +-- |Simplify a discrete distribution by combining equivalent events (the new+-- event will have a weight equal to the sum of all the originals).+-- The comparator function is used to identify events to combine. Once chosen,+-- the events and their weights are combined (independently) by the provided+-- weight and event aggregation functions.+collectDiscreteEventsBy :: (e -> e -> Ordering) -> ([p] -> p) -> ([e] -> e)-> Discrete p e -> Discrete p e+collectDiscreteEventsBy compareE sumWeights mergeEvents (Discrete ds) = + Discrete . map ((sumWeights *** mergeEvents) . unzip) . groupEvents . sortEvents $ ds where- groupEvents = groupBy ((==) `on` snd)- sortEvents = sortBy (compare `on` snd)- - -- don't combine negative weights with positive ones.- -- don't error out ether - just leave them alone, it'll- -- all barf when the distribution is sampled.- combine ps (x:_) = case partition (> 0) (filter (/= 0) ps) of- ([], []) -> []- ([], ns) -> (sum ns, x) : []- (ps, []) -> (sum ps, x) : []- (ps, ns) -> (sum ps, x) : (sum ns, x) : []+ groupEvents = groupBy (\x y -> snd x `compareE` snd y == EQ)+ sortEvents = sortBy (compareE `on` snd) weight (p,x) | p < 0 = error "negative probability in discrete distribution" | otherwise = p event ((p,x):_) = x+ + combine (ps, xs) = (sumWeights ps, mergeEvents xs)
src/Data/Random/Distribution/Normal.hs view
@@ -118,7 +118,7 @@ doubleStdNormal :: RVar Double doubleStdNormal = runZiggurat doubleStdNormalZ --- doubleStdNormalC must not be over 12 if using wordToDoubleWithExcess+-- doubleStdNormalC must not be over 2^12 if using wordToDoubleWithExcess doubleStdNormalC :: Int doubleStdNormalC = 512 doubleStdNormalR, doubleStdNormalV :: Double@@ -140,6 +140,7 @@ floatStdNormal :: RVar Float floatStdNormal = runZiggurat floatStdNormalZ +-- floatStdNormalC must not be over 2^41 if using wordToFloatWithExcess floatStdNormalC :: Int floatStdNormalC = 512 floatStdNormalR, floatStdNormalV :: Float@@ -153,9 +154,6 @@ getIU (normalTail floatStdNormalR) where- -- p must not be over 41 if using wordToFloatWithExcess- p = 6- getIU = do w <- getRandomWord let (u,i) = wordToFloatWithExcess w@@ -164,8 +162,8 @@ normalPdf :: Real a => a -> a -> a -> Double normalPdf m s x = recip (realToFrac s * sqrt (2*pi)) * exp (-0.5 * (realToFrac x - realToFrac m)^2 / (realToFrac s)^2) -normalCdf :: (Real a, Erf a) => a -> a -> a -> Double-normalCdf m s x = realToFrac (normcdf ((x-m) / s))+normalCdf :: (Real a) => a -> a -> a -> Double+normalCdf m s x = normcdf ((realToFrac x - realToFrac m) / realToFrac s) data Normal a = StdNormal@@ -185,7 +183,7 @@ x <- floatStdNormal return (x * s + m) -instance (Real a, Erf a, Distribution Normal a) => CDF Normal a where+instance (Real a, Distribution Normal a) => CDF Normal a where cdf StdNormal = normalCdf 0 1 cdf (Normal m s) = normalCdf m s
src/Data/Random/Distribution/Uniform.hs view
@@ -19,10 +19,12 @@ , realFloatUniform , floatUniform , doubleUniform+ , fixedUniform , boundedStdUniform , boundedEnumStdUniform , realFloatStdUniform+ , fixedStdUniform , floatStdUniform , doubleStdUniform @@ -32,11 +34,13 @@ import Data.Random.Internal.TH import Data.Random.Internal.Words+import Data.Random.Internal.Fixed import Data.Random.Source import Data.Random.Distribution import Data.Random.RVar +import Data.Fixed import Data.Word import Data.Int import Data.List@@ -98,6 +102,14 @@ where one = 1 +fixedStdUniform :: HasResolution r => RVar (Fixed r)+fixedStdUniform = x+ where+ res = resolutionOf2 x+ x = do+ u <- uniform 0 (res)+ return (mkFixed u)+ realStdUniformCDF :: Real a => a -> Double realStdUniformCDF x | x <= 0 = 0@@ -122,6 +134,11 @@ x <- realFloatStdUniform return (a + x * (b - a)) +fixedUniform :: HasResolution r => Fixed r -> Fixed r -> RVar (Fixed r)+fixedUniform a b = do+ u <- integralUniform (unMkFixed a) (unMkFixed b)+ return (mkFixed u)+ realUniformCDF :: Real a => a -> a -> a -> Double realUniformCDF a b x | b < a = realUniformCDF b a x@@ -169,16 +186,16 @@ |]) -- Some integral types have specialized StdUniform rvars:-instance Distribution StdUniform Int8 where rvar StdUniform = fmap fromIntegral getRandomByte-instance Distribution StdUniform Word8 where rvar StdUniform = getRandomByte-instance Distribution StdUniform Word64 where rvar StdUniform = getRandomWord+instance Distribution StdUniform Int8 where rvar ~StdUniform = fmap fromIntegral getRandomByte+instance Distribution StdUniform Word8 where rvar ~StdUniform = getRandomByte+instance Distribution StdUniform Word64 where rvar ~StdUniform = getRandomWord -- and Integer has none... $( replicateInstances ''Int (integralTypes \\ [''Int8, ''Word8, ''Word64, ''Integer]) [d|- instance Distribution StdUniform Int where rvar StdUniform = fmap fromIntegral getRandomWord+ instance Distribution StdUniform Int where rvar ~StdUniform = fmap fromIntegral getRandomWord |]) $( replicateInstances ''Int (integralTypes \\ [''Integer]) [d|- instance CDF StdUniform Int where cdf StdUniform = boundedStdUniformCDF+ instance CDF StdUniform Int where cdf ~StdUniform = boundedStdUniformCDF |]) @@ -187,11 +204,20 @@ instance CDF Uniform Float where cdf (Uniform a b) = realUniformCDF a b instance CDF Uniform Double where cdf (Uniform a b) = realUniformCDF a b -instance Distribution StdUniform Float where rvar StdUniform = floatStdUniform-instance Distribution StdUniform Double where rvar StdUniform = doubleStdUniform-instance CDF StdUniform Float where cdf StdUniform = realStdUniformCDF-instance CDF StdUniform Double where cdf StdUniform = realStdUniformCDF+instance Distribution StdUniform Float where rvar ~StdUniform = floatStdUniform+instance Distribution StdUniform Double where rvar ~StdUniform = doubleStdUniform+instance CDF StdUniform Float where cdf ~StdUniform = realStdUniformCDF+instance CDF StdUniform Double where cdf ~StdUniform = realStdUniformCDF +instance HasResolution r => + Distribution Uniform (Fixed r) where rvar (Uniform a b) = fixedUniform a b+instance HasResolution r => + CDF Uniform (Fixed r) where cdf (Uniform a b) = realUniformCDF a b+instance HasResolution r =>+ Distribution StdUniform (Fixed r) where rvar ~StdUniform = fixedStdUniform+instance HasResolution r => + CDF StdUniform (Fixed r) where cdf ~StdUniform = realStdUniformCDF+ instance Distribution Uniform () where rvar (Uniform a b) = return () instance CDF Uniform () where cdf (Uniform a b) = return 1 $( replicateInstances ''Char [''Char, ''Bool, ''Ordering] [d|@@ -200,13 +226,13 @@ |]) -instance Distribution StdUniform () where rvar StdUniform = return ()-instance CDF StdUniform () where cdf StdUniform = return 1-instance Distribution StdUniform Bool where rvar StdUniform = fmap even getRandomByte-instance CDF StdUniform Bool where cdf StdUniform = boundedEnumStdUniformCDF+instance Distribution StdUniform () where rvar ~StdUniform = return ()+instance CDF StdUniform () where cdf ~StdUniform = return 1+instance Distribution StdUniform Bool where rvar ~StdUniform = fmap even getRandomByte+instance CDF StdUniform Bool where cdf ~StdUniform = boundedEnumStdUniformCDF -instance Distribution StdUniform Char where rvar StdUniform = boundedEnumStdUniform-instance CDF StdUniform Char where cdf StdUniform = boundedEnumStdUniformCDF-instance Distribution StdUniform Ordering where rvar StdUniform = boundedEnumStdUniform-instance CDF StdUniform Ordering where cdf StdUniform = boundedEnumStdUniformCDF+instance Distribution StdUniform Char where rvar ~StdUniform = boundedEnumStdUniform+instance CDF StdUniform Char where cdf ~StdUniform = boundedEnumStdUniformCDF+instance Distribution StdUniform Ordering where rvar ~StdUniform = boundedEnumStdUniform+instance CDF StdUniform Ordering where cdf ~StdUniform = boundedEnumStdUniformCDF
+ src/Data/Random/Internal/Fixed.hs view
@@ -0,0 +1,22 @@+module Data.Random.Internal.Fixed where++import Data.Fixed+import Unsafe.Coerce++resolutionOf :: HasResolution r => f r -> Integer+resolutionOf x = resolution (res x)+ where+ res :: HasResolution r => f r -> r+ res = undefined++resolutionOf2 :: HasResolution r => f (g r) -> Integer+resolutionOf2 x = resolution (res x)+ where+ res :: HasResolution r => f (g r) -> r+ res = undefined++mkFixed :: Integer -> Fixed r+mkFixed = unsafeCoerce++unMkFixed :: Fixed r -> Integer+unMkFixed = unsafeCoerce
src/Data/Random/Source/PureMT.hs view
@@ -29,14 +29,14 @@ -- 'RandomSource' instances for mutable references to 'PureMT' states. getRandomWordFromMTRef :: ModifyRef sr m PureMT => sr -> m Word64 getRandomWordFromMTRef ref = do- atomicModifyRef ref (swap . randomWord64)+ atomicModifyReference ref (swap . randomWord64) where swap (a,b) = (b,a) -getRandomByteFromMTRef :: ModifyRef sr m PureMT => sr -> m Word8+getRandomByteFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> m Word8 getRandomByteFromMTRef ref = do- x <- atomicModifyRef ref (swap . randomInt)+ x <- atomicModifyReference ref (swap . randomInt) return (fromIntegral x) where@@ -46,10 +46,10 @@ -- the mersenne random library is using, at least in the version I have. -- if this changes, switch to the commented version. -- Same thing below, in getRandomDoubleFromMTState.-getRandomDoubleFromMTRef :: ModifyRef sr m PureMT => sr -> m Double+getRandomDoubleFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> m Double getRandomDoubleFromMTRef src = liftM wordToDouble (getRandomWordFromMTRef src) -- getRandomDoubleFromMTRef ref = do--- atomicModifyRef ref (swap . randomDouble)+-- atomicModifyReference ref (swap . randomDouble) -- -- where -- swap (a,b) = (b,a)@@ -95,6 +95,11 @@ getRandomWord = getRandomWordFromMTState getRandomDouble = getRandomDoubleFromMTState +instance (Monad m1, ModifyRef (Ref m2 PureMT) m1 PureMT) => RandomSource m1 (Ref m2 PureMT) where+ getRandomByteFrom = getRandomByteFromMTRef+ getRandomWordFrom = getRandomWordFromMTRef+ getRandomDoubleFrom = getRandomDoubleFromMTRef+ instance Monad m => MonadRandom (StateT PureMT m) where getRandomByte = getRandomByteFromMTState getRandomWord = getRandomWordFromMTState@@ -105,21 +110,20 @@ getRandomWord = getRandomWordFromMTState getRandomDouble = getRandomDoubleFromMTState --instance (ModifyRef (IORef PureMT) m PureMT) => RandomSource m (IORef PureMT) where+instance (Monad m, ModifyRef (IORef PureMT) m PureMT) => RandomSource m (IORef PureMT) where {-# SPECIALIZE instance RandomSource IO (IORef PureMT)#-} getRandomByteFrom = getRandomByteFromMTRef getRandomWordFrom = getRandomWordFromMTRef getRandomDoubleFrom = getRandomDoubleFromMTRef -instance (ModifyRef (STRef s PureMT) m PureMT) => RandomSource m (STRef s PureMT) where+instance (Monad m, ModifyRef (STRef s PureMT) m PureMT) => RandomSource m (STRef s PureMT) where {-# SPECIALIZE instance RandomSource (ST s) (STRef s PureMT) #-} {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s PureMT) #-} getRandomByteFrom = getRandomByteFromMTRef getRandomWordFrom = getRandomWordFromMTRef getRandomDoubleFrom = getRandomDoubleFromMTRef -instance (ModifyRef (TVar PureMT) m PureMT) => RandomSource m (TVar PureMT) where+instance (Monad m, ModifyRef (TVar PureMT) m PureMT) => RandomSource m (TVar PureMT) where {-# SPECIALIZE instance RandomSource IO (TVar PureMT) #-} {-# SPECIALIZE instance RandomSource STM (TVar PureMT) #-} getRandomByteFrom = getRandomByteFromMTRef
src/Data/Random/Source/StdGen.hs view
@@ -16,18 +16,23 @@ import Data.StateRef import Data.Word -instance (ModifyRef (IORef StdGen) m StdGen) => RandomSource m (IORef StdGen) where+instance (Monad m1, ModifyRef (Ref m2 StdGen) m1 StdGen) => RandomSource m1 (Ref m2 StdGen) where+ getRandomByteFrom = getRandomByteFromRandomGenRef+ getRandomWordFrom = getRandomWordFromRandomGenRef+ getRandomDoubleFrom = getRandomDoubleFromRandomGenRef++instance (Monad m, ModifyRef (IORef StdGen) m StdGen) => RandomSource m (IORef StdGen) where {-# SPECIALIZE instance RandomSource IO (IORef StdGen) #-} getRandomByteFrom = getRandomByteFromRandomGenRef getRandomWordFrom = getRandomWordFromRandomGenRef getRandomDoubleFrom = getRandomDoubleFromRandomGenRef-instance (ModifyRef (TVar StdGen) m StdGen) => RandomSource m (TVar StdGen) where+instance (Monad m, ModifyRef (TVar StdGen) m StdGen) => RandomSource m (TVar StdGen) where {-# SPECIALIZE instance RandomSource IO (TVar StdGen) #-} {-# SPECIALIZE instance RandomSource STM (TVar StdGen) #-} getRandomByteFrom = getRandomByteFromRandomGenRef getRandomWordFrom = getRandomWordFromRandomGenRef getRandomDoubleFrom = getRandomDoubleFromRandomGenRef-instance (ModifyRef (STRef s StdGen) m StdGen) => RandomSource m (STRef s StdGen) where+instance (Monad m, ModifyRef (STRef s StdGen) m StdGen) => RandomSource m (STRef s StdGen) where {-# SPECIALIZE instance RandomSource (ST s) (STRef s StdGen) #-} {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s StdGen) #-} getRandomByteFrom = getRandomByteFromRandomGenRef@@ -63,19 +68,19 @@ -- away a lot of perfectly good entropy. Better still is to use these 3 functions -- together to create a 'RandomSource' instance for the reference you're using, -- if one does not already exist.-getRandomByteFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>+getRandomByteFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) => sr -> m Word8-getRandomByteFromRandomGenRef g = atomicModifyRef g (swap . randomR (0,255))+getRandomByteFromRandomGenRef g = atomicModifyReference g (swap . randomR (0,255)) where swap :: (Int, a) -> (a, Word8) swap (a,b) = (b,fromIntegral a) -getRandomWordFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>+getRandomWordFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) => sr -> m Word64-getRandomWordFromRandomGenRef g = atomicModifyRef g (swap . randomR (0,0xffffffffffffffff))+getRandomWordFromRandomGenRef g = atomicModifyReference g (swap . randomR (0,0xffffffffffffffff)) where swap (a,b) = (b,fromInteger a) -getRandomDoubleFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>+getRandomDoubleFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) => sr -> m Double getRandomDoubleFromRandomGenRef g = liftM wordToDouble (getRandomWordFromRandomGenRef g) -- getRandomDoubleFromRandomGenRef g = atomicModifyRef g (swap . randomR (0,1))