commit c1b2e366f4928ea9747dce09097f07af9d6d0600
parent e47c49ce4e7d250e6f618d34858beba3fd0b5d9f
Author: Jared Tobin <jared@jtobin.ca>
Date: Mon, 5 Nov 2012 09:18:23 +1300
Add thinning parameter.
Diffstat:
6 files changed, 39 insertions(+), 16 deletions(-)
diff --git a/Examples/code/BNN_Flat.hs b/Examples/code/BNN_Flat.hs
@@ -28,7 +28,7 @@ main = do
config = MarkovChain inits 0
g <- create
- results <- runChain params nepochs config g
+ results <- runChain params nepochs 1 config g
hPutStrLn stderr $
let nAcc = accepts results
diff --git a/Examples/code/Booth_Flat.hs b/Examples/code/Booth_Flat.hs
@@ -28,7 +28,7 @@ main = do
config = MarkovChain inits 0
g <- create
- results <- runChain params nepochs config g
+ results <- runChain params nepochs 1 config g
hPutStrLn stderr $
let nAcc = accepts results
diff --git a/Examples/code/Himmelblau_Flat.hs b/Examples/code/Himmelblau_Flat.hs
@@ -28,7 +28,7 @@ main = do
config = MarkovChain inits 0
g <- create
- results <- runChain params nepochs config g
+ results <- runChain params nepochs 1 config g
hPutStrLn stderr $
let nAcc = accepts results
diff --git a/Examples/code/Rosenbrock_Flat.hs b/Examples/code/Rosenbrock_Flat.hs
@@ -28,7 +28,7 @@ main = do
config = MarkovChain inits 0
g <- create
- results <- runChain params nepochs config g
+ results <- runChain params nepochs 1 config g
hPutStrLn stderr $
let nAcc = accepts results
diff --git a/Examples/code/SPDE_Flat.hs b/Examples/code/SPDE_Flat.hs
@@ -18,23 +18,25 @@ main = do
args <- getArgs
when (args == []) $ do
putStrLn "(flat-mcmc) Stochastic partial differential equation "
- putStrLn "Usage: ./SPDE_Flat <numSteps> <inits> <granularity> "
+ putStrLn "Usage: ./SPDE_Flat <numSteps> <inits> <thinEvery> <granularity> "
putStrLn " "
putStrLn "numSteps : Number of Markov chain iterations to run."
putStrLn "inits : Filepath containing points at which to "
putStrLn " initialize the ensemble. "
+ putStrLn "thinEvery : Print every n^th iteration. "
putStrLn "granularity : Parallel granularity (smaller is finer). "
exitSuccess
inits <- readInits (args !! 1)
- let nepochs = read (head args) :: Int
- gran = read (args !! 2) :: Int
- params = Options target (V.length inits) gran
- config = MarkovChain inits 0
+ let nepochs = read (head args) :: Int
+ thinEvery = read (args !! 2) :: Int
+ gran = read (args !! 3) :: Int
+ params = Options target (V.length inits) gran
+ config = MarkovChain inits 0
g <- create
- results <- runChain params nepochs config g
+ results <- runChain params nepochs thinEvery config g
hPutStrLn stderr $
let nAcc = accepts results
diff --git a/Numeric/MCMC/Flat.hs b/Numeric/MCMC/Flat.hs
@@ -1,4 +1,5 @@
{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE BangPatterns #-}
module Numeric.MCMC.Flat (
MarkovChain(..), Options(..), Ensemble
@@ -15,6 +16,7 @@ import qualified Data.Vector.Unboxed as U
import Control.Monad.Par (NFData)
import Control.Monad.Par.Scheds.Direct
import Control.Monad.Par.Combinator
+import System.IO
-- | Parallel map with a specified granularity.
parMapChunk :: NFData b => Int -> (a -> b) -> [a] -> Par [b]
@@ -53,6 +55,7 @@ symmetricVariate :: PrimMonad m => Gen (PrimState m) -> m Double
symmetricVariate g = do
z <- uniformR (0 :: Double, 1 :: Double) g
return $! 0.5*(z + 1)^(2 :: Int)
+{-# INLINE symmetricVariate #-}
-- | The result of a single-particle Metropolis accept/reject step. This
-- compares a particle's state to a perturbation made by an affine
@@ -66,6 +69,7 @@ metropolisResult w0 w1 z zc target =
let val = target proposal - target w0 + (fromIntegral (length w0) - 1) * log z
proposal = zipWith (+) (map (*z) w0) (map (*(1-z)) w1)
in if zc <= min 1 (exp val) then (proposal, 1) else (w0, 0)
+{-# INLINE metropolisResult #-}
-- | Execute Metropolis steps on the particles of a sub-ensemble by
-- perturbing them with affine transformations based on particles
@@ -92,6 +96,7 @@ executeMoves e0 e1 n g = do
(newstate, nacc) = (V.fromList . map fst &&& sum . map snd) result
return (newstate, nacc)
+{-# INLINE executeMoves #-}
-- | Perform a Metropolis accept/reject step on the ensemble by
-- perturbing each element and accepting/rejecting the perturbation in
@@ -113,23 +118,39 @@ metropolisStep state g = do
return $!
MarkovChain (V.concat $ map fst [result0, result1])
(nacc + snd result0 + snd result1)
+{-# INLINE metropolisStep #-}
-- | Diffuse through states.
runChain :: Options -- Options of the Markov chain
-> Int -- Number of epochs to iterate the chain
+ -> Int -- Print every nth iteration.
-> MarkovChain -- Initial state of the Markov chain
-> Gen RealWorld -- MWC PRNG
-> IO MarkovChain -- End state of the Markov chain, wrapped in IO
-runChain params nepochs initConfig g
- | nepochs == 0 = return initConfig
- | otherwise = do
- result <- runReaderT (metropolisStep initConfig g) params
- print result
- runChain params (nepochs - 1) result g
+runChain opts nepochs thinEvery initConfig g
+ | l == 0
+ = error "runChain: ensemble must contain at least one particle"
+ | l < (length . V.head) (ensemble initConfig)
+ = do hPutStrLn stderr $ "runChain: ensemble should be twice as large as "
+ ++ "the target's dimension. Continuing anyway."
+ go opts nepochs thinEvery initConfig g
+ | otherwise = go opts nepochs thinEvery initConfig g
+ where
+ l = V.length (ensemble initConfig)
+ go o n t !c g0 | n == 0 = return c
+ | n `rem` t /= 0 = do
+ r <- runReaderT (metropolisStep c g0) o
+ go o (n - 1) t r g0
+ | otherwise = do
+ r <- runReaderT (metropolisStep c g0) o
+ print r
+ go o (n - 1) t r g0
+{-# INLINE runChain #-}
-- | A convenience function to read and parse ensemble inits from disk.
-- Assumes a text file with one particle per line, where each particle
-- element is separated by whitespace.
readInits :: FilePath -> IO Ensemble
readInits p = fmap (V.fromList . map (map read . words) . lines) (readFile p)
+{-# INLINE readInits #-}
