diff --git a/random-fu.cabal b/random-fu.cabal --- a/random-fu.cabal +++ b/random-fu.cabal @@ -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 diff --git a/src/Data/Random/Distribution/Discrete.hs b/src/Data/Random/Distribution/Discrete.hs --- a/src/Data/Random/Distribution/Discrete.hs +++ b/src/Data/Random/Distribution/Discrete.hs @@ -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) diff --git a/src/Data/Random/Distribution/Normal.hs b/src/Data/Random/Distribution/Normal.hs --- a/src/Data/Random/Distribution/Normal.hs +++ b/src/Data/Random/Distribution/Normal.hs @@ -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 diff --git a/src/Data/Random/Distribution/Uniform.hs b/src/Data/Random/Distribution/Uniform.hs --- a/src/Data/Random/Distribution/Uniform.hs +++ b/src/Data/Random/Distribution/Uniform.hs @@ -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 diff --git a/src/Data/Random/Internal/Fixed.hs b/src/Data/Random/Internal/Fixed.hs new file mode 100644 --- /dev/null +++ b/src/Data/Random/Internal/Fixed.hs @@ -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 diff --git a/src/Data/Random/Source/PureMT.hs b/src/Data/Random/Source/PureMT.hs --- a/src/Data/Random/Source/PureMT.hs +++ b/src/Data/Random/Source/PureMT.hs @@ -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 diff --git a/src/Data/Random/Source/StdGen.hs b/src/Data/Random/Source/StdGen.hs --- a/src/Data/Random/Source/StdGen.hs +++ b/src/Data/Random/Source/StdGen.hs @@ -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))