hs-carbon-examples (empty) → 0.0.0.1

raw patch · 6 files changed

+233/−0 lines, 6 filesdep +basedep +deepseqdep +glosssetup-changed

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

Files

• LICENSE +7/−0
• Setup.hs +2/−0
• examples-src/Integral.hs +28/−0
• examples-src/Pi.hs +24/−0
• examples-src/Transport/Transport.hs +134/−0
• hs-carbon-examples.cabal +38/−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,38 @@++name:                hs-carbon-examples+version:             0.0.0.1+synopsis:            Example Monte Carlo simulations implemented with Carbon+description:         Example Monte Carlo simulations implemented with Carbon+license:             MIT+license-file:        LICENSE+author:              Casper M. H. Holmgreen+maintainer:          cholmgreen@gmail.com+-- copyright:           +category:            Simulation+build-type:          Simple+cabal-version:       >=1.8++source-repository head+  type:     git+  location: https://github.com/icasperzen/hs-carbon-examples++executable PiExample+  hs-source-dirs: examples-src+  build-depends:+    base == 4.*, hs-carbon, tf-random+  main-is: Pi.hs+  ghc-options: -Wall -threaded -O3++executable IntegralExample+  hs-source-dirs: examples-src+  build-depends:+    base == 4.*, hs-carbon, tf-random+  main-is: Integral.hs+  ghc-options: -Wall -threaded -O3++executable TransportExample+  hs-source-dirs: examples-src+  build-depends:+      base == 4.*, hs-carbon, gloss, tf-random, deepseq, monad-loops, mtl+  main-is: Transport/Transport.hs+  ghc-options: -Wall -threaded -O3