hs-carbon 0.0.0.3 → 0.0.1.0
raw patch · 4 files changed
+1/−217 lines, 4 filesdep −glossdep −monad-loopsdep −tf-randomPVP ok
version bump matches the API change (PVP)
Dependencies removed: gloss, monad-loops, tf-random
API changes (from Hackage documentation)
Files
- examples/Integral.hs +0/−28
- examples/Pi.hs +0/−24
- examples/Transport/Transport.hs +0/−134
- hs-carbon.cabal +1/−31
− examples/Integral.hs
@@ -1,28 +0,0 @@-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
− examples/Pi.hs
@@ -1,24 +0,0 @@-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)
− examples/Transport/Transport.hs
@@ -1,134 +0,0 @@-{-# 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
hs-carbon.cabal view
@@ -1,8 +1,6 @@--- 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.3+version: 0.0.1.0 synopsis: A Haskell framework for parallel monte carlo simulations description: hs-carbon is a PRNG-agnostic Haskell framework for running monte-carlo@@ -17,10 +15,6 @@ build-type: Simple cabal-version: >=1.8 -flag buildExamples- description: Build example executables- default: False- library exposed-modules: Control.Monad.MonteCarlo , Data.Result@@ -35,27 +29,3 @@ source-repository head type: git location: https://github.com/icasperzen/hs-carbon--executable PiExample- hs-source-dirs: src, examples- buildable: False- 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- buildable: False- 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- buildable: False- build-depends:- base == 4.*, random, mtl, tf-random, parallel, gloss, monad-loops, deepseq- main-is: Transport/Transport.hs- ghc-options: -Wall -threaded -O3