hs-carbon (empty) → 0.0.0.2

raw patch · 10 files changed

+471/−0 lines, 10 filesdep +basedep +deepseqdep +glosssetup-changed

Dependencies added: base, deepseq, gloss, monad-loops, mtl, parallel, random, tf-random

Files

• LICENSE +7/−0
• Setup.hs +2/−0
• examples/Integral.hs +28/−0
• examples/Pi.hs +24/−0
• examples/Transport/Transport.hs +134/−0
• hs-carbon.cabal +58/−0
• src/Control/Monad/MonteCarlo.hs +124/−0
• src/Data/Result.hs +51/−0
• src/Data/Summary.hs +13/−0
• src/Data/Summary/Bool.hs +30/−0

@@ -0,0 +1,7 @@+Copyright © 2014 Casper M. H. Holmgreen++Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,28 @@+module Main where++import Control.Monad.MonteCarlo+import Data.Summary.Bool+import System.Random.TF++----------------------------------------------------------------+-- Example: Integrate sin(x) from 0 to pi+----------------++bounds :: ((Double,Double),(Double,Double))+bounds = ((0,pi),(0,1))++isUnderCurve :: RandomGen g => (Double -> Double) -> MonteCarlo g Bool+isUnderCurve f = do+    x <- randomR (fst bounds)+    y <- randomR (snd bounds)+    return $ y <= f x++noRuns :: Int+noRuns = 1000000++main :: IO ()+main = do+    g <- newTFGen+    let s = experimentP (isUnderCurve sin) noRuns 200000 g :: BoolSumm+    let ((_,r),(_,u)) = bounds+    print $ sampleMean s * r * u

@@ -0,0 +1,24 @@+module Main where++import Control.Monad (liftM2)+import Control.Monad.MonteCarlo+import Data.Summary.Bool+import System.Random.TF++mcSquareD :: RandomGen g => MonteCarlo g (Double,Double)+mcSquareD = liftM2 (,) (randomR (-1,1)) (randomR (-1,1))++inUnitCircle :: RandomGen g => MonteCarlo g Bool+inUnitCircle = do+    (x,y) <- mcSquareD+    return $ x*x + y*y <= 1++noRuns :: Int+noRuns = 1000000++main :: IO ()+main = do+    g <- newTFGen+    let s = experimentP inUnitCircle noRuns (noRuns `div` 200) g :: BoolSumm+    let (m,se) = (4*sampleMean s, 4*sampleSE s)+    putStrLn $ "Pi is probably between " ++ show (m-se,m+se)

