hasty-hamiltonian

commit 49f111bb128bb4204af8d73dc6d6888a98bbfab9
parent 42eb070584628f1bb605b65d0332815976d72c05
Author: Jared Tobin <jared@jtobin.ca>
Date:   Mon,  5 Nov 2012 13:27:26 +1300

Serialize to float rather than double.

Diffstat:
M
Examples/code/Rosenbrock_HMC.hs
 | 
2
+-
A
Examples/code/SPDE_HMC.hs
 | 
55
+++++++++++++++++++++++++++++++++++++++++++++++++++++++
M
Numeric/MCMC/Hamiltonian.hs
 | 
3
++-

3 files changed, 58 insertions(+), 2 deletions(-)
diff --git a/Examples/code/Rosenbrock_HMC.hs b/Examples/code/Rosenbrock_HMC.hs
@@ -26,7 +26,7 @@ main = do
         putStrLn  "                   initialize the chain.                                 "
         exitSuccess
 
-    inits <- fmap (map read . words) (readFile (args !! 3)) :: IO [Double]
+    inits <- fmap (map read . words) (readFile (args !! 4)) :: IO [Double]
 
     let nepochs   = read (head args) :: Int
         thinEvery = read (args !! 1) :: Int
diff --git a/Examples/code/SPDE_HMC.hs b/Examples/code/SPDE_HMC.hs
@@ -0,0 +1,55 @@
+{-# LANGUAGE BangPatterns #-}
+
+import Numeric.MCMC.Hamiltonian
+import System.Random.MWC
+import System.Environment
+import System.Exit
+import System.IO
+import Control.Monad
+import Numeric.AD
+
+target :: RealFloat a => [a] -> a
+target xs = go 0 0 xs 
+  where go !t0 !t1 []         = (- t0 / (2*h)) - (0.5 * h * t1)
+        go !t0 !t1 (u:us:uss) = go (t0 + (us - u)^2) (t1 + v (us + u)) uss
+        h   = 1 / fromIntegral (length xs)
+        v x = (1 - x^2)^2
+{-# INLINE target #-}
+
+gTarget :: [Double] -> [Double]
+gTarget = grad target
+{-# INLINE gTarget #-}
+
+main = do
+    args  <- getArgs 
+    when (args == []) $ do
+        putStrLn  "(hasty-hamiltonian) Stochastic partial differential equation "
+        putStrLn  "Usage: ./SPDE_HMC <numSteps> <inits> <thinEvery>             " 
+        putStrLn  "                                                             "
+        putStrLn  "numSteps         : Number of Markov chain iterations to run.             "
+        putStrLn  "thinEvery        : Print every nth iteration.                            "
+        putStrLn  "nDisc            : Number of discretizing steps to take.                 "
+        putStrLn  "stepSize         : Perturbation scaling parameter.                       "
+        putStrLn  "inits            : Filepath containing points at which to                "
+        putStrLn  "                   initialize the chain.                                 "
+        exitSuccess
+
+    inits <- fmap (map read . words) (readFile (args !! 4)) :: IO [Double]
+
+    let nepochs   = read (head args) :: Int
+        thinEvery = read (args !! 1) :: Int
+        nDisc     = read (args !! 2) :: Double
+        eps       = read (args !! 3) :: Double
+        params    = Options target gTarget nDisc eps
+        config    = MarkovChain inits 0
+
+    g       <- create
+    results <- runChain params nepochs thinEvery config g
+
+    hPutStrLn stderr $ 
+        let nAcc  = accepts results
+            total = nepochs 
+        in  show nAcc ++ " / " ++ show total ++ " (" ++ 
+              show ((fromIntegral nAcc / fromIntegral total) :: Float) ++ 
+              ") proposals accepted"
+
diff --git a/Numeric/MCMC/Hamiltonian.hs b/Numeric/MCMC/Hamiltonian.hs
@@ -7,6 +7,7 @@ import Control.Monad.Reader
 import Control.Monad.Primitive
 import System.Random.MWC
 import System.Random.MWC.Distributions
+import GHC.Float
 
 -- | State of the Markov chain.  Current parameter values are held in 'theta', 
 --   while accepts counts the number of proposals accepted.
@@ -25,7 +26,7 @@ type ViewsOptions = ReaderT Options
 
 -- | Display the current state. 
 instance Show MarkovChain where
-    show config = filter (`notElem` "[]") $ show (theta config)
+    show config = filter (`notElem` "[]") $ show (map double2Float (theta config))
 
 -- | The 'leapfrog' or Stormer-Verlet discretizer.
 leapfrog :: Monad m