mcmc 0.7.0.0 → 0.7.0.1

raw patch · 14 files changed

+51/−52 lines, 14 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• ChangeLog.md +1/−1
• bench/Bench.hs +1/−1
• bench/Normal.hs +4/−4
• bench/Poisson.hs +2/−2
• mcmc.cabal +1/−1
• src/Mcmc/Algorithm.hs +2/−2
• src/Mcmc/Algorithm/MC3.hs +13/−18
• src/Mcmc/Algorithm/MHG.hs +9/−5
• src/Mcmc/Chain/Trace.hs +3/−0
• src/Mcmc/Internal/Shuffle.hs +2/−2
• src/Mcmc/MarginalLikelihood.hs +6/−7
• src/Mcmc/Mcmc.hs +1/−2
• src/Mcmc/Proposal.hs +2/−2
• test/Mcmc/SaveSpec.hs +4/−5

@@ -5,7 +5,7 @@ ## Unreleased changes  -## 0.7.0.0+## 0.7.0.1  -   Use random-1.2. This means, that the random number generator type has changed     from `GenIO` (`System.Random.MWC`) to `IOGenM StdGen`

@@ -30,7 +30,7 @@  main :: IO () main = do-  g <- newIOGenM $ mkStdGen 0+  let g = mkStdGen 0   defaultMain     [ bgroup         "Normal"

@@ -42,7 +42,7 @@ mon :: Monitor Double mon = Monitor monStd [] [] -normalSlideBench :: IOGenM StdGen -> IO ()+normalSlideBench :: StdGen -> IO () normalSlideBench g = do   let s =         Settings@@ -64,7 +64,7 @@     [slideSymmetric 1.0 (PName $ "Medium " ++ show i) (pWeight 1) Tune | i <- [0 .. 100 :: Int]]  -- Should have the same run time as 'normalSlide'.-normalLargeCycleBench :: IOGenM StdGen -> IO ()+normalLargeCycleBench :: StdGen -> IO () normalLargeCycleBench g = do   let s =         Settings@@ -83,7 +83,7 @@ ccBactrian :: Cycle Double ccBactrian = cycleFromList [slideBactrian 0.5 1.0 (PName "Bactrian") (pWeight 1) Tune] -normalBactrianBench :: IOGenM StdGen -> IO ()+normalBactrianBench :: StdGen -> IO () normalBactrianBench g = do   let s =         Settings@@ -99,7 +99,7 @@   a <- mhg s noPrior lh ccBactrian mon 0 g   void $ mcmc s a -normalMC3 :: IOGenM StdGen -> Int -> IO ()+normalMC3 :: StdGen -> Int -> IO () normalMC3 g n = do   let mcmcS =         Settings

@@ -86,7 +86,7 @@ mon :: Monitor I mon = Monitor monStd [] [] -poissonBench :: IOGenM StdGen -> IO ()+poissonBench :: StdGen -> IO () poissonBench g = do   let s =         Settings@@ -123,7 +123,7 @@ hmcProposal :: Cycle I hmcProposal = cycleFromList [hamiltonian hparams htconf hstruct target (PName "Hamiltonian") (pWeight 1)] -poissonHamiltonianBench :: IOGenM StdGen -> IO ()+poissonHamiltonianBench :: StdGen -> IO () poissonHamiltonianBench g = do   let s =         Settings

@@ -1,6 +1,6 @@ cabal-version:      3.0 name:               mcmc-version:            0.7.0.0+version:            0.7.0.1 synopsis:           Sample from a posterior using Markov chain Monte Carlo description:   Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme>

@@ -30,7 +30,7 @@    -- | Check if the current state is invalid. At the moment this should just   -- check if the posterior probability is zero or NaN.-  aIsInValidState :: a -> Bool+  aIsInvalidState :: a -> Bool    -- | Sample the next state.   aIterate :: IterationMode -> ParallelizationMode -> a -> IO a@@ -48,7 +48,7 @@   -- the trace, if required.   aCleanAfterBurnIn :: TraceLength -> a -> IO a -  -- | Summarize the cycle.+  -- | Summarize the proposals in the cycle.   aSummarizeCycle :: IterationMode -> a -> BL.ByteString    -- | Open required monitor files and setup corresponding file handles.

@@ -179,7 +179,7 @@ instance ToJSON a => Algorithm (MC3 a) where   aName = const "Metropolis-coupled Markov chain Monte Carlo (MC3)"   aIteration = mc3Iteration-  aIsInValidState = mc3IsInValidState+  aIsInvalidState = mc3IsInvalidState   aIterate = mc3Iterate   aAutoTune = mc3AutoTune   aResetAcceptance = mc3ResetAcceptance@@ -265,10 +265,6 @@     l = link c     t = trace c --- TODO (medium): Splitmix. Initialization of the MC3 algorithm is an IO action--- because the generators have to be split. And also because of the mutable--- trace.- -- | Initialize an MC3 algorithm with a given number of chains. -- -- Call 'error' if:@@ -284,19 +280,21 @@   Cycle a ->   Monitor a ->   InitialState a ->-  IOGenM StdGen ->+  StdGen ->   IO (MC3 a) mc3 sMc3 s pr lh cc mn i0 g   | n < 2 = error "mc3: The number of chains must be two or larger."   | sp < 1 = error "mc3: The swap period must be strictly positive."   | sn < 1 || sn > n - 1 = error "mc3: The number of swaps must be in [1, NChains - 1]."   | otherwise = do-      -- Split random number generators.-      rs <- replicateM n $ splitGenM g-      gs <- V.fromList <$> mapM newIOGenM rs-      cs <- V.mapM (mhg s pr lh cc mn i0) gs+      -- Split random number generator.+      let gs = take (n + 1) $ unfoldr (Just . split) g+      -- Prepare MHG chains.+      cs <- V.mapM (mhg s pr lh cc mn i0) (V.fromList $ tail gs)       hcs <- V.izipWithM (initMHG pr lh) (V.convert bs) cs-      return $ MC3 sMc3 hcs bs 0 (emptyA [0 .. n - 2]) g+      -- Do not reuse the initial generator.+      gm <- newIOGenM $ head gs+      return $ MC3 sMc3 hcs bs 0 (emptyA [0 .. n - 2]) gm   where     n = fromNChains $ mc3NChains sMc3     sp = fromSwapPeriod $ mc3SwapPeriod sMc3@@ -416,16 +414,13 @@   where     g = mc3Generator a -mc3IsInValidState :: ToJSON a => MC3 a -> Bool-mc3IsInValidState a = V.any aIsInValidState mhgs+mc3IsInvalidState :: ToJSON a => MC3 a -> Bool+mc3IsInvalidState a = V.any aIsInvalidState mhgs   where     mhgs = mc3MHGChains a --- TODO (medium): Splitmix. 'mc3Iterate' is actually not parallel, but--- concurrent because of the IO constraint. Use pure parallel code when we have--- a pure generator.------ However, we have to take care of the mutable traces.+-- NOTE: 'mc3Iterate' is actually not parallel, but concurrent because of the IO+-- constraint of the mutable trace. mc3Iterate ::   ToJSON a =>   IterationMode ->

@@ -61,7 +61,7 @@ instance ToJSON a => Algorithm (MHG a) where   aName = const "Metropolis-Hastings-Green (MHG)"   aIteration = iteration . fromMHG-  aIsInValidState = mhgIsInValidState+  aIsInvalidState = mhgIsInvalidState   aIterate = mhgIterate   aAutoTune = mhgAutoTune   aResetAcceptance = mhgResetAcceptance@@ -73,6 +73,9 @@   aCloseMonitors = mhgCloseMonitors   aSave = mhgSave +-- Calculate required length of trace. The length may be larger during burn in,+-- because the tuners of some proposals (e.g., HMC, NUTS) require the states of+-- the last tuning interval. getTraceLength ::   Maybe BurnInSettings ->   TraceLength ->@@ -101,13 +104,14 @@   Cycle a ->   Monitor a ->   InitialState a ->-  IOGenM StdGen ->+  StdGen ->   IO (MHG a) mhg s pr lh cc mn i0 g = do   -- The trace is a mutable vector and the mutable state needs to be handled by   -- a monad.   tr <- replicateT tl l0-  return $ MHG $ Chain Nothing l0 0 tr ac g 0 pr lh cc mn+  gm <- newIOGenM g+  return $ MHG $ Chain Nothing l0 0 tr ac gm 0 pr lh cc mn   where     l0 = Link i0 (pr i0) (lh i0)     ac = emptyA $ ccProposals cc@@ -280,8 +284,8 @@ -- -- At the moment this just checks whether the prior, likelihood, or posterior -- are NaN or infinite.-mhgIsInValidState :: MHG a -> Bool-mhgIsInValidState a = checkSoft p || check l || check (p * l)+mhgIsInvalidState :: MHG a -> Bool+mhgIsInvalidState a = checkSoft p || check l || check (p * l)   where     x = link $ fromMHG a     p = prior x

@@ -27,6 +27,9 @@ import qualified Data.Vector as VB import Mcmc.Chain.Link +-- NOTE: We directly refer to the 'PrimSate' 'RealWorld' because, otherwise, we+-- have two type variables in 'Chain m a'.+ -- | A 'Trace' is a mutable circular stack that passes through a list of states -- with associated priors and likelihoods called 'Link's. newtype Trace a = Trace {fromTrace :: C.MStack VB.Vector RealWorld (Link a)}

@@ -28,8 +28,8 @@ -- cleaner, in my opinion. However, I would like to move away from MWC so I -- leave the custom implementation for now. ----- NOTE: Moving away from MWC is not anymore true because I can use SplitMix--- with it.+-- NOTE: I do not need to move away from MWC anymore true because I can use+-- SplitMix with it.  -- | Shuffle a vector. shuffle :: StatefulGen g m => [a] -> g -> m [a]

@@ -223,7 +223,7 @@   Cycle a ->   Monitor a ->   a ->-  IOGenM StdGen ->+  StdGen ->   -- For each point a vector of likelihoods stored in log domain.   ML [VU.Vector Likelihood] mlRun k xs em vb prf lhf cc mn i0 g = do@@ -270,12 +270,11 @@   Cycle a ->   Monitor a ->   a ->-  IOGenM StdGen ->+  StdGen ->   ML MarginalLikelihood-tiWrapper s prf lhf cc mn i0 g0 = do+tiWrapper s prf lhf cc mn i0 g = do   logInfoB "Path integral (thermodynamic integration)."-  r1 <- splitGenM g0-  g1 <- newIOGenM r1+  let (g0, g1) = split g    -- Parallel execution of both path integrals.   r <- ask@@ -354,7 +353,7 @@   Cycle a ->   Monitor a ->   a ->-  IOGenM StdGen ->+  StdGen ->   ML MarginalLikelihood sssWrapper s prf lhf cc mn i0 g = do   logInfoB "Stepping stone sampling."@@ -378,7 +377,7 @@   Cycle a ->   Monitor a ->   InitialState a ->-  IOGenM StdGen ->+  StdGen ->   IO MarginalLikelihood marginalLikelihood s prf lhf cc mn i0 g = do   -- Initialize.

@@ -31,7 +31,6 @@ import Mcmc.Logger import Mcmc.Proposal (TuningType (LastTuningStep, NormalTuningStep)) import Mcmc.Settings-import System.Exit import System.IO import Prelude hiding (cycle) @@ -109,7 +108,7 @@   logInfoB "Initial state."   logInfoB $ aStdMonitorHeader a   mcmcExecuteMonitors a-  when (aIsInValidState a) (logWarnB "The initial state is invalid!")+  when (aIsInvalidState a) (logWarnB "The initial state is invalid!")   a' <- mcmcBurnIn a   logInfoS $ "Clean chain after burn in."   let tl = sTraceLength s

@@ -199,8 +199,8 @@     -- the ratio of the backward to forward kernels (the 'KernelRatio' or the     -- probability masses or probability densities) and the 'Jacobian'.     ---    -- NOTE: Actually the 'Jacobian' should be part of the 'KernelRatio'. However,-    -- it is more declarative to have them separate. Like so, we are constantly+    -- Note: The 'Jacobian' should be part of the 'KernelRatio'. However, it is+    -- more declarative to have them separate. Like so, we are constantly     -- reminded: Is the Jacobian modifier different from 1.0?     Propose !a !KernelRatio !Jacobian   deriving (Show, Eq)

@@ -66,7 +66,7 @@   describe "save and load" $     it "doesn't change the MCMC chain" $       do-        gen <- newIOGenM $ mkStdGen 0+        let gen = mkStdGen 0         a <- mhg settings noPrior lh proposals mon 0 gen         c <- fromMHG <$> mcmc settings a         savedChain <- toSavedChain c@@ -95,12 +95,11 @@   describe "mhContinue" $     it "mcmc 50 + mcmcContinue 50 == mcmc 100" $       do-        gen1 <- newIOGenM $ mkStdGen 0-        a1 <- mhg settings noPrior lh proposals mon 0 gen1+        let gen = mkStdGen 0+        a1 <- mhg settings noPrior lh proposals mon 0 gen         r1 <- fromMHG <$> mcmc settings a1-        gen2 <- newIOGenM $ mkStdGen 0         let settings' = settings {sIterations = Iterations 50}-        a2 <- mhg settings' noPrior lh proposals mon 0 gen2+        a2 <- mhg settings' noPrior lh proposals mon 0 gen         r2' <- mcmc settings' a2         r2 <- fromMHG <$> mcmcContinue (Iterations 50) settings' r2'         link r1 `shouldBe` link r2