@@ -0,0 +1,134 @@+{-# LANGUAGE BangPatterns #-}++module Main where++import Control.Monad.State+import Control.Monad.Writer+import Control.Monad.Reader+import Control.Monad.MonteCarlo+import Control.Monad.Loops+import Control.DeepSeq+import Control.Exception+import System.Random.TF+import Transport.NISTData+import Data.List (foldl')++import Graphics.Gloss hiding (Point, rotate)++----------------------------------------------------------------+-- Datatypes+----------------++data ParticleState = PS+                     {+                       noColls   :: !Int+                     , remEnergy :: !Energy+                     , curPos    :: !Point+                     , curDir    :: !Angle+                     , path      :: [(Point,Energy)]+                     } deriving (Show)+psInit :: ParticleState+psInit = PS 0 5000 (0,0) (0,1) []++type Energy = Float+type Point = (Float,Float)+type Angle = (Float,Float)++----------------------------------------------------------------+-- MonteCarlo+----------------++type Simulation = ReaderT (Float -> (Float,Float)) (StateT ParticleState (MonteCarlo TFGen))++-- Helper function for getting the cross-section data for the current energy+getMu :: Simulation (Float,Float)+getMu = do+    en <- gets remEnergy+    (t,a) <- asks (\f -> f en)+    return (rho*t,rho*a)+  where rho = 1 -- g/cm^3 (water)++-- Helper functions for sampling random numbers+uniform :: Simulation Float+uniform = lift (lift random)+uniformR :: (Float,Float) -> Simulation Float+uniformR bounds = lift (lift (randomR bounds))++-- Flies the particle some random distance with prob. according to+--  cross-section data+fly :: Simulation ()+fly = do+    (PS i en (x,y) (ux,uy) ps) <- get+    (mu_t,_) <- getMu+    !eta <- uniform+    let s = -(log eta / mu_t)+    put (PS i en (x+ux*s,y+uy*s) (ux,uy) ps)++-- The main loop responsible for a single photon's lifetime+loop :: Simulation [(Point,Energy)]+loop = do+    untilM_ (fly >> scatter) isBelowCutoff+    exit++-- Terminates a particle if its energy is below the cutoff+isBelowCutoff :: Simulation Bool+isBelowCutoff = do+    en <- gets remEnergy+    return $ en < 10++-- Returns the path stored in the ParticleState+exit :: Simulation [(Point,Energy)]+exit = do+    ps <- gets path+    return $ ps++-- Randomly determines whether the scattering event scatters left or right+_scatterDir :: Simulation Float+_scatterDir = do+    eta <- uniform+    return $ if eta >= 0.5 then 1 else (-1)++-- Scattering event; compute scattering angle, record collision site+scatter :: Simulation ()+scatter = do+    (PS i en (x,y) (ux,uy) ps) <- get+    (mu_t,mu_en) <- getMu+    let deltaW = mu_en * en / mu_t+    dir <- _scatterDir+    let angle = diffAngle en (en-deltaW) * dir+    let (ux',uy') = rotate (ux,uy) angle+    put (PS (i+1) (en-deltaW) (x,y) (ux',uy') (ps++[((x,y),deltaW)]))++-- Computes the angle to rotate based on energy exchanged in coll.+diffAngle :: Energy -> Energy -> Float+diffAngle en en' = acos $ 1 - 0.511 * (1/en' - 1/en) -- Knuth++-- Rotates a vector+rotate :: Point -> Float -> Point+rotate (x,y) th = (x*cos th - y*sin th, x*sin th + y*cos th)++----------------------------------------------------------------+-- Main+----------------++noRuns :: Int+noRuns = 10000++main :: IO ()+main = do+    g <- newTFGen+    fnist <- loadData "water.dat"+    let unrolled = evalStateT (runReaderT loop fnist) psInit+    let bs = experimentP (unrolled)+                         noRuns (noRuns `div` 200) g :: [[(Point,Energy)]]+    evaluate (rnf bs)+    let lengthF = fromIntegral . length+    let avgCol = (foldl' (+) 0 (map lengthF bs)) / lengthF bs :: Double+    putStrLn $ "Average number of collisions: " ++ show avgCol+    displayResults bs++displayResults :: [[(Point, Energy)]] -> IO ()+displayResults res = display (InWindow "Sim." (800,800) (200,200))+                         white (results `mappend` Color white (Line [(0,-100),(0,0)]))+  where color' = makeColor8 0 0 0 100+        results = mconcat $ map (\p -> Color color' $ Line (map fst p)) res

@@ -0,0 +1,58 @@+-- Initial hs-carbon.cabal generated by cabal init.  For further +-- documentation, see http://haskell.org/cabal/users-guide/++name:                hs-carbon+version:             0.0.0.2+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.+license:             MIT+license-file:        LICENSE+author:              Casper M. H. Holmgreen+maintainer:          cholmgreen@gmail.com+-- copyright:           +category:            Simulation+build-type:          Simple+cabal-version:       >=1.8++flag buildExamples+    description: Build example executables+    default: False++library+  exposed-modules:     Control.Monad.MonteCarlo+                     , Data.Result+                     , Data.Summary+                     , Data.Summary.Bool+  -- other-modules:       +  build-depends:+    base == 4.*, mtl, random, parallel, deepseq+  hs-source-dirs:      src+  ghc-options:         -Wall++source-repository head+  type:     git+  location: https://github.com/icasperzen/hs-carbon++executable PiExample+  hs-source-dirs: src, examples+  build-depends:+    base == 4.*, random, mtl, tf-random, parallel, deepseq+  main-is: Pi.hs+  ghc-options: -Wall -threaded -O3++executable IntegralExample+  hs-source-dirs: src, examples+  build-depends:+    base == 4.*, random, mtl, tf-random, parallel, deepseq+  main-is: Integral.hs+  ghc-options: -Wall -threaded -O3++executable TransportExample+  hs-source-dirs: src, examples+  build-depends:+      base == 4.*, random, mtl, tf-random, parallel, gloss, monad-loops, deepseq+  main-is: Transport/Transport.hs+  ghc-options: -Wall -threaded -O3

@@ -0,0 +1,124 @@+{-# LANGUAGE TypeFamilies, BangPatterns #-}++-----------------------------------------------------------------------------+-- |+-- Module      :  Control.Monad.MonteCarlo+-- Copyright   :  (c) Casper M. H. Holmgreen+-- License     :  MIT+-- +-- Maintainer  :  cholmgreen@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- The goal of Carbon is to capture idiomatic usage patterns for Monte Carlo simulations in an easy to use interface.+-- Currently, only one pattern is implemented but more are planned.+--+-- There are two main parts to the library:+--+-- * The 'MonteCarlo' monad for building simulations+--+-- * The 'Data.Result' type family for describing how to aggregate results+--+-- Running a simulation will require a basic understanding of how to use both.+-- Control.Monad.MonteCarlo exports the 'System.Random.RandomGen' typeclass for convenience.+--+-----------------------------------------------------------------------------++module Control.Monad.MonteCarlo+  (+    -- * Running Simulations+    -- $runningsims+    MonteCarlo+  , experimentS+  , experimentP+    -- * Building Simulations+    -- $buildingsims+  , runMC+  , evalMC+  , random+  , randomR+  , R.RandomGen -- export RandomGen for convenience+  )+where++import Control.Monad.State+import Control.Parallel.Strategies+import Data.List (foldl')+import Data.Result+import qualified System.Random as R++{- $runningsims #runningsims#+Running simulations should be done using the high level functions provided here.+Monte Carlo simulations are specified by two parts: a 'MonteCarlo' action and a type annotation describing how to aggregate the results.+See 'Data.Result' for more information.+-}++{- $buildingsims #buildingsims#+Simulations are built by structuring 'MonteCarlo' actions over top of eachother, in a manner similar to parser combinators.+This compositional approach makes it very easy to express complex simulations without sacrificing readability.+Additionally, more complex simulations can layer monad transformers overtop of 'MonteCarlo' to add support for State, RO-environments, etc.+-}++-- | The 'MonteCarlo' monad is just a 'Control.Monad.State' monad in disguise.+-- This allows us to thread the internal PRNG state through the simulation as we sample new numbers.+--+type MonteCarlo g = State g++-- | This is a high level function for running a full Monte Carlo simulation.+-- It takes a 'MonteCarlo' action, the number of observations to aggregate, and an instance of 'System.Random.RandomGen'.+-- The return value is dictated by the type family 'Data.Result'; a type annotation is required to specify how observations should be aggregated.+--+-- For example, given a 'MonteCarlo' action, mySim, with type:+--+-- > mySim :: RandomGen g => MonteCarlo g Bool+--+-- We can get radically different results with just a type annotation:+--+-- > experimentS mySimulation 100 g :: [Bool]+-- > experimentS mySimulation 100 g :: BoolSumm+--+experimentS :: (R.RandomGen g, Result s)+            => MonteCarlo g (Obs s) -> Int -> g -> s+experimentS m n g = let xs = evalMC (replicateM n m) g+                     in foldl' addObs rzero xs++-- | This is a high level function for running a full Monte Carlo simulation in parallel.+-- It is identical to 'experimentS' except that it takes an additional Int argument representing the chunk size.+-- Determining a good chunk size is an art, but a good starting point may be to divide the number of runs by 200.+--+-- Note: you must compile an executable with the -threaded flag in order for sparks to run in parallel.+--+experimentP :: (R.RandomGen g, Result s)+            => 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)++-- | 'runMC' is an alias for 'runState'.+runMC :: R.RandomGen g => MonteCarlo g a -> g -> (a,g)+runMC = runState++-- | 'evalMC' is an alias for 'evalState'.+evalMC :: R.RandomGen g => MonteCarlo g a -> g -> a+evalMC = evalState++-- | 'mcNext' is a higher-order function which runs typical "System.Random" functions and updates the internal state.+mcNext :: R.RandomGen g => (g -> (a,g)) -> MonteCarlo g a+mcNext f = do+    !g <- get+    let !(!x,!g') = f g+    put g'+    return x++-- | 'random' calls 'System.Random.random' and updates the internal state+random :: (R.RandomGen g, R.Random a) => MonteCarlo g a+random = mcNext R.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)

@@ -0,0 +1,51 @@+{-# LANGUAGE TypeFamilies #-}++module Data.Result+  (+    -- * Describing Simulation Results with Type Annotations+    -- $description+    Result(..)+  )+where++import Control.DeepSeq+import Data.Monoid++{- $description #description#+Carbon simulations are built up from 'Control.Monad.MonteCarlo' actions.+A 'MonteCarlo' action describes how to arrive at an observation, but not how to aggregate observations.+This functionality is specified with a type annotation telling Haskell which instance of the type family 'Result' should be used.++For example, given a 'MonteCarlo' action, mySim, with type:++> mySim :: RandomGen g => MonteCarlo g Bool++We get different results based on the instance of 'Result' chosen:++> experimentS mySimulation 100 g :: [Bool]+> experimentS mySimulation 100 g :: BoolSumm+-}++-- | Result is the type family used to describe the aggregation techniques to be used in a Monte Carlo simulation.+-- Instances of Result should specify the type of a single observation and how to include one.+-- The value of a 'Result' without any observations should be specified.+-- Additionally, 'Result's should be joinable.+--+-- Note that almost all instances of 'Result' will be monoidal.+class Result r where+    -- | The type of a single observation+    type Obs r :: *+    -- | How to add a single observation to the result+    addObs :: r -> Obs r -> r+    -- | How to join two results together+    rjoin  :: r -> r -> r+    -- | The value of a result without any observations+    rzero  :: r++-- | One common aggregation technique is appending observations to a list.+-- In the interest of preventing thunks from building up, we force deep evaluation of observations.+instance NFData a => Result [a] where+    type Obs [a] = a+    addObs r o   = o `deepseq` (o : r)+    rjoin        = mappend+    rzero        = mempty

@@ -0,0 +1,13 @@+{-# LANGUAGE TypeFamilies #-}++module Data.Summary where++-- | Many Monte Carlo simulations require statistical analysis of the results.+-- Any 'Data.Result' instances which can be described statistically should be made instances of 'Summary'.+class Summary s where+    -- | Compute the mean of the aggregated observations+    sampleMean :: s -> Double+    -- | Compute the std. error of the aggregated observations+    sampleSE   :: s -> Double+    -- | Return the number of observations aggregated+    sampleSize :: s -> Int

@@ -0,0 +1,30 @@+{-# LANGUAGE TypeFamilies #-}++module Data.Summary.Bool+  (BoolSumm, Summary(..))+  where++import Data.Result (Result(..))+import Data.Summary (Summary(..))++-- | A 'BoolSumm' counts the number of True and all events observed.+data BoolSumm = BoolSumm+                 {+                   _noSuccess :: !Int+                 , _noTotal   :: !Int+                 }++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))+    rjoin (BoolSumm s t) (BoolSumm s' t') = BoolSumm (s+s') (t+t')+    rzero = BoolSumm 0 0++instance Summary BoolSumm where+    sampleMean (BoolSumm s t) = fromIntegral s / fromIntegral t+    sampleSE s = sqrt (p * (1 - p) / n)+      where+        p = sampleMean s+        n = fromIntegral $ sampleSize s+    sampleSize (BoolSumm _ t) = t