mcmc 0.8.1.0 → 0.8.2.0

raw patch · 2 files changed

+87/−27 lines, 2 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• mcmc.cabal +1/−1
• src/Mcmc/MarginalLikelihood.hs +86/−26

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

@@ -21,6 +21,7 @@   ) where +import Control.Concurrent (getNumCapabilities) import Control.Concurrent.Async hiding (link) import Control.Monad import Control.Monad.IO.Class@@ -87,13 +88,15 @@   { mlAnalysisName :: AnalysisName,     mlAlgorithm :: MLAlgorithm,     mlNPoints :: NPoints,-    -- | Initial burn in at the starting point of the path.+    -- | Initial burn in at the starting point of the path (or each segment if+    -- running in parallel).     mlInitialBurnIn :: BurnInSettings,     -- | Repetitive burn in at each point on the path.     mlPointBurnIn :: BurnInSettings,     -- | The number of iterations performed at each point.     mlIterations :: Iterations,     mlExecutionMode :: ExecutionMode,+    mlParallelizationMode :: ParallelizationMode,     mlLogMode :: LogMode,     mlVerbosity :: Verbosity   }@@ -139,23 +142,27 @@  sampleAtPoint ::   ToJSON a =>+  Bool ->   Point ->   Settings ->   LikelihoodFunction a ->   MHG a ->   ML (MHG a)-sampleAtPoint x ss lhf a = do+sampleAtPoint isInitialBurnIn x ss lhf a = do   a'' <- liftIO $ mcmc ss' a'   let ch'' = fromMHG a''       ac = acceptances ch''       mAr = sequence $ acceptanceRates ac   logDebugB "sampleAtPoint: Summarize cycle."   logDebugB $ summarizeCycle AllProposals ac $ cycle ch''-  case mAr of-    Nothing -> logWarnB "Some acceptance rates are unavailable."-    Just ar -> do-      unless (M.null $ M.filter (<= 0.1) ar) $ logWarnB "Some acceptance rates are below 0.1."-      unless (M.null $ M.filter (>= 0.9) ar) $ logWarnB "Some acceptance rates are above 0.9."+  unless+    isInitialBurnIn+    ( case mAr of+        Nothing -> logWarnB "Some acceptance rates are unavailable."+        Just ar -> do+          unless (M.null $ M.filter (<= 0.1) ar) $ logWarnB "Some acceptance rates are below 0.1."+          unless (M.null $ M.filter (>= 0.9) ar) $ logWarnB "Some acceptance rates are above 0.9."+    )   return a''   where     -- For debugging set a proper analysis name.@@ -181,19 +188,18 @@  traversePoints ::   ToJSON a =>-  -- Current point.-  Int ->   NPoints ->-  [Point] ->+  [(Int, Point)] ->   Settings ->   LikelihoodFunction a ->   MHG a ->   -- For each point a vector of obtained likelihoods stored in the log domain.   ML [VU.Vector Likelihood]-traversePoints _ _ [] _ _ _ = return []-traversePoints i k (b : bs) ss lhf a = do-  logInfoS $ "Point " <> show i <> " of " <> show k' <> ": " <> show b <> "."-  a' <- sampleAtPoint b ss lhf a+traversePoints _ [] _ _ _ = return []+traversePoints k ((idb, b) : bs) ss lhf a = do+  let msg = printf "Point %4d of %4d: %.12f." idb k' b+  logInfoS msg+  a' <- sampleAtPoint False b ss lhf a   -- Get the links samples at this point.   ls <- liftIO $ takeT n $ trace $ fromMHG a'   -- Extract the likelihoods.@@ -205,16 +211,67 @@   -- resulting in a severe memory leak.   let !lhs = VU.convert $ VB.map (lhf . state) ls   -- Sample the other points.-  lhss <- traversePoints (i + 1) k bs ss lhf a'+  lhss <- traversePoints k bs ss lhf a'   return $ lhs : lhss   where     n = fromIterations $ sIterations ss     (NPoints k') = k +nChunks :: Int -> [a] -> [[a]]+nChunks k xs = chop (chunks k l) xs+  where+    l = length xs++chunks :: Int -> Int -> [Int]+chunks c n = filter (> 0) ns+  where+    n' = n `div` c+    r = n `mod` c+    ns = replicate r (n' + 1) ++ replicate (c - r) n'++chop :: [Int] -> [a] -> [[a]]+chop [] [] = []+chop (n : ns) xs+  | n > 0 = take n xs : chop ns (drop n xs)+  | otherwise = error "chop: n negative or zero"+chop _ _ = error "chop: not all list elements handled"++mlRunPar ::+  ToJSON a =>+  ParallelizationMode ->+  NPoints ->+  [(Int, Point)] ->+  ExecutionMode ->+  Verbosity ->+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  a ->+  StdGen ->+  ML [VU.Vector Likelihood]+mlRunPar pm k xs em vb prf lhf cc mn i0 g = do+  nThreads <- case pm of+    Sequential -> do+      logInfoB "Sequential execution."+      pure 1+    Parallel -> do+      n <- liftIO getNumCapabilities+      logInfoS $ "Parallel execution with " <> show n <> " cores."+      pure n+  let xsChunks = nChunks nThreads xs+  r <- ask+  xss <-+    liftIO $+      mapConcurrently+        (\thesePoints -> runReaderT (mlRun k thesePoints em vb prf lhf cc mn i0 g) r)+        xsChunks+  pure $ concat xss+ mlRun ::   ToJSON a =>   NPoints ->-  [Point] ->+  [(Int, Point)] ->   ExecutionMode ->   Verbosity ->   PriorFunction a ->@@ -241,12 +298,13 @@       ssP = Settings nm biP is trLen em Sequential NoSave LogFileOnly vb'   logDebugB "mlRun: Initialize MHG algorithm."   a0 <- liftIO $ mhg ssI prf lhf cc mn i0 g-  logDebugS $ "mlRun: Perform initial burn in at first point " <> show x0 <> "."-  a1 <- sampleAtPoint x0 ssI lhf a0+  let msg = printf "Initial burn in at point %.12f with ID %4d." x0 id0+  logInfoS msg+  a1 <- sampleAtPoint True x0 ssI lhf a0   logDebugB "mlRun: Traverse points."-  traversePoints 1 k xs ssP lhf a1+  traversePoints k xs ssP lhf a1   where-    x0 = head xs+    (id0, x0) = head xs  -- Use lists since the number of points is expected to be low. integrateSimpsonTriangle ::@@ -280,8 +338,8 @@   (lhssForward, lhssBackward) <-     lift $       concurrently-        (runReaderT (mlRun k bsForward em vb prf lhf cc mn i0 g0) r)-        (runReaderT (mlRun k bsBackward em vb prf lhf cc mn i0 g1) r)+        (runReaderT (mlRunPar pm k (zip [1 ..] bsForward) em vb prf lhf cc mn i0 g0) r)+        (runReaderT (mlRunPar pm k (zip [1 ..] bsBackward) em vb prf lhf cc mn i0 g1) r)   logInfoEndTime    logDebugB "tiWrapper: Calculate mean log likelihoods."@@ -300,6 +358,7 @@     bsForward = getPoints k     bsBackward = reverse bsForward     em = mlExecutionMode s+    pm = mlParallelizationMode s     vb = mlVerbosity s  -- Helper function to exponentiate log domain values with a double value.@@ -356,7 +415,7 @@   ML MarginalLikelihood sssWrapper s prf lhf cc mn i0 g = do   logInfoB "Stepping stone sampling."-  logLhss <- mlRun k bsForward' em vb prf lhf cc mn i0 g+  logLhss <- mlRunPar pm k (zip [1 ..] bsForward') em vb prf lhf cc mn i0 g   logInfoB "The last point does not need to be sampled with stepping stone sampling."   logDebugB "sssWrapper: Calculate marginal likelihood."   return $ sssCalculateMarginalLikelihood bsForward logLhss@@ -365,6 +424,7 @@     bsForward = getPoints k     bsForward' = init bsForward     em = mlExecutionMode s+    pm = mlParallelizationMode s     vb = mlVerbosity s  -- | Estimate the marginal likelihood.@@ -382,9 +442,9 @@   -- Initialize.   e <- initializeEnvironment s -  when (mlVerbosity s == Debug) $ do-    let n = fromAnalysisName $ mlAnalysisName s-    createDirectoryIfMissing True n+  -- Create marginal likelihood analysis directory.+  let n = fromAnalysisName $ mlAnalysisName s+  createDirectoryIfMissing True n    -- Run.   runReaderT