monte-carlo 0.2 → 0.3
raw patch · 23 files changed
+1043/−654 lines, 23 filesdep +vectordep −arraydep −uvectordep ~basedep ~gsl-randomdep ~mtl
Dependencies added: vector
Dependencies removed: array, uvector
Dependency ranges changed: base, gsl-random, mtl
Files
- LICENSE +1/−1
- NEWS +40/−0
- examples/Binomial.hs +55/−0
- examples/Pi.hs +0/−82
- examples/Pi.lhs +113/−0
- examples/Poker.hs +23/−24
- examples/Queue.hs +192/−0
- examples/Sampling.hs +0/−58
- lib/Control/Monad/MC.hs +6/−2
- lib/Control/Monad/MC/Base.hs +36/−9
- lib/Control/Monad/MC/Class.hs +7/−14
- lib/Control/Monad/MC/GSL.hs +11/−3
- lib/Control/Monad/MC/GSLBase.hs +139/−121
- lib/Control/Monad/MC/Repeat.hs +14/−32
- lib/Control/Monad/MC/Sample.hs +106/−75
- lib/Control/Monad/MC/Summary.hs +0/−96
- lib/Control/Monad/MC/Walker.hs +31/−30
- lib/Data/Summary.hs +15/−0
- lib/Data/Summary/Bool.hs +86/−0
- lib/Data/Summary/Double.hs +107/−0
- lib/Data/Summary/Utils.hs +36/−0
- monte-carlo.cabal +17/−13
- tests/Main.hs +8/−94
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) Patrick Perry <patperry@stanford.edu> 2008+Copyright (c) Patrick Perry <patperry@gmail.com> 2010 All rights reserved.
+ NEWS view
@@ -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+
+ examples/Binomial.hs view
@@ -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)
− examples/Pi.hs
@@ -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
+ examples/Pi.lhs view
@@ -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}
examples/Poker.hs view
@@ -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"
+ examples/Queue.hs view
@@ -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 ""
− examples/Sampling.hs
@@ -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)
lib/Control/Monad/MC.hs view
@@ -1,10 +1,14 @@ ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com> -- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Maintainer : Patrick Perry <patperry@gmail.com> -- 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 (
lib/Control/Monad/MC/Base.hs view
@@ -1,20 +1,17 @@-{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, FlexibleContexts #-}+{-# LANGUAGE TypeFamilies #-} ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Base--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com> -- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Maintainer : Patrick Perry <patperry@gmail.com> -- 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
lib/Control/Monad/MC/Class.hs view
@@ -1,32 +1,25 @@ ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Class--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com> -- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Maintainer : Patrick Perry <patperry@gmail.com> -- 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
lib/Control/Monad/MC/GSL.hs view
@@ -1,11 +1,13 @@ ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.GSL--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com> -- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Maintainer : Patrick Perry <patperry@gmail.com> -- 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 )
lib/Control/Monad/MC/GSLBase.hs view
@@ -1,10 +1,10 @@-{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-} ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.GSLBase--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com> -- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Maintainer : Patrick Perry <patperry@gmail.com> -- 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
lib/Control/Monad/MC/Repeat.hs view
@@ -1,49 +1,31 @@ ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Repeat--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com> -- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Maintainer : Patrick Perry <patperry@gmail.com> -- 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 []
lib/Control/Monad/MC/Sample.hs view
@@ -1,10 +1,9 @@-{-# LANGUAGE ScopedTypeVariables #-} ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Sample--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com> -- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Maintainer : Patrick Perry <patperry@gmail.com> -- 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' #-}
− lib/Control/Monad/MC/Summary.hs
@@ -1,96 +0,0 @@--------------------------------------------------------------------------------- |--- Module : Control.Monad.MC.Summary--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>--- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>--- 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
lib/Control/Monad/MC/Walker.hs view
@@ -1,10 +1,9 @@-{-# LANGUAGE TypeOperators #-} ----------------------------------------------------------------------------- -- | -- Module : Control.Monad.MC.Walker--- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com> -- License : BSD3--- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Maintainer : Patrick Perry <patperry@gmail.com> -- 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)
+ lib/Data/Summary.hs view
@@ -0,0 +1,15 @@+-----------------------------------------------------------------------------+-- |+-- Module : Data.Summary+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@gmail.com>+-- Stability : experimental+--+-- Summary Statistics+--+module Data.Summary (+ module Data.Summary.Double+ ) where++import Data.Summary.Double
+ lib/Data/Summary/Bool.hs view
@@ -0,0 +1,86 @@+-----------------------------------------------------------------------------+-- |+-- Module : Data.Summary.Bool+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@gmail.com>+-- 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)
+ lib/Data/Summary/Double.hs view
@@ -0,0 +1,107 @@+-----------------------------------------------------------------------------+-- |+-- Module : Data.Summary.Double+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@gmail.com>+-- 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
+ lib/Data/Summary/Utils.hs view
@@ -0,0 +1,36 @@+-----------------------------------------------------------------------------+-- |+-- Module : Data.Summary.Utils+-- Copyright : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@gmail.com>+-- 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
monte-carlo.cabal view
@@ -1,10 +1,10 @@ name: monte-carlo-version: 0.2+version: 0.3 license: BSD3 license-file: LICENSE author: Patrick Perry-maintainer: Patrick Perry <patperry@stanford.edu>-homepage: http://github.com/patperry/monte-carlo+maintainer: Patrick Perry <patperry@gmail.com>+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
tests/Main.hs view
@@ -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 ]