diff --git a/hs-carbon.cabal b/hs-carbon.cabal --- a/hs-carbon.cabal +++ b/hs-carbon.cabal @@ -1,11 +1,13 @@ name: hs-carbon -version: 0.0.1.0 +version: 0.1.0.0 synopsis: A Haskell framework for parallel monte carlo simulations description: - hs-carbon is a PRNG-agnostic Haskell framework for running monte-carlo - simulations. The library will provide several "skeletons" for abstracting - away common usage patterns. + Carbon is an open-source, Haskell framework aiming to provide easy access to + parallel Monte Carlo simulations by providing a simple, but powerful + compositional method for building simulations and high-level functions for + running them. + Examples can be found at https://github.com/icasperzen/hs-carbon-examples license: MIT license-file: LICENSE author: Casper M. H. Holmgreen @@ -20,11 +22,19 @@ , Data.Result , Data.Summary , Data.Summary.Bool + , Data.Summary.Double -- other-modules: build-depends: base == 4.*, mtl, random, parallel, deepseq hs-source-dirs: src ghc-options: -Wall + +Test-Suite tests + type: exitcode-stdio-1.0 + main-is: Summary.hs + hs-source-dirs: spec + build-depends: base, HUnit, hs-carbon + ghc-options: -Wall source-repository head type: git diff --git a/spec/Summary.hs b/spec/Summary.hs new file mode 100644 --- /dev/null +++ b/spec/Summary.hs @@ -0,0 +1,50 @@ +module Main where + +import Test.HUnit +import System.Exit + +import Data.Summary.Bool +import Data.Summary.Double + +main :: IO () +main = do + c <- runTestTT allTests + case failures c + errors c of + 0 -> exitSuccess + _ -> exitFailure + +allTests :: Test +allTests = TestList [boolTests, doubleTests] + +t :: String -> Assertion -> Test +t cs a = TestLabel cs $ TestCase a + +---------------------------------------------------------------- +-- Data.Summary.Bool +---------------- + +bs :: BoolSumm +bs = boolSumm [True, False, True, False] + +boolTests :: Test +boolTests = TestLabel "Data.Summary.Bool" $ TestList [ + t "sampleMean" $ sampleMean bs @?= 0.5 + , t "sampleSE" $ sampleSE bs @?= 0.25 + , t "sampleSize" $ sampleSize bs @?= 4 + ] + +---------------------------------------------------------------- +-- Data.Summary.Double +---------------- + +ds :: DoubleSumm +ds = doubleSumm [1..5] + +doubleTests :: Test +doubleTests = TestLabel "Data.Summary.Double" $ TestList [ + TestCase $ sampleMean ds @?= 3 + , TestCase $ sampleVar ds @?= 2.5 + , TestCase $ sampleSD ds @?= sqrt 2.5 + , TestCase $ sampleSE ds @?= sqrt 2.5 / sqrt 5 + , TestCase $ sampleSize ds @?= 5 + ] diff --git a/src/Control/Monad/MonteCarlo.hs b/src/Control/Monad/MonteCarlo.hs --- a/src/Control/Monad/MonteCarlo.hs +++ b/src/Control/Monad/MonteCarlo.hs @@ -92,12 +92,11 @@ => MonteCarlo g (Obs s) -> Int -> Int -> g -> s experimentP m n c g | c <= 0 = error "Chunk size must be positive" - | n <= c = experimentS m n g - | otherwise = runEval $ do - let !(!g1,!g2) = R.split g - s <- rpar $ experimentS m c g1 - ss <- rpar $ experimentP m (n-c) c g2 - return (s `rjoin` ss) + | otherwise = let n' = n `div` c + f = experimentS m c + mkGens no seed = (take no $ tail $ map fst $ iterate (\(_,g') -> R.split g') (undefined,seed)) + es = (map f (mkGens n' g) `using` parList rseq) + in foldl' rjoin rzero es -- | 'runMC' is an alias for 'runState'. runMC :: R.RandomGen g => MonteCarlo g a -> g -> (a,g) @@ -114,11 +113,14 @@ let !(!x,!g') = f g put g' return x +{-# INLINE mcNext #-} -- | 'random' calls 'System.Random.random' and updates the internal state random :: (R.RandomGen g, R.Random a) => MonteCarlo g a random = mcNext R.random +{-# INLINE random #-} -- | 'randomR' calls 'System.Random.randomR' and updates the internal state randomR :: (R.RandomGen g, R.Random a) => (a,a) -> MonteCarlo g a -randomR !bounds = mcNext (R.randomR bounds) +randomR !(!l,!u) = mcNext (R.randomR (l,u)) +{-# INLINE randomR #-} diff --git a/src/Data/Summary.hs b/src/Data/Summary.hs --- a/src/Data/Summary.hs +++ b/src/Data/Summary.hs @@ -9,5 +9,9 @@ sampleMean :: s -> Double -- | Compute the std. error of the aggregated observations sampleSE :: s -> Double + -- | Compute the variance of the aggregated observations + sampleVar :: s -> Double + -- | Compute the standard deviation of the aggregated observations + sampleSD :: s -> Double -- | Return the number of observations aggregated sampleSize :: s -> Int diff --git a/src/Data/Summary/Bool.hs b/src/Data/Summary/Bool.hs --- a/src/Data/Summary/Bool.hs +++ b/src/Data/Summary/Bool.hs @@ -1,23 +1,29 @@ {-# LANGUAGE TypeFamilies #-} module Data.Summary.Bool - (BoolSumm, Summary(..)) + (BoolSumm, Summary(..), boolSumm) where import Data.Result (Result(..)) import Data.Summary (Summary(..)) +import Data.List (foldl') +import Control.DeepSeq (NFData(..)) -- | A 'BoolSumm' counts the number of True and all events observed. -data BoolSumm = BoolSumm - { - _noSuccess :: !Int - , _noTotal :: !Int - } +data BoolSumm = BoolSumm { + _noSuccess :: !Int + , _noTotal :: !Int + } deriving (Show) +instance NFData BoolSumm + +boolSumm :: [Bool] -> BoolSumm +boolSumm = foldl' addObs rzero + instance Result BoolSumm where type Obs BoolSumm = Bool - addObs (BoolSumm s t) True = (BoolSumm (s+1) (t+1)) - addObs (BoolSumm s t) False = (BoolSumm s (t+1)) + addObs (BoolSumm s t) True = BoolSumm (s+1) (t+1) + addObs (BoolSumm s t) False = BoolSumm s (t+1) rjoin (BoolSumm s t) (BoolSumm s' t') = BoolSumm (s+s') (t+t') rzero = BoolSumm 0 0 @@ -28,3 +34,9 @@ p = sampleMean s n = fromIntegral $ sampleSize s sampleSize (BoolSumm _ t) = t + sampleSD = error ("sampleSD" ++ undefBinObs) + sampleVar = error ("sampleVar" ++ undefBinObs) + +undefBinObs :: String +undefBinObs = " is undefined for binary observations. Please contact" + ++ " the package maintainer if you can define it." diff --git a/src/Data/Summary/Double.hs b/src/Data/Summary/Double.hs new file mode 100644 --- /dev/null +++ b/src/Data/Summary/Double.hs @@ -0,0 +1,44 @@ +{-# LANGUAGE TypeFamilies #-} + +module Data.Summary.Double where + +import Data.Result (Result(..)) +import Data.Summary (Summary(..)) +import Data.List (foldl') +import Control.DeepSeq (NFData(..)) + +-- | Computes running stats as demonstrated by +-- http://www.johndcook.com/skewness_kurtosis.html +data DoubleSumm = DoubleSumm { + _m1 :: !Double + , _m2 :: !Double + , _size :: !Int +} deriving (Show) + +instance NFData DoubleSumm + +doubleSumm :: [Double] -> DoubleSumm +doubleSumm = foldl' addObs rzero + +instance Result DoubleSumm where + type Obs DoubleSumm = Double + addObs (DoubleSumm m v n) x = DoubleSumm m' v' n' + where + delta = x - m + n' = n + 1 + m' = m + delta/fromIntegral n' + v' = v + delta*(delta/fromIntegral n')*fromIntegral n + rjoin (DoubleSumm m1 v1 n1) (DoubleSumm m2 v2 n2) = DoubleSumm m' v' n' + where + delta = m2 - m1 + n' = n1 + n2 + m' = (fromIntegral n1*m1 + fromIntegral n2*m2)/fromIntegral n' + v' = v1+v2 + delta*delta*fromIntegral n1/fromIntegral n'*fromIntegral n2 + rzero = DoubleSumm 0 0 0 + +instance Summary DoubleSumm where + sampleMean = _m1 + sampleVar ds = _m2 ds / fromIntegral (_size ds - 1) + sampleSD = sqrt . sampleVar + sampleSE ds = sampleSD ds / sqrt (fromIntegral (_size ds)) + sampleSize = _size