diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) Patrick Perry 2008 +Copyright (c) Patrick Perry 2010 All rights reserved. diff --git a/NEWS b/NEWS new file mode 100644 --- /dev/null +++ b/NEWS @@ -0,0 +1,40 @@ + +Changes in 0.3: + +* Add strict versions of sampleSubset, sampleIntSubset, and shuffleInt. + +* Port to vector-0.6.0. + +* Add Exponential and Levy alpha-Stable distributions. + +* Add Summary.Bool for indicators. + +* Move Summary to Data.Summary + +* Introduce `repeatMC`, which produces an infinite (lazy) stream of values, and + `replicateMC`, which produces a lazy list of specified length. + +* Remove `repeatMC/repeatMCWith`. + +* Build fix for 6.8.2 from Robert Gunst. + +* The function `sample`, `sampleWithWeights`, `sampleSubset`, and + `shuffle` no longer require that you explicitly pass in the length. + +* The pure RNG is now a newtype, so you can't use the functions from + GLS.Random.Gen on it anymore. + +* The internals of the monad have been cleaned up. IO is used internally + instead of `seq` calls and `unsafePerformIO` everywhere. This results in + a modest performance boost. + + +Changes in 0.2: + +* More general type class, MonadMC, which allows all the functions to work + in both MC and MCT monads. + +* Functions to sample from discrete distributions. + +* Functions to sample subsets + diff --git a/examples/Binomial.hs b/examples/Binomial.hs new file mode 100644 --- /dev/null +++ b/examples/Binomial.hs @@ -0,0 +1,55 @@ +module Main + where + +import Control.Monad +import Data.List( foldl' ) +import Text.Printf( printf ) + +import Control.Monad.MC +import Data.Summary +import Data.Summary.Utils( inInterval ) + +-- | Sample from a binomial distribution with the given parameters. +binomial :: (MonadMC m) => Int -> Double -> m Int +binomial n p = let + q = 1 - p + probs = map (\i -> (fromIntegral $ n `choose` i) * p^^i * q^^(n-i)) [0..n] + in sampleIntWithWeights probs (n+1) + +-- | Get a sample confidence interval for the mean after @reps@ replications of +-- a binomial with the given parameters. +binomialMean :: (MonadMC m) => Int -> Double -> Int -> m (Double,Double) +binomialMean n p reps = + liftM (sampleCI 0.95 . summary . map fromIntegral) $ + replicateMC reps (binomial n p) + +-- | Compute @reps@ 95% confidence intervals for the mean of an @(n,p)@ +-- binormal based on samples of the given size, and record the number +-- of intervals that contain the true mean. +coverage :: (MonadMC m) => Int -> Double -> Int -> Int -> m Int +coverage n p size reps = + liftM (length . filter (mu `inInterval`)) $ + replicateMC reps $ + binomialMean n p size + where + mu = fromIntegral n * p + +main = + let seed = 0 + reps = 100 + n = 10 + p = 0.2 + size = 500 + c = coverage n p size reps `evalMC` mt19937 seed + in + printf "\nOf %d 95%%-intervals, %d contain the true value.\n" reps c + + +--------------------------- Utility functions ----------------------------- + +factorial :: Int -> Int +factorial n | n <= 0 = 1 + | otherwise = n * factorial (n-1) + +choose :: Int -> Int -> Int +choose n k = factorial n `div` (factorial (n-k) * factorial k) diff --git a/examples/Pi.hs b/examples/Pi.hs deleted file mode 100644 --- a/examples/Pi.hs +++ /dev/null @@ -1,82 +0,0 @@ - -import Control.Monad.MC -import Control.Monad -import Data.List( foldl' ) -import System.Environment( getArgs ) -import Text.Printf( printf ) - --- | Generate a point in the box [-1,1) x [-1,1) -unitBox :: MC (Double,Double) -unitBox = liftM2 (,) (uniform (-1) 1) - (uniform (-1) 1) - --- | Indicates whether or not a point is in the unit circle -inUnitCircle :: (Double,Double) -> Bool -inUnitCircle (x,y) = x*x + y*y <= 1 - --- | Given a list of indicators, return the sample mean and standard --- error. -average :: [Bool] -> (Double,Double) -average is = let - (t,n) = foldl' count (0,0) is - p = toDouble t / toDouble n - se = sqrt (p * (1 - p) / toDouble n) - in (p, se) - where - count (t,n) i = let - t' = if i then t+1 else t - n' = n+1 - in t' `seq` n' `seq` (t',n') - - toDouble = realToFrac . toInteger - --- | Compute a Monte Carlo estimate of pi based on @n@ samples. Return --- the estimate and the standard error of the estimate. -computePi :: Int -> MC (Double,Double) -computePi n = do - is <- liftM (map inUnitCircle) (unsafeInterleaveMC $ replicateM n unitBox) - let (mu ,se ) = average is - (mu',se') = (4*mu,4*se) - return (mu',se') - --- | Given an estimate and standard error, produce a 99% confidence --- interval based on the Central Limit Theorem -interval :: Double -> Double -> (Double,Double) -interval mu se = let - delta = 2.575*se - in (mu-delta, mu+delta) - --- | Tests if the value is in the interval [a,b] -inInterval :: Double -> (Double,Double) -> Bool -x `inInterval` (a,b) = x >= a && x <= b - --- | Compute an estimate of pi based on @n@ points and see if the true --- value is in the confidence interval -covers :: Int -> MC Bool -covers n = do - (mu,se) <- computePi n - return $ pi `inInterval` (interval mu se) - --- | Compute @r@ estimates of pi based on @n@ samples each, and count --- how many times the true values is included in the 99% confidence --- inverval -coverage :: Int -> Int -> MC Int -coverage r n = do - liftM count $ replicateM r (covers n) - where - count = length . filter id - -main = do - [n] <- map read `fmap` getArgs - main' n - -main' n = let - seed = 0 - (mu,se) = evalMC (computePi n) $ mt19937 seed - (l,u) = interval mu se - r = 500 - c = evalMC (coverage r n) $ mt19937 seed - in do - printf "Estimate from one simulation: %g\n" mu - printf "99%% Confidence Interval: (%g,%g)\n" l u - printf "\nOf %d intervals, %d contain the true value.\n" r c diff --git a/examples/Pi.lhs b/examples/Pi.lhs new file mode 100644 --- /dev/null +++ b/examples/Pi.lhs @@ -0,0 +1,113 @@ + +In this example, we compute a Monte Carlo estimate of pi by +generating random points in the unit box, and counting how many +of them fall in the unit circle. + +\begin{code} +import Control.Monad( liftM, liftM2 ) +import Control.Monad.MC( MC, uniform, replicateMC, evalMC, mt19937 ) +import Data.Summary.Bool( summary, sampleMean, sampleSE ) +import Data.Summary.Utils( interval ) +import Text.Printf +\end{code} + +First, we need a function to test whether or not a point is in the +unit circle. We define + +\begin{code} +inUnitCircle :: (Double,Double) -> Bool +inUnitCircle (x,y) = x*x + y*y <= 1 +\end{code} + +The first line is the type signature, which tells us that "inUnitCircle" +is a function which takes a pair of `Double`s and returns a `Bool`. In +English, "::" means "has type". The second line is the one that defines +the function. + +\begin{code} +estimatePi :: [(Double,Double)] -> (Double,Double) +estimatePi xs = + let s = summary $ map inUnitCircle xs + (mu,se) = (sampleMean s, sampleSE s) in + (4*mu,4*se) +\end{code} + +Next, we need to generate a random point in the unit box. In the +Control.Monad.MC module, there is a function for generating uniform +values in an interval, called "uniform". This function has a +funny-looking type, but you can think of it as: + + uniform :: Double -> Double -> MC Double + +This type means that the function takes the two endpoints of the +interval as arguments, and returns a Monte-Carlo action which produces +a Double. + +You can think of the type "MC Double" as a random number generator. +For general types, "MC a" is a generator for values of type "a". In +fact, "MC" is one of a general class of objects called a Monad. We +use Monads in Haskell for making sure events happen in the right order. +That is, if we have three parts of a simulation, say A, B, and C, and +we want them to happen in the order + + A ====> B ====> C + +then we would like to make sure that A is done consuming random numbers +before B consumes anything. Likewise, we want B to finish before C +starts. Monads are the magic that enable us to ensure this. + +There are a number of functions in the standard library for working with +monads. The first we will use is + + liftM2 :: (Monad m) => (a1 -> a2 -> r) -> m a1 -> m a2 -> m r + +This function works on *any* Monad. When we use it on the MC monad, it +will have type + + liftM2 :: (a1 -> a2 -> r) -> MC a1 -> MC a2 -> MC r + +What liftM2 does is it takes a function of two arguments and two Monte +Carlo actions. It returns a new Monte Carlo action that does the following: + + 1. generate a random value of type a1 using the first action + 2. generate a random value of type a2 using the second action + 3. apply a function to the two values and return the result + +We do not need to write liftM2 ourselves, since it is provided in the +"Control.Monad" module. But, if we did have to define it, the code would +look like: + +liftM2 f ma1 ma2 = do + a1 <- ma1 + a2 <- ma2 + return (f a1 a2) + +This code for this uses the "do" notation of Haskell, which allows us +to specify a series of actions in sequential order. + +\begin{code} +unitBox :: MC (Double,Double) +unitBox = liftM2 (,) (uniform (-1) 1) + (uniform (-1) 1) +\end{code} + +-- | Compute a Monte Carlo estimate of pi based on @n@ samples. Return +-- the sample mean and standard error. + +\begin{code} +simulation :: Int -> MC (Double,Double) +simulation n = + estimatePi `fmap` replicateMC n unitBox +\end{code} + +\begin{code} +main = + let seed = 0 + n = 1000000 + (mu,se) = simulation n `evalMC` mt19937 seed + (l,u) = interval 0.95 mu se + in do + printf "\nEstimate: %g" mu + printf "\n99%% Confidence Interval: (%g, %g)" l u + printf "\n" +\end{code} diff --git a/examples/Poker.hs b/examples/Poker.hs --- a/examples/Poker.hs +++ b/examples/Poker.hs @@ -5,7 +5,6 @@ import Data.List import Data.Map( Map ) import qualified Data.Map as Map -import System.Environment import Text.Printf -- | Data types for representing cards. An Ace has 'number' equal to @1@. @@ -23,6 +22,12 @@ queen = 12 king = 13 +-- | Get a list of cards that make up a 52-card deck. +deck :: [Card] +deck = [ Card i s + | i <- [ ace..king ] + , s <- [ Club, Diamond, Heart, Spade ] ] + -- | A type for the various poker hands. data Hand = HighCard | Pair | TwoPair | ThreeOfAKind | Straight | Flush | FullHouse | FourOfAKind | StraightFlush @@ -54,13 +59,9 @@ matches = (sort . map length . group) (x:xs) --- | Get a list of cards that make up a 52-card deck. -deck :: [Card] -deck = [ Card i s | i <- [ 1..13 ], s <- [ Club, Diamond, Heart, Spade ] ] - -- | Deal a five-card hand by choosing a random subset of the deck. deal :: (MonadMC m) => m [Card] -deal = sampleSubset 5 52 deck +deal = sampleSubset 5 deck -- | A type for storing the frequencies of the various hands. type HandCounts = Map Hand Int @@ -74,22 +75,20 @@ updateCounts counts cs = Map.insertWith' (+) (hand cs) 1 counts -main = do - [reps] <- map read `fmap` getArgs - main' reps - -main' reps = +main = let seed = 0 - counts = repeatMCWith updateCounts emptyCounts reps deal - `evalMC` mt19937 seed in do - printf "\n" - printf " Hand Count Probability 99%% Interval \n" - printf "-------------------------------------------------------\n" - forM_ ((reverse . Map.toAscList) counts) $ \(h,c) -> - let n = fromIntegral reps :: Double - p = fromIntegral c / n - se = sqrt (p * (1 - p) / n) - delta = 2.575829 * se - (l,u) = (p-delta, p+delta) in - printf "%-13s %7d %.6f (%.6f,%.6f)\n" (show h) c p l u - printf "\n" + reps = 100000 + counts = foldl' updateCounts emptyCounts $ + replicateMC reps deal `evalMC` mt19937 seed + in do + printf "\n" + printf " Hand Count Probability 99%% Interval \n" + printf "-------------------------------------------------------\n" + forM_ ((reverse . Map.toAscList) counts) $ \(h,c) -> + let n = fromIntegral reps :: Double + p = fromIntegral c / n + se = sqrt (p * (1 - p) / n) + delta = 2.575829 * se + (l,u) = (p-delta, p+delta) in + printf "%-13s %7d %.6f (%.6f,%.6f)\n" (show h) c p l u + printf "\n" diff --git a/examples/Queue.hs b/examples/Queue.hs new file mode 100644 --- /dev/null +++ b/examples/Queue.hs @@ -0,0 +1,192 @@ + +import Control.Monad +import Control.Monad.MC +import Data.List( foldl' ) +import Data.Summary +import Text.Printf( printf ) + +-- | There a three items on the menu. +data Item = Cheeseburger | Fries | Milkshake + +-- | A customer orders some number of items +data Customer = Customer { orderOf :: [Item] } + +-- | The order size is a Poisson random variable with mean 2. +orderSize :: MC Int +orderSize = liftM (1+) $ poisson 2 + +-- | The items are sampled with the given weights. +item :: MC Item +item = sampleWithWeights [ (4, Cheeseburger), (2, Fries), (1, Milkshake) ] + +-- | Generate a random order. +order :: MC [Item] +order = do + n <- orderSize + replicateM n item + +-- | Generate a random customer. +customer :: MC Customer +customer = liftM Customer order + +-- | A customer event. The interarrival time is the time that elapeses +-- between when the previous customer arrives and when the current customer +-- arrives. +data CustomerEvent = CustomerEvent { customerOf :: !Customer + , interarrivalTime :: !Double + } + +-- | Generate a random customer event. The interarrival time distribution +-- is exponential with mean 1. +customerEvent :: MC CustomerEvent +customerEvent = do + c <- customer + delta <- exponential 10 + return $ CustomerEvent c delta + +-- | The time it takes to make an item. +cook :: Item -> MC Double +cook Cheeseburger = exponential 3 +cook Fries = exponential 1 +cook Milkshake = exponential 2 + +-- | The time it takes to cook all of the items in the list is equal +-- to the maximum time. +cookAll :: [Item] -> MC Double +cookAll items = do + ts <- mapM cook items + return $ foldl' max 0 ts + +-- | A customer in line, along with how long they have been waiting. +data Waiting = Waiting { waiting :: !Customer + , hasBeenWaiting :: !Double + } + +-- | A customer, along with how long it takes to prepare the customer's order +-- and how long the customer has to wait. +data Service = Service { serving :: !Customer + , waitingTime :: !Double + , serviceTime :: !Double + } + +-- | Given a customer who has been wating in line, provide them with service. +-- If the customer has been waiting for longer than 5 minutes, work twice as +-- fast to cook the food. +serveWaiting :: Waiting -> MC Service +serveWaiting (Waiting c w) = do + t <- cookAll $ orderOf c + let t' = if w > 5 then 0.5*t else t + return $ Service c w t' + +-- | A resturant has one server, who may be busy. There is a list of +-- customers wating in line. +data Restaurant = Restaurant { inProgress :: Maybe InProgress + , waitingLine :: [Waiting] + } + +-- | An in-progress service event. +data InProgress = InProgress { service :: !Service + , timeToFinish :: !Double + } + +-- | Update the amount of time the customers have been waiting by adding +-- the given amount. +addToWait :: Double -> [Waiting] -> [Waiting] +addToWait delta = map (\(Waiting w t) -> Waiting w (t+delta)) + +-- | Serve customers in the restaurant for the given amount of time. +serveForTime :: Double -> Restaurant -> MC ([Service], Restaurant) +serveForTime = + let serveForTimeHelp ss t r = case r of + -- When no one is being served and no one is in line, do nothing. + Restaurant Nothing [] -> + return $ (ss, r) + + -- When no one is being served, take the first person in line + -- and start cooking their order. + Restaurant Nothing (x:xs) -> do + s <- serveWaiting x + let y = Just $ InProgress s $ serviceTime s + serveForTimeHelp ss t $ Restaurant y xs + + -- When somone is being served, serve them for the given amount + -- of time. If we have enough time, finish serving them and + -- update the amount of time everyone else has had to wait. + -- Otherwise, just update the time to finish serving and + -- update the waiting times of the customers in line. + Restaurant (Just (InProgress s delta)) xs -> + if delta <= t then let t' = t - delta + xs' = addToWait delta xs + r' = Restaurant Nothing xs' in + serveForTimeHelp (ss ++ [s]) t' r' + else let delta' = delta - t + y' = Just $ InProgress s delta' + xs' = addToWait t xs + r' = Restaurant y' xs' in + return (ss,r') + in serveForTimeHelp [] + +-- | Given a new customer arrival event, produce a list of all of the new +-- service events that happen before the customer gets there, and return +-- the updated restaurant state at the time immediately after the customer +-- arrives. +processEvent :: CustomerEvent + -> Restaurant + -> MC ([Service], Restaurant) +processEvent (CustomerEvent c t) r = do + (ss,(Restaurant y xs)) <- serveForTime t r + return $ (ss, (Restaurant y $ xs ++ [Waiting c 0])) + +-- | Finish serving all of the customers in line. +finishServing :: Restaurant -> MC [Service] +finishServing r = do + (ss,_) <- serveForTime infinity r + return ss + where + infinity = 1/0 + +-- | A restaurant takes a list of customer events and generates a random +-- list of service events. The reason for the call to "unsafeInterleaveMC" +-- is that we want to make sure that we return a lazy list. Without it, +-- the function will return only after it has consumed all of the random +-- numbers it needs. This is problemeatic if the input list is large or +-- or infinite. +restaurant :: [CustomerEvent] -> MC [Service] +restaurant = + let restaurantHelp r [] = finishServing r + restaurantHelp r (c:cs) = unsafeInterleaveMC $ do + (ss,r') <- processEvent c r + ss' <- restaurantHelp r' cs + return $ ss ++ ss' + in restaurantHelp (Restaurant Nothing []) + +-- | An infinite stream of customerEvents. This stream uses its own private +-- random number generator (mt19937 is the Mersenne-Twister algorithm). +customerEvents :: Seed -> [CustomerEvent] +customerEvents seed = repeatMC customerEvent `evalMC` mt19937 seed + +-- | Given a seed for the customers and a seed for the restaurant, run the +-- simulation. +simulation :: Seed -> Seed -> [Service] +simulation customerSeed restaurantSeed = + restaurant (customerEvents customerSeed) `evalMC` mt19937 restaurantSeed + +-- | Compute a summary of the total waitings time for each customer. +summarize :: [Service] -> Summary +summarize = summary . map totalTime + where + totalTime (Service _ w s) = w+s + +-- | Run the program +main = + let customerSeed = 0 + restaurantSeed = 100 + numTransactions = 100000 + results = summarize $ take numTransactions $ + simulation customerSeed restaurantSeed + in do + putStrLn "" + putStrLn "Total Service Time:" + putStrLn "-------------------" + putStrLn $ show $ results + putStrLn "" diff --git a/examples/Sampling.hs b/examples/Sampling.hs deleted file mode 100644 --- a/examples/Sampling.hs +++ /dev/null @@ -1,58 +0,0 @@ -module Main - where - -import Control.Monad.MC -import Control.Monad -import Data.List( foldl' ) -import System.Environment( getArgs ) -import Text.Printf( printf ) - - --- | Sample from a binomial distribution with the given parameters. -binomial :: (MonadMC m) => Int -> Double -> m Int -binomial n p = let - q = 1 - p - probs = map (\i -> (fromIntegral $ n `choose` i) * p^^i * q^^(n-i)) [0..n] - in sampleIntWithWeights probs (n+1) - --- | Get a sample confidence interval for the mean after @reps@ replications of --- a binomial with the given parameters. -binomialMean :: (MonadMC m) => Int -> Double -> Int -> m (Double,Double) -binomialMean n p reps = - liftM (sampleCI 0.95) $ repeatMC reps $ liftM fromIntegral (binomial n p) - --- | Compute @reps@ 95% confidence intervals for the mean of an @(n,p)@ --- binormal based on samples of the given size, and record the number --- of intervals that contain the true mean. -coverage :: (MonadMC m) => Int -> Double -> Int -> Int -> m Int -coverage n p size reps = - repeatMCWith - (\c ci -> if mu `inInterval` ci then c+1 else c) - 0 - reps - (binomialMean n p size) - where - mu = fromIntegral n * p - x `inInterval` (l,h) = x > l && x < h - -main = do - [reps] <- map read `fmap` getArgs - main' reps - -main' reps = - let seed = 0 - n = 10 - p = 0.2 - size = 500 - c = evalMC (coverage n p size reps) $ mt19937 seed in - printf "\nOf %d 95%%-intervals, %d contain the true value.\n" reps c - - ---------------------------- Utility functions ----------------------------- - -factorial :: Int -> Int -factorial n | n <= 0 = 1 - | otherwise = n * factorial (n-1) - -choose :: Int -> Int -> Int -choose n k = factorial n `div` (factorial (n-k) * factorial k) diff --git a/lib/Control/Monad/MC.hs b/lib/Control/Monad/MC.hs --- a/lib/Control/Monad/MC.hs +++ b/lib/Control/Monad/MC.hs @@ -1,10 +1,14 @@ ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC --- Copyright : Copyright (c) , Patrick Perry +-- Copyright : Copyright (c) 2010, Patrick Perry -- License : BSD3 --- Maintainer : Patrick Perry +-- Maintainer : Patrick Perry -- Stability : experimental +-- +-- A monad and monad transformer for monte carlo computations. Currently, +-- the default is the GNU Scientific Library-based implementation, but this +-- may change in the future. -- module Control.Monad.MC ( diff --git a/lib/Control/Monad/MC/Base.hs b/lib/Control/Monad/MC/Base.hs --- a/lib/Control/Monad/MC/Base.hs +++ b/lib/Control/Monad/MC/Base.hs @@ -1,20 +1,17 @@ -{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, FlexibleContexts #-} +{-# LANGUAGE TypeFamilies #-} ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Base --- Copyright : Copyright (c) , Patrick Perry +-- Copyright : Copyright (c) 2010, Patrick Perry -- License : BSD3 --- Maintainer : Patrick Perry +-- Maintainer : Patrick Perry -- Stability : experimental -- -module Control.Monad.MC.Base ( - -- * MonadMC type classes - MonadMC(..), - HasRNG(..), - - ) where +module Control.Monad.MC.Base + where +import Control.Monad import qualified Control.Monad.MC.GSLBase as GSL class HasRNG m where @@ -38,6 +35,19 @@ -- | @normal mu sigma@ generates a Normal random variable with mean -- @mu@ and standard deviation @sigma@. normal :: Double -> Double -> m Double + + -- | @exponential mu@ generates an Exponential variate with mean @mu@. + exponential :: Double -> m Double + + -- | @levy c alpha@ gets a Levy alpha-stable variate with scale @c@ and + -- exponent @alpha@. The algorithm only works for @0 < alpha <= 2@. + levy :: Double -> Double -> m Double + + -- | @levySkew c alpha beta @ gets a skew Levy alpha-stable variate + -- with scale @c@, exponent @alpha@, and skewness @beta@. The skew + -- parameter must lie in the range @[-1,1]@. The algorithm only works + -- for @0 < alpha <= 2@. + levySkew :: Double -> Double -> Double -> m Double -- | @poisson mu@ generates a Poisson random variable with mean @mu@. poisson :: Double -> m Int @@ -47,6 +57,11 @@ unsafeInterleaveMC :: m a -> m a +-- | Generate 'True' events with the given probability +bernoulli :: (MonadMC m) => Double -> m Bool +bernoulli p = liftM (< p) $ uniform 0 1 +{-# INLINE bernoulli #-} + ------------------------------- Instances ----------------------------------- instance HasRNG GSL.MC where @@ -63,6 +78,12 @@ {-# INLINE uniformInt #-} normal = GSL.normal {-# INLINE normal #-} + exponential = GSL.exponential + {-# INLINE exponential #-} + levy = GSL.levy + {-# INLINE levy #-} + levySkew = GSL.levySkew + {-# INLINE levySkew #-} poisson = GSL.poisson {-# INLINE poisson #-} unsafeInterleaveMC = GSL.unsafeInterleaveMC @@ -82,6 +103,12 @@ {-# INLINE uniformInt #-} normal mu sigma = GSL.liftMCT $ GSL.normal mu sigma {-# INLINE normal #-} + exponential mu = GSL.liftMCT $ GSL.exponential mu + {-# INLINE exponential #-} + levy c alpha = GSL.liftMCT $ GSL.levy c alpha + {-# INLINE levy #-} + levySkew c alpha beta = GSL.liftMCT $ GSL.levySkew c alpha beta + {-# INLINE levySkew #-} poisson mu = GSL.liftMCT $ GSL.poisson mu {-# INLINE poisson #-} unsafeInterleaveMC = GSL.unsafeInterleaveMCT diff --git a/lib/Control/Monad/MC/Class.hs b/lib/Control/Monad/MC/Class.hs --- a/lib/Control/Monad/MC/Class.hs +++ b/lib/Control/Monad/MC/Class.hs @@ -1,32 +1,25 @@ ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Class --- Copyright : Copyright (c) , Patrick Perry +-- Copyright : Copyright (c) 2010, Patrick Perry -- License : BSD3 --- Maintainer : Patrick Perry +-- Maintainer : Patrick Perry -- Stability : experimental -- +-- The abstract MonadMC interface and utility functions for Monte Carlo +-- computations. +-- module Control.Monad.MC.Class ( -- * The Monte Carlo monad type class HasRNG(..), - MonadMC, - - -- * Getting and setting the generator - getRNG, - setRNG, + MonadMC(..), -- * Random distributions - uniform, - uniformInt, - normal, - poisson, + bernoulli, module Control.Monad.MC.Sample, module Control.Monad.MC.Repeat, - - -- * Interleaving computations - unsafeInterleaveMC ) where import Control.Monad.MC.Base diff --git a/lib/Control/Monad/MC/GSL.hs b/lib/Control/Monad/MC/GSL.hs --- a/lib/Control/Monad/MC/GSL.hs +++ b/lib/Control/Monad/MC/GSL.hs @@ -1,11 +1,13 @@ ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.GSL --- Copyright : Copyright (c) , Patrick Perry +-- Copyright : Copyright (c) 2010, Patrick Perry -- License : BSD3 --- Maintainer : Patrick Perry +-- Maintainer : Patrick Perry -- Stability : experimental -- +-- A monad and monad transformer for monte carlo computations built on top +-- of the functions in the GNU Scientific Library. module Control.Monad.MC.GSL ( -- * The Monte Carlo monad @@ -22,7 +24,12 @@ -- * Pure random number generator creation RNG, + Seed, mt19937, + mt19937WithState, + rngName, + rngSize, + rngState, -- * Overloaded Monte Carlo monad interface module Control.Monad.MC.Class, @@ -30,5 +37,6 @@ ) where import Control.Monad.MC.GSLBase ( MC, runMC, evalMC, execMC, - MCT, runMCT, evalMCT, execMCT, RNG, mt19937 ) + MCT, runMCT, evalMCT, execMCT, RNG, Seed, mt19937, mt19937WithState, + rngName, rngSize, rngState ) import Control.Monad.MC.Class hiding ( RNG ) diff --git a/lib/Control/Monad/MC/GSLBase.hs b/lib/Control/Monad/MC/GSLBase.hs --- a/lib/Control/Monad/MC/GSLBase.hs +++ b/lib/Control/Monad/MC/GSLBase.hs @@ -1,10 +1,10 @@ -{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-} +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-} ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.GSLBase --- Copyright : Copyright (c) , Patrick Perry +-- Copyright : Copyright (c) 2010, Patrick Perry -- License : BSD3 --- Maintainer : Patrick Perry +-- Maintainer : Patrick Perry -- Stability : experimental -- @@ -26,7 +26,12 @@ -- * Pure random number generator creation RNG, + Seed, mt19937, + mt19937WithState, + rngName, + rngSize, + rngState, -- * Getting and setting the random number generator getRNG, @@ -36,6 +41,9 @@ uniform, uniformInt, normal, + exponential, + levy, + levySkew, poisson, ) where @@ -47,21 +55,21 @@ import Control.Monad.Writer ( MonadWriter(..) ) import Control.Monad.Trans ( MonadTrans(..), MonadIO(..) ) import Data.Word -import System.IO.Unsafe ( unsafePerformIO ) +import System.IO.Unsafe ( unsafePerformIO, unsafeInterleaveIO ) -import GSL.Random.Gen hiding ( mt19937 ) -import qualified GSL.Random.Gen as Gen +import qualified GSL.Random.Gen as GSL import GSL.Random.Dist -- | A Monte Carlo monad with an internal random number generator. -newtype MC a = MC (RNG -> (a,RNG)) +newtype MC a = MC (GSL.RNG -> IO a) -- | Run this Monte Carlo monad with the given initial random number generator, -- getting the result and the new random number generator. runMC :: MC a -> RNG -> (a, RNG) -runMC (MC g) r = - let r' = unsafePerformIO $ cloneRNG r - in r' `seq` g r' +runMC (MC g) (RNG r) = unsafePerformIO $ do + r' <- GSL.cloneRNG r + a <- g r' + return (a,RNG r') {-# NOINLINE runMC #-} -- | Evaluate this Monte Carlo monad and throw away the final random number @@ -75,36 +83,36 @@ execMC g r = snd $ runMC g r unsafeInterleaveMC :: MC a -> MC a -unsafeInterleaveMC (MC m) = MC $ \r -> let - (a,_) = m r - in (a,r) - +unsafeInterleaveMC (MC m) = MC $ \r -> + unsafeInterleaveIO (m r) instance Functor MC where - fmap f (MC m) = MC $ \r -> let - (a,r') = m r - in (f a, r') + fmap f (MC m) = MC $ \r -> + fmap f (m r) instance Monad MC where - return a = MC $ \r -> (a,r) + return a = MC $ \_ -> return a {-# INLINE return #-} (MC m) >>= k = - MC $ \r -> let - (a, r') = m r - (MC m') = k a - in m' r' + MC $ \r -> m r >>= \a -> + let (MC m') = k a + in m' r {-# INLINE (>>=) #-} + + fail s = MC $ \_ -> fail s + {-# INLINE fail #-} -- | A parameterizable Monte Carlo monad for encapsulating an inner -- monad. -newtype MCT m a = MCT (RNG -> m (a,RNG)) +newtype MCT m a = MCT (GSL.RNG -> IO (m a)) -- | Similar to 'runMC'. runMCT :: (Monad m) => MCT m a -> RNG -> m (a,RNG) -runMCT (MCT g) r = - let r' = unsafePerformIO $ cloneRNG r - in r' `seq` g r' +runMCT (MCT g) (RNG r) = unsafePerformIO $ do + r' <- GSL.cloneRNG r + ma <- g r' + return (ma >>= \a -> return (a, RNG r')) {-# NOINLINE runMCT #-} -- | Similar to 'evalMC'. @@ -121,60 +129,69 @@ -- | Take a Monte Carlo computations and lift it to an MCT computation. liftMCT :: (Monad m) => MC a -> MCT m a -liftMCT (MC m) = MCT $ return . m +liftMCT (MC g) = MCT $ \r -> do + a <- g r + return (return a) {-# INLINE liftMCT #-} unsafeInterleaveMCT :: (Monad m) => MCT m a -> MCT m a -unsafeInterleaveMCT (MCT g) = MCT $ \r -> do - ~(a,_) <- g r - return (a,r) +unsafeInterleaveMCT (MCT g) = MCT $ \r -> + unsafeInterleaveIO (g r) {-# INLINE unsafeInterleaveMCT #-} instance (Monad m) => Functor (MCT m) where - fmap f (MCT m) = MCT $ \r -> do - ~(x, r') <- m r - return (f x, r') + fmap f (MCT g) = MCT $ \r -> do + ma <- g r + return (ma >>= return . f) {-# INLINE fmap #-} instance (Monad m) => Monad (MCT m) where - return a = MCT $ \r -> return (a,r) + return a = MCT $ \_ -> return (return a) {-# INLINE return #-} - (MCT m) >>= k = + (MCT g) >>= k = MCT $ \r -> do - ~(a,r') <- m r - let (MCT m') = k a - m' r' - {-# INLINE (>>=) #-} + ma <- g r + return $ ma >>= \a -> + let (MCT m') = k a + in unsafePerformIO $ m' r + {-# NOINLINE (>>=) #-} fail str = MCT $ \_ -> fail str + {-# INLINE fail #-} instance (MonadPlus m) => MonadPlus (MCT m) where mzero = MCT $ \_ -> mzero {-# INLINE mzero #-} (MCT m) `mplus` (MCT n) = - MCT $ \r -> - let r' = unsafePerformIO $ cloneRNG r - in r' `seq` (m r `mplus` n r') - {-# NOINLINE mplus #-} + MCT $ \r -> do + r' <- GSL.cloneRNG r + mr <- m r + nr <- n r' + return (mr `mplus` nr) instance MonadTrans MCT where - lift m = MCT $ \r -> do - a <- m - return (a,r) + lift m = MCT $ \_ -> return m {-# INLINE lift #-} instance (MonadCont m) => MonadCont (MCT m) where callCC f = MCT $ \r -> - callCC $ \c -> - let (MCT m) = (f (\a -> MCT $ \r' -> c (a, r'))) - in m r + return $ callCC $ \k -> + let (MCT m) = f (\a -> MCT $ \_ -> return (k a)) + in unsafePerformIO (m r) + {-# NOINLINE callCC #-} instance (MonadError e m) => MonadError e (MCT m) where - throwError = lift . throwError - (MCT m) `catchError` h = MCT $ \r -> - m r `catchError` \e -> let (MCT m') = h e in m' r + throwError = lift . throwError + {-# INLINE throwError #-} + + (MCT g) `catchError` h = MCT $ \r -> do + ma <- g r + return $ ma `catchError` \e -> + let (MCT m') = h e + in unsafePerformIO (m' r) + {-# NOINLINE catchError #-} instance (MonadIO m) => MonadIO (MCT m) where liftIO = lift . liftIO @@ -182,104 +199,105 @@ instance (MonadReader r m) => MonadReader r (MCT m) where ask = lift ask - local f (MCT m) = MCT $ \r -> - local f (m r) + {-# INLINE ask #-} + + local f (MCT g) = MCT $ \r -> do + ma <- g r + return $ local f ma + {-# INLINE local #-} instance (MonadState s m) => MonadState s (MCT m) where get = lift get + {-# INLINE get #-} + put = lift . put + {-# INLINE put #-} instance (MonadWriter w m) => MonadWriter w (MCT m) where - tell = lift . tell - listen (MCT m) = MCT $ \r -> do - ~((a,r'),w) <- listen (m r) - return ((a,w),r') - pass (MCT m) = MCT $ \r -> pass $ do - ~((a,f),r') <- m r - return ((a,r'),f) + tell = lift . tell + {-# INLINE tell #-} + + listen (MCT g) = MCT $ \r -> do + ma <- g r + return (listen ma) + {-# INLINE listen #-} + + pass (MCT g) = MCT $ \r -> do + maf <- g r + return (pass maf) + {-# INLINE pass #-} ---------------------------- Random Number Generators ----------------------- +-- | The random number generator type associated with 'MC' and 'MCT'. +newtype RNG = RNG GSL.RNG + +-- | The seed type for the random number generators. +type Seed = Word64 + +-- | Get the name of the random number generator algorithm. +rngName :: RNG -> String +rngName (RNG r) = unsafePerformIO $ GSL.getName r +{-# NOINLINE rngName #-} + +-- | Get the size of the generator state, in bytes. +rngSize :: RNG -> Int +rngSize (RNG r) = fromIntegral $ unsafePerformIO $ GSL.getSize r +{-# NOINLINE rngSize #-} + +-- | Get the state of the generator. +rngState :: RNG -> [Word8] +rngState (RNG r) = unsafePerformIO $ GSL.getState r +{-# NOINLINE rngState #-} + getRNG :: MC RNG -getRNG = MC $ getHelp +getRNG = MC (\r -> liftM RNG $ GSL.cloneRNG r) {-# INLINE getRNG #-} -getHelp :: RNG -> (RNG,RNG) -getHelp r = unsafePerformIO $ do - r' <- cloneRNG r - r' `seq` return (r',r) -{-# NOINLINE getHelp #-} - setRNG :: RNG -> MC () -setRNG r' = MC $ setHelp r' +setRNG (RNG r') = MC $ \r -> GSL.copyRNG r r' {-# INLINE setRNG #-} -setHelp :: RNG -> RNG -> ((),RNG) -setHelp r' r = unsafePerformIO $ do - io <- copyRNG r r' - io `seq` return ((),r) -{-# NOINLINE setHelp #-} - -- | Get a Mersenne Twister random number generator seeded with the given -- value. -mt19937 :: Word64 -> RNG +mt19937 :: Seed -> RNG mt19937 s = unsafePerformIO $ do - r <- newRNG Gen.mt19937 - setSeed r s - return r + r <- GSL.newRNG GSL.mt19937 + GSL.setSeed r s + return (RNG r) {-# NOINLINE mt19937 #-} +-- | Get a Mersenne Twister seeded with the given state. +mt19937WithState :: [Word8] -> RNG +mt19937WithState xs = unsafePerformIO $ do + r <- GSL.newRNG GSL.mt19937 + GSL.setState r xs + return (RNG r) +{-# NOINLINE mt19937WithState #-} -------------------------- Random Number Distributions ---------------------- uniform :: Double -> Double -> MC Double -uniform a b = MC $ uniformHelp a b -{-# INLINE uniform #-} - -uniformHelp :: Double -> Double -> RNG -> (Double,RNG) -uniformHelp 0 1 r = unsafePerformIO $ do - x <- getUniform r - x `seq` return (x,r) -uniformHelp a b r = unsafePerformIO $ do - x <- getFlat r a b - x `seq` return (x,r) -{-# NOINLINE uniformHelp #-} +uniform 0 1 = MC $ \r -> GSL.getUniform r +uniform a b = MC $ \r -> getFlat r a b uniformInt :: Int -> MC Int -uniformInt n = MC $ uniformIntHelp n -{-# INLINE uniformInt #-} - -uniformIntHelp :: Int -> RNG -> (Int,RNG) -uniformIntHelp n r = unsafePerformIO $ do - x <- getUniformInt r n - x `seq` return (x,r) -{-# NOINLINE uniformIntHelp #-} +uniformInt n = MC $ \r -> GSL.getUniformInt r n normal :: Double -> Double -> MC Double -normal mu sigma = MC $ normalHelp mu sigma -{-# INLINE normal #-} +normal 0 1 = MC $ \r -> getUGaussianRatioMethod r +normal mu 1 = MC $ \r -> liftM (mu +) (getUGaussianRatioMethod r) +normal 0 sigma = MC $ \r -> getGaussianRatioMethod r sigma +normal mu sigma = MC $ \r -> liftM (mu +) (getGaussianRatioMethod r sigma) -normalHelp :: Double -> Double -> RNG -> (Double,RNG) -normalHelp 0 1 r = unsafePerformIO $ do - x <- getUGaussianRatioMethod r - x `seq` return (x,r) -normalHelp mu 1 r = unsafePerformIO $ do - x <- liftM (mu +) $ getUGaussianRatioMethod r - x `seq` return (x,r) -normalHelp 0 sigma r = unsafePerformIO $ do - x <- getGaussianRatioMethod r sigma - x `seq` return (x,r) -normalHelp mu sigma r = unsafePerformIO $ do - x <- liftM (mu +) $ getGaussianRatioMethod r sigma - x `seq` return (x,r) -{-# NOINLINE normalHelp #-} +exponential :: Double -> MC Double +exponential mu = MC $ \r -> getExponential r mu poisson :: Double -> MC Int -poisson mu = MC $ poissonHelp mu -{-# INLINE poisson #-} +poisson mu = MC $ \r -> getPoisson r mu -poissonHelp :: Double -> RNG -> (Int,RNG) -poissonHelp mu r = unsafePerformIO $ do - x <- getPoisson r mu - x `seq` return (x,r) -{-# NOINLINE poissonHelp #-} +levy :: Double -> Double -> MC Double +levy c alpha = MC $ \r -> getLevy r c alpha + +levySkew :: Double -> Double -> Double -> MC Double +levySkew c alpha beta = MC $ \r -> getLevySkew r c alpha beta diff --git a/lib/Control/Monad/MC/Repeat.hs b/lib/Control/Monad/MC/Repeat.hs --- a/lib/Control/Monad/MC/Repeat.hs +++ b/lib/Control/Monad/MC/Repeat.hs @@ -1,49 +1,31 @@ ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Repeat --- Copyright : Copyright (c) , Patrick Perry +-- Copyright : Copyright (c) 2010, Patrick Perry -- License : BSD3 --- Maintainer : Patrick Perry +-- Maintainer : Patrick Perry -- Stability : experimental -- module Control.Monad.MC.Repeat ( - -- * Averaging functions + -- * Repeating computations repeatMC, - repeatMCWith, - - module Control.Monad.MC.Summary, + replicateMC, ) where -import Control.Monad import Control.Monad.MC.Base -import Control.Monad.MC.Summary -import Data.List( foldl' ) --- | Repeat a Monte Carlo generator the given number of times and return --- the sample summary statistics. Note that this only works with --- @Double@s. -repeatMC :: (MonadMC m) - => Int - -> m Double - -> m Summary -repeatMC = repeatMCWith update summary +-- | Produce a lazy infinite list of values from the given Monte Carlo +-- generator. +repeatMC :: (MonadMC m) => m a -> m [a] +repeatMC = interleaveSequence . repeat {-# INLINE repeatMC #-} - --- | Generalized version of 'repeatMC'. Run a Monte Carlo generator --- the given number of times and accumulate the results. The accumulator --- is strictly evaluated. -repeatMCWith :: (MonadMC m) - => (a -> b -> a) -- ^ accumulator - -> a -- ^ initial value - -> Int -- ^ number of repetitions - -> m b -- ^ generator - -> m a -repeatMCWith f a n mb = do - bs <- interleaveSequence $ replicate n mb - return $! foldl' f a bs -{-# INLINE repeatMCWith #-} - + +-- | Produce a lazy list of the given length using the specified +-- generator. +replicateMC :: (MonadMC m) => Int -> m a -> m [a] +replicateMC n = interleaveSequence . replicate n +{-# INLINE replicateMC #-} interleaveSequence :: (MonadMC m) => [m a] -> m [a] interleaveSequence [] = return [] diff --git a/lib/Control/Monad/MC/Sample.hs b/lib/Control/Monad/MC/Sample.hs --- a/lib/Control/Monad/MC/Sample.hs +++ b/lib/Control/Monad/MC/Sample.hs @@ -1,10 +1,9 @@ -{-# LANGUAGE ScopedTypeVariables #-} ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Sample --- Copyright : Copyright (c) , Patrick Perry +-- Copyright : Copyright (c) 2010, Patrick Perry -- License : BSD3 --- Maintainer : Patrick Perry +-- Maintainer : Patrick Perry -- Stability : experimental -- @@ -13,15 +12,18 @@ sample, sampleWithWeights, sampleSubset, + sampleSubset', -- * Sampling @Int@s sampleInt, sampleIntWithWeights, sampleIntSubset, + sampleIntSubset', -- * Shuffling shuffle, shuffleInt, + shuffleInt', ) where import Control.Monad @@ -29,65 +31,78 @@ import Control.Monad.MC.Base import Control.Monad.MC.Walker -import Data.Array.Base -import Data.Array.IArray -import Data.Array.ST -import Data.Array.Vector +import Data.Vector.Unboxed( MVector, Unbox ) +import qualified Data.Vector as BV +import qualified Data.Vector.Mutable as BMV +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector.Generic.Mutable as MV --- | @sample n xs@ samples a value uniformly from @take n xs@. The results --- are undefined if @length xs@ is less than @n@. -sample :: (MonadMC m) => Int -> [a] -> m a -sample n xs = - sampleHelp n xs $ sampleInt n +-- | @sample xs@ samples a value uniformly from the elements of @xs@. The +-- results are undefined if @length xs@ is zero. +sample :: (MonadMC m) => [a] -> m a +sample xs = let + n = length xs + in sampleHelp n xs $ sampleInt n {-# INLINE sample #-} --- | @sampleWithWeights ws n xs@ samples a value from @take n xs@, putting --- weight @ws !! i@ on element @xs !! i@. The results --- are undefined if @length xs@ or @length ws@ is less than @n@. -sampleWithWeights :: (MonadMC m) => [Double] -> Int -> [a] -> m a -sampleWithWeights ws n xs = - sampleHelp n xs $ sampleIntWithWeights ws n +-- | @sampleWithWeights wxs@ samples a value from the list with the given +-- weight. +sampleWithWeights :: (MonadMC m) => [(Double, a)] -> m a +sampleWithWeights wxs = let + (ws,xs) = unzip wxs + n = length xs + in sampleHelp n xs $ sampleIntWithWeights ws n {-# INLINE sampleWithWeights #-} --- | @sampleSubset k n xs@ samples a subset of size @k@ from @take n xs@ by +-- | @sampleSubset k xs@ samples a subset of size @k@ from @xs@ by -- sampling without replacement. The return value is a list of length @k@ -- with the elements in the subset in the order that they were sampled. Note --- also that the elements are lazily generated. The results are undefined --- if @k > n@ or if @length xs < n@. -sampleSubset :: (MonadMC m) => Int -> Int -> [a] -> m [a] -sampleSubset k n xs = - sampleListHelp n xs $ sampleIntSubset k n +-- also that the elements are lazily generated. +sampleSubset :: (MonadMC m) => Int -> [a] -> m [a] +sampleSubset k xs = let + n = length xs + in sampleListHelp n xs $ sampleIntSubset k n {-# INLINE sampleSubset #-} +-- | Strict version of 'sampleSubset'. +sampleSubset' :: (MonadMC m) => Int -> [a] -> m [a] +sampleSubset' k xs = do + s <- sampleSubset k xs + length s `seq` return s +{-# INLINE sampleSubset' #-} + sampleHelp :: (Monad m) => Int -> [a] -> m Int -> m a -sampleHelp n (xs :: [a]) f = let - arr = listArray (0,n-1) xs :: Array Int a - in liftM (unsafeAt arr) f +sampleHelp _n xs f = let + arr = BV.fromList xs + in liftM (BV.unsafeIndex arr) f +{-# INLINE sampleHelp #-} -sampleHelpUA :: (UA a, Monad m) => Int -> [a] -> m Int -> m a -sampleHelpUA n xs f = let - arr = newU n (\marr -> zipWithM_ (writeMU marr) [0..n-1] xs) - in liftM (indexU arr) f +sampleHelpU :: (Unbox a, Monad m) => Int -> [a] -> m Int -> m a +sampleHelpU _n xs f = let + arr = V.fromList xs + in liftM (V.unsafeIndex arr) f +{-# INLINE sampleHelpU #-} {-# RULES "sampleHelp/Double" forall n xs f. - sampleHelp n (xs :: [Double]) f = sampleHelpUA n xs f #-} + sampleHelp n (xs :: [Double]) f = sampleHelpU n xs f #-} {-# RULES "sampleHelp/Int" forall n xs f. - sampleHelp n (xs :: [Int]) f = sampleHelpUA n xs f #-} + sampleHelp n (xs :: [Int]) f = sampleHelpU n xs f #-} sampleListHelp :: (Monad m) => Int -> [a] -> m [Int] -> m [a] -sampleListHelp n (xs :: [a]) f = let - arr = listArray (0,n-1) xs :: Array Int a - in liftM (map $ unsafeAt arr) f +sampleListHelp _n xs f = let + arr = BV.fromList xs + in liftM (map $ BV.unsafeIndex arr) f +{-# INLINE sampleListHelp #-} -sampleListHelpUA :: (UA a, Monad m) => Int -> [a] -> m [Int] -> m [a] -sampleListHelpUA n xs f = let - arr = newU n (\marr -> zipWithM_ (writeMU marr) [0..n-1] xs) - in liftM (map $ indexU arr) f +sampleListHelpU :: (Unbox a, Monad m) => Int -> [a] -> m [Int] -> m [a] +sampleListHelpU _n xs f = let + arr = V.fromList xs + in liftM (map $ V.unsafeIndex arr) f {-# RULES "sampleListHelp/Double" forall n xs f. - sampleListHelp n (xs :: [Double]) f = sampleListHelpUA n xs f #-} + sampleListHelp n (xs :: [Double]) f = sampleListHelpU n xs f #-} {-# RULES "sampleListHelp/Int" forall n xs f. - sampleListHelp n (xs :: [Int]) f = sampleListHelpUA n xs f #-} + sampleListHelp n (xs :: [Int]) f = sampleListHelpU n xs f #-} -- | @sampleInt n@ samples integers uniformly from @[ 0..n-1 ]@. It is an -- error to call this function with a non-positive @n@. @@ -116,8 +131,8 @@ | otherwise = do us <- randomIndices k n return $ runST $ do - ints <- newMU n - sequence_ [ writeMU ints i i | i <- [0 .. n-1] ] + ints <- MV.new n :: ST s (MVector s Int) + sequence_ [ MV.unsafeWrite ints i i | i <- [0 .. n-1] ] sampleIntSubsetHelp ints us (n-1) where randomIndices k' n' | k' == 0 = return [] @@ -128,51 +143,60 @@ sampleIntSubsetHelp _ [] _ = return [] sampleIntSubsetHelp ints (u:us) n' = unsafeInterleaveST $ do - i <- readMU ints u - writeMU ints u =<< readMU ints n' + i <- MV.unsafeRead ints u + MV.unsafeWrite ints u =<< MV.unsafeRead ints n' is <- sampleIntSubsetHelp ints us (n'-1) return (i:is) {-# INLINE sampleIntSubset #-} --- | @shuffle n xs@ randomly permutes the list @take n xs@ and returns +-- | Strict version of 'sampleIntSubset'. +sampleIntSubset' :: (MonadMC m) => Int -> Int -> m [Int] +sampleIntSubset' k n = do + s <- sampleIntSubset k n + length s `seq` return s +{-# INLINE sampleIntSubset' #-} + +-- | @shuffle xs@ randomly permutes the list @xs@ and returns -- the result. All permutations of the elements of @xs@ are equally --- likely. The results are undefined if @length xs@ is less than @n@. -shuffle :: (MonadMC m) => Int -> [a] -> m [a] -shuffle n (xs :: [a]) = - shuffleInt n >>= \swaps -> (return . runST) $ do - marr <- newListArray (0,n-1) xs :: ST s (STArray s Int a) - mapM_ (swap marr) swaps - getElems marr +-- likely. +shuffle :: (MonadMC m) => [a] -> m [a] +shuffle xs = let + n = length xs + in shuffleInt n >>= \swaps -> (return . BV.toList . BV.create) $ do + marr <- MV.new n :: ST s (BMV.MVector s a) + zipWithM_ (MV.unsafeWrite marr) [0 .. n-1] xs + mapM_ (swap marr) swaps + return marr where swap marr (i,j) | i == j = return () | otherwise = do - x <- unsafeRead marr i - y <- unsafeRead marr j - unsafeWrite marr i y - unsafeWrite marr j x + x <- MV.unsafeRead marr i + y <- MV.unsafeRead marr j + MV.unsafeWrite marr i y + MV.unsafeWrite marr j x {-# INLINE shuffle #-} -shuffleUA :: (UA a, MonadMC m) => Int -> [a] -> m [a] -shuffleUA n (xs :: [a]) = - shuffleInt n >>= \swaps -> (return . runST) $ do - marr <- newMU n - zipWithM_ (writeMU marr) [0 .. n-1] xs - mapM_ (swap marr) swaps - arr <- unsafeFreezeAllMU marr - return $ fromU arr +shuffleU :: (Unbox a, MonadMC m) => [a] -> m [a] +shuffleU xs = let + n = length xs + in shuffleInt n >>= \swaps -> (return . V.toList . V.create) $ do + marr <- MV.new n + zipWithM_ (MV.unsafeWrite marr) [0 .. n-1] xs + mapM_ (swap marr) swaps + return marr where swap marr (i,j) | i == j = return () | otherwise = do - x <- readMU marr i - y <- readMU marr j - writeMU marr i y - writeMU marr j x -{-# INLINE shuffleUA #-} + x <- MV.unsafeRead marr i + y <- MV.unsafeRead marr j + MV.unsafeWrite marr i y + MV.unsafeWrite marr j x +{-# INLINE shuffleU #-} -{-# RULES "shuffle/Double" forall n xs. - shuffle n (xs :: [Double]) = shuffleUA n xs #-} -{-# RULES "shuffle/Int" forall n xs. - shuffle n (xs :: [Int]) = shuffleUA n xs #-} +{-# RULES "shuffle/Double" forall xs. + shuffle (xs :: [Double]) = shuffleU xs #-} +{-# RULES "shuffle/Int" forall xs. + shuffle (xs :: [Int]) = shuffleU xs #-} -- | @shuffleInt n@ generates a sequence of swaps equivalent to a @@ -187,3 +211,10 @@ return $ (i-1,j):ijs in shuffleIntHelp n {-# INLINE shuffleInt #-} + +-- | Strict version of 'shuffleInt'. +shuffleInt' :: (MonadMC m) => Int -> m [(Int,Int)] +shuffleInt' n = do + ss <- shuffleInt n + length ss `seq` return ss +{-# INLINE shuffleInt' #-} diff --git a/lib/Control/Monad/MC/Summary.hs b/lib/Control/Monad/MC/Summary.hs deleted file mode 100644 --- a/lib/Control/Monad/MC/Summary.hs +++ /dev/null @@ -1,96 +0,0 @@ ------------------------------------------------------------------------------ --- | --- Module : Control.Monad.MC.Summary --- Copyright : Copyright (c) , Patrick Perry --- License : BSD3 --- Maintainer : Patrick Perry --- Stability : experimental --- - -module Control.Monad.MC.Summary ( - -- * Summary statistics - -- ** The @Summary@ data type - Summary, - summary, - update, - - -- ** @Summary@ properties - sampleSize, - sampleMean, - sampleVar, - sampleSD, - sampleSE, - sampleCI, - sampleMin, - sampleMax, - - ) where - -import GSL.Random.Dist( ugaussianPInv ) - --- | A type for storing summary statistics for a data set including --- sample size, min and max values, and first and second moments. -data Summary = S {-# UNPACK #-} !Int -- sample size - {-# UNPACK #-} !Double -- sample mean - {-# UNPACK #-} !Double -- sum of squares - {-# UNPACK #-} !Double -- sample min - {-# UNPACK #-} !Double -- sample max - --- | Get an empty summary. -summary :: Summary -summary = S 0 0 0 (1/0) (-1/0) - --- | Update the summary with a data point. --- Running mean and variance computed as in Knuth, Vol 2, page 232, --- 3rd edition, see http://www.johndcook.com/standard_deviation.html for --- a description. -update :: Summary -> Double -> Summary -update (S n m s l h) x = - let n' = n+1 - delta = x - m - m' = m + delta / fromIntegral n' - s' = s + delta*(x - m') - l' = if x < l then x else l - h' = if x > h then x else h - in S n' m' s' l' h' - --- | Get the sample size. -sampleSize :: Summary -> Int -sampleSize (S n _ _ _ _) = n - --- | Get the sample mean. -sampleMean :: Summary -> Double -sampleMean (S _ m _ _ _) = m - --- | Get the sample variance. -sampleVar :: Summary -> Double -sampleVar (S n _ s _ _) = s / fromIntegral (n - 1) - --- | Get the sample standard deviation. -sampleSD :: Summary -> Double -sampleSD s = sqrt (sampleVar s) - --- | Get the sample standard error. -sampleSE :: Summary -> Double -sampleSE s = sqrt (sampleVar s / fromIntegral (sampleSize s)) - --- | Get a Central Limit Theorem-based confidence interval for the mean --- with the specified coverage level. The level must be in the range @(0,1)@. -sampleCI :: Double -> Summary -> (Double,Double) -sampleCI level s | not (level > 0 && level < 1) = - error "level must be between 0 and 1" - | otherwise = - let alpha = (0.5 - level) + 0.5 - z = -(ugaussianPInv (0.5*alpha)) - se = sampleSE s - delta = z*se - xbar = sampleMean s - in (xbar-delta, xbar+delta) - --- | Get the minimum of the sample. -sampleMin :: Summary -> Double -sampleMin (S _ _ _ l _) = l - --- | Get the maximum of the sample. -sampleMax :: Summary -> Double -sampleMax (S _ _ _ _ h) = h diff --git a/lib/Control/Monad/MC/Walker.hs b/lib/Control/Monad/MC/Walker.hs --- a/lib/Control/Monad/MC/Walker.hs +++ b/lib/Control/Monad/MC/Walker.hs @@ -1,10 +1,9 @@ -{-# LANGUAGE TypeOperators #-} ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Walker --- Copyright : Copyright (c) , Patrick Perry +-- Copyright : Copyright (c) 2010, Patrick Perry -- License : BSD3 --- Maintainer : Patrick Perry +-- Maintainer : Patrick Perry -- Stability : experimental -- -- An implementation of Walker's Alias method for sampling from discrete @@ -21,53 +20,55 @@ import Control.Monad import Control.Monad.ST -import Data.Array.Vector +import Data.Vector.Unboxed( Vector, MVector ) +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector.Generic.Mutable as MV -- | The table, which represents an equiprobable mixture of two-point -- distributions. The @l@th entry of the table represents a mixture -- distribution with weight @q[l]@ on @l@ and weight @(1-q[l])@ on @j[l]@. -- The @l@th element of the table stores the pair @q[l] :*: j[l]@. -newtype Table = T (UArr (Double :*: Int)) +newtype Table = T (Vector (Double, Int)) -- | Get the @i@th mixture component. That is, return @q[i]@ and @j[i]@, -- where the @i@th mixture component puts mass @q[i]@ on @i@ and mass -- @1 - q[i]@ on @j[i]@. component :: Table -> Int -> (Double,Int) component (T qjs) i = let - (q' :*: j) = indexU qjs i + (q', j) = V.unsafeIndex qjs i q = q' - fromIntegral i in (q,j) -- | Compute the table for use in Walker's aliasing method. computeTable :: Int -> [Double] -> Table -computeTable n ws = runST $ do +computeTable n ws = T $ V.create $ do (qjs, sets) <- initTable n ws breakLarger qjs sets scaleTable qjs - liftM T $ unsafeFreezeAllMU qjs + return qjs -- | Given an alias table and a number in the range [0,1), -- get the corresponding sample in the table. indexTable :: Table -> Double -> Int indexTable (T qjs) u = let - n = lengthU qjs + n = V.length qjs nu = u * fromIntegral n l = floor nu - (ql :*: jl) = indexU qjs l + (ql,jl) = V.unsafeIndex qjs l in if nu < ql then l else jl -- | Get the size of the table tableSize :: Table -> Int -tableSize (T qjs) = lengthU qjs +tableSize (T qjs) = V.length qjs -- | An intermediate result for use in computing a Table. -type STTable s = MUArr (Double :*: Int) s +type STTable s = MVector s (Double, Int) -- | A partition of indices into the sets /Greater/ and /Smaller/. The -- indices of the /Smaller/ set are stored in positions @0, ..., numSmall - 1@, -- and the indices of the /Greater/ set are stored in positions -- @numSmall, ..., n-1@, where @n@ is the size of the underlying array. -data STPartition s = P !(MUArr Int s) +data STPartition s = P !(MVector s Int) !Int -- | Given a list of weights, @ws@, compute corresponding probabilities, @ps@, @@ -77,15 +78,15 @@ initTable :: Int -> [Double] -> ST s (STTable s, STPartition s) initTable n ws = do when (n < 0) $ fail "negative table size" - sets <- newMU n :: ST s (MUArr Int s) - qjs <- newMU n :: ST s (MUArr (Double :*: Int) s) + sets <- MV.new n :: ST s (MVector s Int) + qjs <- MV.new n :: ST s (MVector s (Double, Int)) -- Store the weights in the table and compute their total. total <- foldM (\current (i,w) -> do if w >= 0 then do - writeMU qjs i (w :*: i) + MV.unsafeWrite qjs i (w,i) return $! current + w else fail $ "negative probability" ) @@ -99,15 +100,15 @@ let scale = fromIntegral n / total nsmall <- liftM fst $ foldM (\(smaller,greater) i -> do - p <- liftM fstS $ readMU qjs i + p <- liftM fst $ MV.unsafeRead qjs i let q = scale*p - writeMU qjs i (q :*: i) + MV.unsafeWrite qjs i (q,i) if q < 1 then do - writeMU sets smaller i + MV.unsafeWrite sets smaller i return (smaller+1,greater) else do - writeMU sets greater i + MV.unsafeWrite sets greater i return (smaller,greater-1) ) (0,n-1) [0 .. n-1] @@ -121,24 +122,24 @@ breakLarger :: STTable s -> STPartition s -> ST s () breakLarger qjs (P sets nsmall) | nsmall == 0 = return () | otherwise = let - n = lengthMU qjs + n = MV.length qjs breakLargerHelp nsmall' i | nsmall' == n = return () | i == n = return () | otherwise = do -- while Greater is not empty -- choose k from Greater, l from Smaller - k <- readMU sets $ nsmall' - l <- readMU sets $ i - qk <- liftM fstS $ readMU qjs k - ql <- liftM fstS $ readMU qjs l + k <- MV.unsafeRead sets $ nsmall' + l <- MV.unsafeRead sets $ i + qk <- liftM fst $ MV.unsafeRead qjs k + ql <- liftM fst $ MV.unsafeRead qjs l -- set jl := k, finalize (ql,jl) let jl = k - writeMU qjs l (ql :*: jl) + MV.unsafeWrite qjs l (ql,jl) -- set qk := qk - (1-ql) let qk' = qk - (1-ql) - writeMU qjs k (qk' :*: k) + MV.unsafeWrite qjs k (qk',k) -- if qk' < 1, move k from Greater to Smaller let nsmall'' = if qk' < 1 then nsmall'+1 else nsmall' @@ -152,8 +153,8 @@ -- from the table. scaleTable :: STTable s -> ST s () scaleTable qjs = let - n = lengthMU qjs in + n = MV.length qjs in forM_ [ 0..(n-1) ] $ \l -> do - (ql :*: jl) <- readMU qjs l - writeMU qjs l ((ql + fromIntegral l) :*: jl) + (ql, jl) <- MV.unsafeRead qjs l + MV.unsafeWrite qjs l ((ql + fromIntegral l), jl) diff --git a/lib/Data/Summary.hs b/lib/Data/Summary.hs new file mode 100644 --- /dev/null +++ b/lib/Data/Summary.hs @@ -0,0 +1,15 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Data.Summary +-- Copyright : Copyright (c) 2010, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- +-- Summary Statistics +-- +module Data.Summary ( + module Data.Summary.Double + ) where + +import Data.Summary.Double diff --git a/lib/Data/Summary/Bool.hs b/lib/Data/Summary/Bool.hs new file mode 100644 --- /dev/null +++ b/lib/Data/Summary/Bool.hs @@ -0,0 +1,86 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Data.Summary.Bool +-- Copyright : Copyright (c) 2010, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- +-- Summary statistics for @Bool@s. +-- + +module Data.Summary.Bool ( + -- * The @Summary@ data type + Summary, + summary, + update, + + -- * @Summary@ properties + sampleSize, + count, + sampleMean, + sampleSE, + sampleCI, + + ) where + +import Data.List( foldl' ) +import Text.Printf + +import Data.Summary.Utils + + +-- | A type for storing summary statistics for a data set of +-- booleans. Specifically, this just keeps track of the number +-- of 'True' events and gives estimates for the success +-- probability. 'True' is interpreted as a one, and 'False' +-- is interpreted as a zero. +data Summary = S {-# UNPACK #-} !Int -- sample size + {-# UNPACK #-} !Int -- number of successes + +instance Show Summary where + show s@(S n c) = + printf " sample size: %d" n + ++ printf "\n successes: %g" c + ++ printf "\n proportion: %g" (sampleMean s) + ++ printf "\n SE: %g" (sampleSE s) + ++ printf "\n 99%% CI: (%g, %g)" c1 c2 + where (c1,c2) = sampleCI 0.99 s + +-- | Get a summary of a list of values. +summary :: [Bool] -> Summary +summary = foldl' update empty + +-- | Get an empty summary. +empty :: Summary +empty = S 0 0 + +-- | Update the summary with a data point. +update :: Summary -> Bool -> Summary +update (S n c) i = + let n' = n+1 + c' = if i then c+1 else c + in S n' c' + +-- | Get the sample size. +sampleSize :: Summary -> Int +sampleSize (S n _) = n + +-- | Get the number of 'True' values +count :: Summary -> Int +count (S _ c) = c + +-- | Get the proportion of 'True' events. +sampleMean :: Summary -> Double +sampleMean (S n c) = fromIntegral c / fromIntegral n + +-- | Get the standard error for the sample proportion. +sampleSE :: Summary -> Double +sampleSE s = sqrt (p*(1-p) / n) + where p = sampleMean s + n = fromIntegral $ sampleSize s + +-- | Get a Central Limit Theorem-based confidence interval for the mean +-- with the specified coverage level. The level must be in the range @(0,1)@. +sampleCI :: Double -> Summary -> (Double,Double) +sampleCI level s = interval level (sampleMean s) (sampleSE s) diff --git a/lib/Data/Summary/Double.hs b/lib/Data/Summary/Double.hs new file mode 100644 --- /dev/null +++ b/lib/Data/Summary/Double.hs @@ -0,0 +1,107 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Data.Summary.Double +-- Copyright : Copyright (c) 2010, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- +-- Summary statistics for @Double@s. +-- + +module Data.Summary.Double ( + -- * The @Summary@ data type + Summary, + summary, + update, + + -- * @Summary@ properties + sampleSize, + sampleMin, + sampleMax, + sampleMean, + sampleSE, + sampleVar, + sampleSD, + sampleCI, + + ) where + +import Data.List( foldl' ) +import Text.Printf + +import Data.Summary.Utils + + +-- | A type for storing summary statistics for a data set including +-- sample size, min and max values, and first and second moments. +data Summary = S {-# UNPACK #-} !Int -- sample size + {-# UNPACK #-} !Double -- sample mean + {-# UNPACK #-} !Double -- sum of squares + {-# UNPACK #-} !Double -- sample min + {-# UNPACK #-} !Double -- sample max + +instance Show Summary where + show s@(S n mu _ l h) = + printf " sample size: %d" n + ++ printf "\n min: %g" l + ++ printf "\n max: %g" h + ++ printf "\n mean: %g" mu + ++ printf "\n SE: %g" (sampleSE s) + ++ printf "\n 99%% CI: (%g, %g)" c1 c2 + where (c1,c2) = sampleCI 0.99 s + +-- | Get a summary of a list of values. +summary :: [Double] -> Summary +summary = foldl' update empty + +-- | Get an empty summary. +empty :: Summary +empty = S 0 0 0 (1/0) (-1/0) + +-- | Update the summary with a data point. +-- Running mean and variance computed as in Knuth, Vol 2, page 232, +-- 3rd edition, see http://www.johndcook.com/standard_deviation.html for +-- a description. +update :: Summary -> Double -> Summary +update (S n m s l h) x = + let n' = n+1 + delta = x - m + m' = m + delta / fromIntegral n' + s' = s + delta*(x - m') + l' = if x < l then x else l + h' = if x > h then x else h + in S n' m' s' l' h' + +-- | Get the sample size. +sampleSize :: Summary -> Int +sampleSize (S n _ _ _ _) = n + +-- | Get the sample mean. +sampleMean :: Summary -> Double +sampleMean (S _ m _ _ _) = m + +-- | Get the sample variance. +sampleVar :: Summary -> Double +sampleVar (S n _ s _ _) = s / fromIntegral (n - 1) + +-- | Get the sample standard deviation. +sampleSD :: Summary -> Double +sampleSD s = sqrt (sampleVar s) + +-- | Get the sample standard error. +sampleSE :: Summary -> Double +sampleSE s = sqrt (sampleVar s / fromIntegral (sampleSize s)) + +-- | Get a Central Limit Theorem-based confidence interval for the mean +-- with the specified coverage level. The level must be in the range @(0,1)@. +sampleCI :: Double -> Summary -> (Double,Double) +sampleCI level s = interval level (sampleMean s) (sampleSE s) + +-- | Get the minimum of the sample. +sampleMin :: Summary -> Double +sampleMin (S _ _ _ l _) = l + +-- | Get the maximum of the sample. +sampleMax :: Summary -> Double +sampleMax (S _ _ _ _ h) = h diff --git a/lib/Data/Summary/Utils.hs b/lib/Data/Summary/Utils.hs new file mode 100644 --- /dev/null +++ b/lib/Data/Summary/Utils.hs @@ -0,0 +1,36 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Data.Summary.Utils +-- Copyright : Copyright (c) 2010, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- +-- Utilities for data summaries. +-- + +module Data.Summary.Utils ( + interval, + inInterval, + ) where + +import GSL.Random.Dist( ugaussianPInv ) + +-- | Get a Central Limit Theorem-based confidence interval for the +-- population mean with the specified coverage level. The level must +-- be in the range @(0,1)@. +interval :: Double -- ^ the confidence level + -> Double -- ^ the sample mean + -> Double -- ^ the sample standard error + -> (Double,Double) +interval level xbar se | not (level > 0 && level < 1) = + error "level must be between 0 and 1" + | otherwise = + let alpha = (0.5 - level) + 0.5 + z = -(ugaussianPInv (0.5*alpha)) + delta = z*se + in (xbar-delta, xbar+delta) + +-- | Tests if the value is in the open interval (a,b) +inInterval :: Double -> (Double,Double) -> Bool +x `inInterval` (a,b) = x > a && x < b diff --git a/monte-carlo.cabal b/monte-carlo.cabal --- a/monte-carlo.cabal +++ b/monte-carlo.cabal @@ -1,10 +1,10 @@ name: monte-carlo -version: 0.2 +version: 0.3 license: BSD3 license-file: LICENSE author: Patrick Perry -maintainer: Patrick Perry -homepage: http://github.com/patperry/monte-carlo +maintainer: Patrick Perry +homepage: http://github.com/patperry/hs-monte-carlo category: Math synopsis: A monad and transformer for Monte Carlo calculations. description: A monad and transformer for Monte Carlo calculations. The @@ -16,33 +16,37 @@ build-type: Simple stability: experimental cabal-version: >= 1.2.3 -extra-source-files: examples/Pi.hs, examples/Sampling.hs examples/Poker.hs - tests/Main.hs tests/Makefile +extra-source-files: NEWS examples/Binomial.hs examples/Pi.lhs + examples/Poker.hs examples/Queue.hs tests/Main.hs + tests/Makefile library - build-depends: array, base, mtl, gsl-random >=0.2.3, uvector - exposed-modules: Control.Monad.MC Control.Monad.MC.Class + Control.Monad.MC.GSL + Data.Summary + Data.Summary.Bool + Data.Summary.Double + Data.Summary.Utils other-modules: Control.Monad.MC.Base - Control.Monad.MC.GSL Control.Monad.MC.GSLBase Control.Monad.MC.Repeat Control.Monad.MC.Sample - Control.Monad.MC.Summary Control.Monad.MC.Walker extensions: - FlexibleContexts, - FlexibleInstances, + FlexibleInstances, MultiParamTypeClasses, - ScopedTypeVariables, TypeFamilies, - TypeOperators, UndecidableInstances + + build-depends: base >= 4 && < 5, + gsl-random >= 0.3.1, + mtl >= 1.1 && < 1.2, + vector >= 0.6 && < 0.7 hs-source-dirs: lib ghc-options: -Wall diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -10,6 +10,8 @@ import System.Random import Text.Printf import Test.QuickCheck +import Test.Framework +import Test.Framework.Providers.QuickCheck2 import Control.Monad.MC.Walker @@ -25,9 +27,10 @@ i = indexTable table u in i >= 0 && i < n && (ws !! i > 0) -tests_Walker = [ ("table probabilities", mytest prop_table_probs) - , ("table indexing" , mytest prop_table_index) - ] +tests_Walker = testGroup "Walker" + [ testProperty "table probabilities" prop_table_probs + , testProperty "table indexing" prop_table_index + ] probOf table i = (((sum . map ((1-) . fst) . filter ((==i) . snd)) @@ -70,103 +73,14 @@ data Weights = Weights Int [Double] deriving Show instance Arbitrary Weights where arbitrary = do - n <- posInt + n <- choose (1, 500) ws <- weights n return $ Weights n ws - coarbitrary (Weights n ws) = - coarbitrary (n,ws) - data Unif = Unif Double deriving Show instance Arbitrary Unif where arbitrary = liftM Unif unif - coarbitrary (Unif u) = coarbitrary u ------------------------------------------------------------------------- --- --- QC driver ( taken from xmonad-0.6 ) --- -debug = False - -mytest :: Testable a => a -> Int -> IO (Bool, Int) -mytest a n = mycheck defaultConfig - { configMaxTest=n - , configEvery = \n args -> let s = show n in s ++ [ '\b' | _ <- s ] } a - -- , configEvery= \n args -> if debug then show n ++ ":\n" ++ unlines args else [] } a - -mycheck :: Testable a => Config -> a -> IO (Bool, Int) -mycheck config a = do - rnd <- newStdGen - mytests config (evaluate a) rnd 0 0 [] - -mytests :: Config -> Gen Result -> StdGen -> Int -> Int -> [[String]] -> IO (Bool, Int) -mytests config gen rnd0 ntest nfail stamps - | ntest == configMaxTest config = done "OK," ntest stamps >> return (True, ntest) - | nfail == configMaxFail config = done "Arguments exhausted after" ntest stamps >> return (True, ntest) - | otherwise = - do putStr (configEvery config ntest (arguments result)) >> hFlush stdout - case ok result of - Nothing -> - mytests config gen rnd1 ntest (nfail+1) stamps - Just True -> - mytests config gen rnd1 (ntest+1) nfail (stamp result:stamps) - Just False -> - putStr ( "Falsifiable after " - ++ show ntest - ++ " tests:\n" - ++ unlines (arguments result) - ) >> hFlush stdout >> return (False, ntest) - where - result = generate (configSize config ntest) rnd2 gen - (rnd1,rnd2) = split rnd0 - -done :: String -> Int -> [[String]] -> IO () -done mesg ntest stamps = putStr ( mesg ++ " " ++ show ntest ++ " tests" ++ table ) - where - table = display - . map entry - . reverse - . sort - . map pairLength - . group - . sort - . filter (not . null) - $ stamps - - display [] = ".\n" - display [x] = " (" ++ x ++ ").\n" - display xs = ".\n" ++ unlines (map (++ ".") xs) - - pairLength xss@(xs:_) = (length xss, xs) - entry (n, xs) = percentage n ntest - ++ " " - ++ concat (intersperse ", " xs) - - percentage n m = show ((100 * n) `div` m) ++ "%" - ------------------------------------------------------------------------- - - - main :: IO () -main = do - args <- getArgs - let n = if null args then 100 else read (head args) - - (results, passed) <- liftM unzip $ - foldM ( \prev (name,subtests) -> do - printf "\n%s\n" name - printf "%s\n" $ replicate (length name) '-' - cur <- mapM (\(s,a) -> printf "%-30s: " s >> a n) subtests - return (prev ++ cur) - ) - [] - tests - - printf "\nPassed %d tests!\n\n" (sum passed) - when (not . and $ results) $ fail "\nNot all tests passed!" - where - - tests = [ ("Walker" , tests_Walker) - ] +main = defaultMain [ tests_Walker ]