mcmc 0.8.2.0 → 0.8.3.0

raw patch · 21 files changed

+152/−144 lines, 21 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• ChangeLog.md +6/−1
• bench/Bench.hs +1/−1
• mcmc.cabal +1/−1
• src/Mcmc/Acceptance.hs +6/−6
• src/Mcmc/Algorithm/MC3.hs +17/−14
• src/Mcmc/Algorithm/MHG.hs +5/−5
• src/Mcmc/Chain/Link.hs +3/−3
• src/Mcmc/Cycle.hs +1/−1
• src/Mcmc/Environment.hs +3/−3
• src/Mcmc/Internal/Shuffle.hs +2/−2
• src/Mcmc/Internal/SpecFunctions.hs +19/−19
• src/Mcmc/Likelihood.hs +1/−1
• src/Mcmc/Logger.hs +1/−1
• src/Mcmc/MarginalLikelihood.hs +7/−7
• src/Mcmc/Mcmc.hs +53/−53
• src/Mcmc/Monitor/ParameterBatch.hs +4/−4
• src/Mcmc/Monitor/Time.hs +1/−1
• src/Mcmc/Prior.hs +14/−14
• src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs +2/−2
• src/Mcmc/Proposal/Hamiltonian/Internal.hs +3/−3
• src/Mcmc/Proposal/Hamiltonian/Nuts.hs +2/−2

@@ -5,6 +5,11 @@ ## Unreleased changes  +## 0.8.3.0++-   Fix auto tuning with MC3 algorithm.++ ## 0.8.1.0  -   Automatic intermediate tuning for HMC and NUTS.@@ -161,7 +166,7 @@ ## 0.2.4  -   **Change order of arguments for proposals**.--   &rsquo;slideStem&rsquo; was renamed to &rsquo;slideBranch&rsquo;.+-   'slideStem' was renamed to 'slideBranch'. -   Change ProposalSimple from newtype to type. -   Contravariant instances of parameter and batch monitors. Use `(>$<)` instead     of `(@.)` and `(@#)`.

@@ -21,7 +21,7 @@ import Poisson import System.Random.Stateful -gammaBenchG :: RealFloat a => (a -> a) -> [a] -> a+gammaBenchG :: (RealFloat a) => (a -> a) -> [a] -> a gammaBenchG f = foldl' (\acc x -> acc + f x) 0 {-# SPECIALIZE gammaBenchG :: (Double -> Double) -> [Double] -> Double #-} 

@@ -1,6 +1,6 @@ cabal-version:      3.0 name:               mcmc-version:            0.8.2.0+version:            0.8.3.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>

@@ -78,7 +78,7 @@ newtype Acceptances k = Acceptances {fromAcceptances :: M.Map k Acceptance}   deriving (Eq, Show) -instance ToJSONKey k => ToJSON (Acceptances k) where+instance (ToJSONKey k) => ToJSON (Acceptances k) where   toJSON (Acceptances m) = toJSON m   toEncoding (Acceptances m) = toEncoding m @@ -88,17 +88,17 @@ -- | In the beginning there was the Word. -- -- Initialize an empty storage of accepted/rejected values.-emptyA :: Ord k => [k] -> Acceptances k+emptyA :: (Ord k) => [k] -> Acceptances k emptyA ks = Acceptances $ M.fromList [(k, A noCounts Nothing) | k <- ks]   where     noCounts = AcceptanceCounts 0 0  -- | For key @k@, add an accept.-pushAccept :: Ord k => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k+pushAccept :: (Ord k) => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k pushAccept mr k = Acceptances . M.adjust (addAccept mr) k . fromAcceptances  -- | For key @k@, add a reject.-pushReject :: Ord k => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k+pushReject :: (Ord k) => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k pushReject mr k = Acceptances . M.adjust (addReject mr) k . fromAcceptances  -- | Reset acceptance specification.@@ -109,7 +109,7 @@     ResetExpectedRatesOnly  -- | Reset acceptance counts.-resetA :: Ord k => ResetAcceptance -> Acceptances k -> Acceptances k+resetA :: (Ord k) => ResetAcceptance -> Acceptances k -> Acceptances k resetA ResetEverything = emptyA . M.keys . fromAcceptances resetA ResetExpectedRatesOnly = Acceptances . M.map f . fromAcceptances   where@@ -133,7 +133,7 @@ -- Return 'Nothing' if no proposals have been accepted or rejected (division by -- zero). acceptanceRate ::-  Ord k =>+  (Ord k) =>   k ->   Acceptances k ->   -- | (nAccepts, nRejects, actualRate, expectedRate)

@@ -176,7 +176,7 @@     mc3Generator :: IOGenM StdGen   } -instance ToJSON a => Algorithm (MC3 a) where+instance (ToJSON a) => Algorithm (MC3 a) where   aName = const "Metropolis-coupled Markov chain Monte Carlo (MC3)"   aIteration = mc3Iteration   aIsInvalidState = mc3IsInvalidState@@ -311,7 +311,7 @@  -- | Save an MC3 algorithm. mc3Save ::-  ToJSON a =>+  (ToJSON a) =>   AnalysisName ->   MC3 a ->   IO ()@@ -325,7 +325,7 @@ -- -- See 'Mcmc.Mcmc.mcmcContinue'. mc3Load ::-  FromJSON a =>+  (FromJSON a) =>   PriorFunction a ->   LikelihoodFunction a ->   Cycle a ->@@ -419,7 +419,7 @@   where     g = mc3Generator a -mc3IsInvalidState :: ToJSON a => MC3 a -> Bool+mc3IsInvalidState :: (ToJSON a) => MC3 a -> Bool mc3IsInvalidState a = V.any aIsInvalidState mhgs   where     mhgs = mc3MHGChains a@@ -427,7 +427,7 @@ -- NOTE: 'mc3Iterate' is actually not parallel, but concurrent because of the IO -- constraint of the mutable trace. mc3Iterate ::-  ToJSON a =>+  (ToJSON a) =>   IterationMode ->   ParallelizationMode ->   MC3 a ->@@ -477,7 +477,7 @@     rNew = (brOld / blOld) ** xi     brNew = blNew * rNew -mc3AutoTune :: ToJSON a => TuningType -> Int -> MC3 a -> IO (MC3 a)+mc3AutoTune :: (ToJSON a) => TuningType -> Int -> MC3 a -> IO (MC3 a) mc3AutoTune b l a = do   -- 1. Auto tune all chains.   mhgs' <- V.mapM (aAutoTune b l) $ mc3MHGChains a@@ -510,9 +510,12 @@             (setReciprocalTemperature coldPrF coldLhF)             (V.convert $ U.tail bs')             (V.tail mhgs')-  return $ a {mc3MHGChains = mhgs'', mc3ReciprocalTemperatures = bs'}+  return $+    if b == NormalTuningFastProposalsOnly || b == NormalTuningAllProposals+      then a {mc3MHGChains = mhgs'', mc3ReciprocalTemperatures = bs'}+      else a {mc3MHGChains = mhgs''} -mc3ResetAcceptance :: ToJSON a => ResetAcceptance -> MC3 a -> MC3 a+mc3ResetAcceptance :: (ToJSON a) => ResetAcceptance -> MC3 a -> MC3 a mc3ResetAcceptance x a = a'   where     -- 1. Reset acceptance of all chains.@@ -522,7 +525,7 @@     --     a' = a {mc3MHGChains = mhgs', mc3SwapAcceptances = ac'} -mc3CleanAfterBurnIn :: ToJSON a => TraceLength -> MC3 a -> IO (MC3 a)+mc3CleanAfterBurnIn :: (ToJSON a) => TraceLength -> MC3 a -> IO (MC3 a) mc3CleanAfterBurnIn tl a = do   cs' <- V.mapM (aCleanAfterBurnIn tl) cs   pure $ a {mc3MHGChains = cs'}@@ -536,7 +539,7 @@ -- - The combined acceptance rate of proposals within the hot chains. -- -- - The temperatures of the chains and the acceptance rates of the state swaps.-mc3SummarizeCycle :: ToJSON a => IterationMode -> MC3 a -> BL.ByteString+mc3SummarizeCycle :: (ToJSON a) => IterationMode -> MC3 a -> BL.ByteString mc3SummarizeCycle m a =   BL.intercalate "\n" $     [ "MC3: Cycle of cold chain.",@@ -585,7 +588,7 @@     proposalHLine = BL.replicate (BL.length proposalHeader) '-'  -- No extra monitors are opened.-mc3OpenMonitors :: ToJSON a => AnalysisName -> ExecutionMode -> MC3 a -> IO (MC3 a)+mc3OpenMonitors :: (ToJSON a) => AnalysisName -> ExecutionMode -> MC3 a -> IO (MC3 a) mc3OpenMonitors nm em a = do   mhgs' <- V.imapM mhgOpenMonitors (mc3MHGChains a)   return $ a {mc3MHGChains = mhgs'}@@ -598,7 +601,7 @@         suf = printf "%02d" i  mc3ExecuteMonitors ::-  ToJSON a =>+  (ToJSON a) =>   Verbosity ->   -- Starting time.   UTCTime ->@@ -613,10 +616,10 @@     -- All other chains are to be quiet.     f _ = aExecuteMonitors Quiet t0 iTotal -mc3StdMonitorHeader :: ToJSON a => MC3 a -> BL.ByteString+mc3StdMonitorHeader :: (ToJSON a) => MC3 a -> BL.ByteString mc3StdMonitorHeader = aStdMonitorHeader . V.head . mc3MHGChains -mc3CloseMonitors :: ToJSON a => MC3 a -> IO (MC3 a)+mc3CloseMonitors :: (ToJSON a) => MC3 a -> IO (MC3 a) mc3CloseMonitors a = do   mhgs' <- V.mapM aCloseMonitors $ mc3MHGChains a   return $ a {mc3MHGChains = mhgs'}

@@ -57,7 +57,7 @@ -- | The MHG algorithm. newtype MHG a = MHG {fromMHG :: Chain a} -instance ToJSON a => Algorithm (MHG a) where+instance (ToJSON a) => Algorithm (MHG a) where   aName = const "Metropolis-Hastings-Green (MHG)"   aIteration = iteration . fromMHG   aIsInvalidState = mhgIsInvalidState@@ -121,7 +121,7 @@  -- | Save an MHG algorithm. mhgSave ::-  ToJSON a =>+  (ToJSON a) =>   AnalysisName ->   MHG a ->   IO ()@@ -135,7 +135,7 @@ -- -- See 'Mcmc.Mcmc.mcmcContinue'. mhgLoad ::-  FromJSON a =>+  (FromJSON a) =>   PriorFunction a ->   LikelihoodFunction a ->   Cycle a ->@@ -151,7 +151,7 @@ -- Useful when restarting a run with changed prior function, likelihood function -- or proposals. Use with care! mhgLoadUnsafe ::-  FromJSON a =>+  (FromJSON a) =>   PriorFunction a ->   LikelihoodFunction a ->   Cycle a ->@@ -162,7 +162,7 @@  -- Nice type :-). mhgLoadWith ::-  FromJSON a =>+  (FromJSON a) =>   (PriorFunction a -> LikelihoodFunction a -> Cycle a -> Monitor a -> SavedChain a -> IO (Chain a)) ->   PriorFunction a ->   LikelihoodFunction a ->

@@ -34,16 +34,16 @@   }   deriving (Eq, Ord, Show, Read) -instance ToJSON a => ToJSON (Link a) where+instance (ToJSON a) => ToJSON (Link a) where   toJSON (Link x (Exp p) (Exp l)) = object ["s" .= x, "p" .= p, "l" .= l]   toEncoding (Link x (Exp p) (Exp l)) = pairs ("s" .= x <> "p" .= p <> "l" .= l) -link :: FromJSON a => Object -> Parser (Link a)+link :: (FromJSON a) => Object -> Parser (Link a) link v = do   s <- v .: "s"   p <- v .: "p"   l <- v .: "l"   return $ Link s (Exp p) (Exp l) -instance FromJSON a => FromJSON (Link a) where+instance (FromJSON a) => FromJSON (Link a) where   parseJSON = withObject "Link" link

@@ -126,7 +126,7 @@   deriving (Eq)  -- | Replicate 'Proposal's according to their weights and possibly shuffle them.-prepareProposals :: StatefulGen g m => IterationMode -> Cycle a -> g -> m [Proposal a]+prepareProposals :: (StatefulGen g m) => IterationMode -> Cycle a -> g -> m [Proposal a] prepareProposals m (Cycle xs o _ _) g =   if null ps     then

@@ -39,7 +39,7 @@   }   deriving (Eq) -instance HasExecutionMode s => HasExecutionMode (Environment s) where+instance (HasExecutionMode s) => HasExecutionMode (Environment s) where   getExecutionMode = getExecutionMode . settings  instance HasLock (Environment s) where@@ -51,10 +51,10 @@ instance HasStartingTime (Environment s) where   getStartingTime = startingTime -instance HasLogMode s => HasLogMode (Environment s) where+instance (HasLogMode s) => HasLogMode (Environment s) where   getLogMode = getLogMode . settings -instance HasVerbosity s => HasVerbosity (Environment s) where+instance (HasVerbosity s) => HasVerbosity (Environment s) where   getVerbosity = getVerbosity . settings  -- | Initialize the environment.

@@ -26,12 +26,12 @@ -- vectors.  -- | Shuffle a vector.-shuffle :: StatefulGen g m => [a] -> g -> m [a]+shuffle :: (StatefulGen g m) => [a] -> g -> m [a] shuffle xs = grabble xs (length xs)  -- @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@ -- elements from @xs@, without replacement, and that @m@ times.-grabble :: StatefulGen g m => [a] -> Int -> g -> m [a]+grabble :: (StatefulGen g m) => [a] -> Int -> g -> m [a] grabble xs m gen = do   swaps <- forM [0 .. min (l - 1) m] $ \i -> do     j <- uniformRM (i, l) gen

@@ -25,10 +25,10 @@ import Numeric.SpecFunctions import Unsafe.Coerce -mSqrtEpsG :: RealFloat a => a+mSqrtEpsG :: (RealFloat a) => a mSqrtEpsG = 1.4901161193847656e-8 -mEulerMascheroniG :: RealFloat a => a+mEulerMascheroniG :: (RealFloat a) => a mEulerMascheroniG = 0.5772156649015328606065121  -- | Generalized version of the log gamma distribution. See@@ -40,7 +40,7 @@ {-# SPECIALIZE logGammaG :: Double -> Double #-}  -- See 'Numeric.SpecFunctions.logGamma'.-logGammaNonDouble :: RealFloat a => a -> a+logGammaNonDouble :: (RealFloat a) => a -> a logGammaNonDouble z   | z <= 0 = 1 / 0   | z < mSqrtEpsG = log (1 / z - mEulerMascheroniG)@@ -51,7 +51,7 @@   | z < 15 = lgammaSmallG z   | otherwise = lanczosApproxG z -lgamma1_15G :: RealFloat a => a -> a -> a+lgamma1_15G :: (RealFloat a) => a -> a -> a lgamma1_15G zm1 zm2 =   r * y     + r@@ -62,7 +62,7 @@     r = zm1 * zm2     y = 0.52815341949462890625 -tableLogGamma_1_15PG :: RealFloat a => VB.Vector a+tableLogGamma_1_15PG :: (RealFloat a) => VB.Vector a tableLogGamma_1_15PG =   VB.fromList     [ 0.490622454069039543534e-1,@@ -75,7 +75,7 @@     ] {-# NOINLINE tableLogGamma_1_15PG #-} -tableLogGamma_1_15QG :: RealFloat a => VB.Vector a+tableLogGamma_1_15QG :: (RealFloat a) => VB.Vector a tableLogGamma_1_15QG =   VB.fromList     [ 1,@@ -88,7 +88,7 @@     ] {-# NOINLINE tableLogGamma_1_15QG #-} -lgamma15_2G :: RealFloat a => a -> a -> a+lgamma15_2G :: (RealFloat a) => a -> a -> a lgamma15_2G zm1 zm2 =   r * y     + r@@ -99,7 +99,7 @@     r = zm1 * zm2     y = 0.452017307281494140625 -tableLogGamma_15_2PG :: RealFloat a => VB.Vector a+tableLogGamma_15_2PG :: (RealFloat a) => VB.Vector a tableLogGamma_15_2PG =   VB.fromList     [ -0.292329721830270012337e-1,@@ -111,7 +111,7 @@     ] {-# NOINLINE tableLogGamma_15_2PG #-} -tableLogGamma_15_2QG :: RealFloat a => VB.Vector a+tableLogGamma_15_2QG :: (RealFloat a) => VB.Vector a tableLogGamma_15_2QG =   VB.fromList     [ 1,@@ -124,7 +124,7 @@     ] {-# NOINLINE tableLogGamma_15_2QG #-} -lgammaSmallG :: RealFloat a => a -> a+lgammaSmallG :: (RealFloat a) => a -> a lgammaSmallG = go 0   where     go acc z@@ -133,7 +133,7 @@       where         zm1 = z - 1 -lgamma2_3G :: RealFloat a => a -> a+lgamma2_3G :: (RealFloat a) => a -> a lgamma2_3G z =   r * y     + r@@ -145,7 +145,7 @@     zm2 = z - 2     y = 0.158963680267333984375e0 -tableLogGamma_2_3PG :: RealFloat a => VB.Vector a+tableLogGamma_2_3PG :: (RealFloat a) => VB.Vector a tableLogGamma_2_3PG =   VB.fromList     [ -0.180355685678449379109e-1,@@ -158,7 +158,7 @@     ] {-# NOINLINE tableLogGamma_2_3PG #-} -tableLogGamma_2_3QG :: RealFloat a => VB.Vector a+tableLogGamma_2_3QG :: (RealFloat a) => VB.Vector a tableLogGamma_2_3QG =   VB.fromList     [ 1,@@ -172,14 +172,14 @@     ] {-# NOINLINE tableLogGamma_2_3QG #-} -lanczosApproxG :: RealFloat a => a -> a+lanczosApproxG :: (RealFloat a) => a -> a lanczosApproxG z =   (log (z + g - 0.5) - 1) * (z - 0.5)     + log (evalRatioG tableLanczosG z)   where     g = 6.024680040776729583740234375 -tableLanczosG :: RealFloat a => VB.Vector (a, a)+tableLanczosG :: (RealFloat a) => VB.Vector (a, a) tableLanczosG =   VB.fromList     [ (56906521.91347156388090791033559122686859, 0),@@ -200,7 +200,7 @@  data LG a = LG !a !a -evalRatioG :: RealFloat a => VB.Vector (a, a) -> a -> a+evalRatioG :: (RealFloat a) => VB.Vector (a, a) -> a -> a evalRatioG coef x   | x > 1 = fini $ VB.foldl' stepL (LG 0 0) coef   | otherwise = fini $ VB.foldr' stepR (LG 0 0) coef@@ -228,13 +228,13 @@     stirling = (x - 0.5) * log x - x + mLnSqrt2Pi     x = fromIntegral n + 1     rx = recip x-{-# SPECIALIZE logFactorialNonDouble :: RealFloat a => Int -> a #-}+{-# SPECIALIZE logFactorialNonDouble :: (RealFloat a) => Int -> a #-} -mLnSqrt2Pi :: RealFloat a => a+mLnSqrt2Pi :: (RealFloat a) => a mLnSqrt2Pi = 0.9189385332046727417803297364056176398613974736377834128171 {-# INLINE mLnSqrt2Pi #-} -factorialTable :: RealFloat a => VB.Vector a+factorialTable :: (RealFloat a) => VB.Vector a {-# NOINLINE factorialTable #-} factorialTable =   VB.fromListN

@@ -29,6 +29,6 @@ type LikelihoodFunctionG a b = a -> Log b  -- | Flat likelihood function. Useful for testing and debugging.-noLikelihood :: RealFloat b => LikelihoodFunctionG a b+noLikelihood :: (RealFloat b) => LikelihoodFunctionG a b noLikelihood = const 1.0 {-# SPECIALIZE noLikelihood :: LikelihoodFunction Double #-}

@@ -88,7 +88,7 @@ msgPrepare pref msg = BL.intercalate "\n" $ map (BL.append pref) $ BL.lines msg  -- Make sure that concurrent output is not scrambled.-atomicAction :: HasLock e => IO () -> Logger e ()+atomicAction :: (HasLock e) => IO () -> Logger e () atomicAction a = do   l <- reader getLock   liftIO $ withMVar l (const a)

@@ -141,7 +141,7 @@     f j = fromIntegral j / fromIntegral k1  sampleAtPoint ::-  ToJSON a =>+  (ToJSON a) =>   Bool ->   Point ->   Settings ->@@ -187,7 +187,7 @@     a' = MHG ch'  traversePoints ::-  ToJSON a =>+  (ToJSON a) =>   NPoints ->   [(Int, Point)] ->   Settings ->@@ -237,7 +237,7 @@ chop _ _ = error "chop: not all list elements handled"  mlRunPar ::-  ToJSON a =>+  (ToJSON a) =>   ParallelizationMode ->   NPoints ->   [(Int, Point)] ->@@ -269,7 +269,7 @@   pure $ concat xss  mlRun ::-  ToJSON a =>+  (ToJSON a) =>   NPoints ->   [(Int, Point)] ->   ExecutionMode ->@@ -320,7 +320,7 @@     go _ _ = 0  tiWrapper ::-  ToJSON a =>+  (ToJSON a) =>   MLSettings ->   PriorFunction a ->   LikelihoodFunction a ->@@ -404,7 +404,7 @@ --                 llhsNormed = VU.map (\x -> x - llhMax) llhs --                 lhsNormedPowered = VU.map (\x -> exp $ dbeta * x) llhsNormed sssWrapper ::-  ToJSON a =>+  (ToJSON a) =>   MLSettings ->   PriorFunction a ->   LikelihoodFunction a ->@@ -429,7 +429,7 @@  -- | Estimate the marginal likelihood. marginalLikelihood ::-  ToJSON a =>+  (ToJSON a) =>   MLSettings ->   PriorFunction a ->   LikelihoodFunction a ->

@@ -40,7 +40,7 @@ -- transforming the state @a@. type MCMC = ReaderT (Environment Settings) IO -mcmcExecute :: Algorithm a => a -> MCMC a+mcmcExecute :: (Algorithm a) => a -> MCMC a mcmcExecute a = do   logDebugB "Executing MCMC run."   s <- reader settings@@ -51,12 +51,12 @@   logDebugB "Executed MCMC run."   return a' -mcmcResetAcceptance :: Algorithm a => a -> MCMC a+mcmcResetAcceptance :: (Algorithm a) => a -> MCMC a mcmcResetAcceptance a = do   logDebugB "Reset acceptance rates."   return $ aResetAcceptance ResetEverything a -mcmcExceptionHandler :: Algorithm a => Environment Settings -> a -> AsyncException -> IO b+mcmcExceptionHandler :: (Algorithm a) => Environment Settings -> a -> AsyncException -> IO b mcmcExceptionHandler e a err = do   _ <- runReaderT action e   putStrLn "Graceful termination successful."@@ -69,7 +69,7 @@       logWarnS "Press CTRL-C (again) to terminate now."       mcmcClose a -mcmcExecuteMonitors :: Algorithm a => a -> MCMC ()+mcmcExecuteMonitors :: (Algorithm a) => a -> MCMC () mcmcExecuteMonitors a = do   e <- ask   let s = settings e@@ -87,46 +87,46 @@   | IntermediateTuningOff   deriving (Eq) -mcmcIterate :: Algorithm a => IntermediateTuningSpec -> IterationMode -> Int -> a -> MCMC a-mcmcIterate t m n a-  | n < 0 = error "mcmcIterate: Number of iterations is negative."-  | n == 0 = return a-  | otherwise = do-      e <- ask-      let p = sParallelizationMode $ settings e-      -- NOTE: Handle interrupts during iterations, before writing monitors,-      -- using the old algorithm state @a@.-      let handlerOld = mcmcExceptionHandler e a-          maybeIntermediateAutoTune x =-            -- Do not perform intermediate tuning at the last step, because a-            -- normal tuning will be performed.-            case t of-              IntermediateTuningFastProposalsOnlyOn-                | n > 1 ->-                    aAutoTune IntermediateTuningFastProposalsOnly 1 x-                      <&> aResetAcceptance ResetExpectedRatesOnly-              IntermediateTuningAllProposalsOn-                | n > 1 ->-                    aAutoTune IntermediateTuningAllProposals 1 x-                      <&> aResetAcceptance ResetExpectedRatesOnly-              _ -> pure x-          actionIterate = aIterate m p a >>= maybeIntermediateAutoTune-      a' <- liftIO $ actionIterate `catch` handlerOld-      -- NOTE: Mask asynchronous exceptions while writing monitor files. Handle-      -- interrupts after writing monitors; use the new state @a'@.-      ---      -- The problem that arises using this method is: What if executing the-      -- monitors actually throws an error (and not the user or the operating-      -- system that want to stop the chain). In this case, the chain is left in-      -- an undefined state because the monitor files are partly written; the-      -- new state is saved by the handler. However, I do not think I can-      -- recover from partly written monitor files.-      let handlerNew = mcmcExceptionHandler e a'-          actionWrite = runReaderT (mcmcExecuteMonitors a') e-      liftIO $ uninterruptibleMask_ actionWrite `catch` handlerNew-      mcmcIterate t m (n - 1) a'+mcmcIterate :: (Algorithm a) => IntermediateTuningSpec -> IterationMode -> Int -> a -> MCMC a+mcmcIterate t m n a = case n `compare` 0 of+  LT -> error "mcmcIterate: Number of iterations is negative."+  EQ -> return a+  GT -> do+    e <- ask+    let p = sParallelizationMode $ settings e+    -- NOTE: Handle interrupts during iterations, before writing monitors,+    -- using the old algorithm state @a@.+    let handlerOld = mcmcExceptionHandler e a+        maybeIntermediateAutoTune x =+          -- Do not perform intermediate tuning at the last step, because a+          -- normal tuning will be performed.+          case t of+            IntermediateTuningFastProposalsOnlyOn+              | n > 1 ->+                  aAutoTune IntermediateTuningFastProposalsOnly 1 x+                    <&> aResetAcceptance ResetExpectedRatesOnly+            IntermediateTuningAllProposalsOn+              | n > 1 ->+                  aAutoTune IntermediateTuningAllProposals 1 x+                    <&> aResetAcceptance ResetExpectedRatesOnly+            _otherTuningSpecs -> pure x+        actionIterate = aIterate m p a >>= maybeIntermediateAutoTune+    a' <- liftIO $ actionIterate `catch` handlerOld+    -- NOTE: Mask asynchronous exceptions while writing monitor files. Handle+    -- interrupts after writing monitors; use the new state @a'@.+    --+    -- The problem that arises using this method is: What if executing the+    -- monitors actually throws an error (and not the user or the operating+    -- system that want to stop the chain). In this case, the chain is left in+    -- an undefined state because the monitor files are partly written; the+    -- new state is saved by the handler. However, I do not think I can+    -- recover from partly written monitor files.+    let handlerNew = mcmcExceptionHandler e a'+        actionWrite = runReaderT (mcmcExecuteMonitors a') e+    liftIO $ uninterruptibleMask_ actionWrite `catch` handlerNew+    mcmcIterate t m (n - 1) a' -mcmcNewRun :: Algorithm a => a -> MCMC a+mcmcNewRun :: (Algorithm a) => a -> MCMC a mcmcNewRun a = do   s <- reader settings   logInfoB "Starting new MCMC sampler."@@ -145,7 +145,7 @@   logInfoB $ aStdMonitorHeader a''   mcmcIterate IntermediateTuningOff AllProposals i a'' -mcmcContinueRun :: Algorithm a => a -> MCMC a+mcmcContinueRun :: (Algorithm a) => a -> MCMC a mcmcContinueRun a = do   s <- reader settings   let iBurnIn = burnInIterations (sBurnIn s)@@ -164,7 +164,7 @@   logInfoB $ aStdMonitorHeader a   mcmcIterate IntermediateTuningOff AllProposals di a -mcmcBurnIn :: Algorithm a => a -> MCMC a+mcmcBurnIn :: (Algorithm a) => a -> MCMC a mcmcBurnIn a = do   s <- reader settings   case sBurnIn s of@@ -213,7 +213,7 @@       return a''  -- Auto tune the proposals.-mcmcAutotune :: Algorithm a => TuningType -> Int -> a -> MCMC a+mcmcAutotune :: (Algorithm a) => TuningType -> Int -> a -> MCMC a mcmcAutotune t n a = do   case t of     NormalTuningFastProposalsOnly -> logDebugB "Normal auto tune; fast proposals only."@@ -224,7 +224,7 @@     LastTuningAllProposals -> logDebugB "Last auto tune; all proposals."   liftIO $ aAutoTune t n a -mcmcBurnInWithAutoTuning :: Algorithm a => IterationMode -> [Int] -> a -> MCMC a+mcmcBurnInWithAutoTuning :: (Algorithm a) => IterationMode -> [Int] -> a -> MCMC a mcmcBurnInWithAutoTuning _ [] _ = error "mcmcBurnInWithAutoTuning: Empty list." mcmcBurnInWithAutoTuning m [x] a = do   -- Last round.@@ -248,7 +248,7 @@   a''' <- mcmcResetAcceptance a''   mcmcBurnInWithAutoTuning m xs a''' -mcmcInitialize :: Algorithm a => a -> MCMC a+mcmcInitialize :: (Algorithm a) => a -> MCMC a mcmcInitialize a = do   logInfoS $ aName a ++ " algorithm."   s <- settings <$> ask@@ -258,7 +258,7 @@   return a'  -- Save the MCMC run.-mcmcSave :: Algorithm a => a -> MCMC ()+mcmcSave :: (Algorithm a) => a -> MCMC () mcmcSave a = do   s <- reader settings   case sSaveMode s of@@ -273,7 +273,7 @@       logInfoB "Markov chain saved. Analysis can be continued."  -- Report and finish up.-mcmcClose :: Algorithm a => a -> MCMC a+mcmcClose :: (Algorithm a) => a -> MCMC a mcmcClose a = do   logInfoS "Closing monitors."   a' <- liftIO $ aCloseMonitors a@@ -285,7 +285,7 @@   return a'  -- Initialize the run, execute the run, and close the run.-mcmcRun :: Algorithm a => a -> MCMC a+mcmcRun :: (Algorithm a) => a -> MCMC a mcmcRun a = do   -- Header.   logInfoHeader@@ -304,7 +304,7 @@   mcmcClose a''  -- | Run an MCMC algorithm with given settings.-mcmc :: Algorithm a => Settings -> a -> IO a+mcmc :: (Algorithm a) => Settings -> a -> IO a mcmc s a = do   settingsCheck s $ aIteration a   e <- initializeEnvironment s@@ -322,7 +322,7 @@ -- - 'Mcmc.Algorithm.MHG.mhgLoad' -- -- - 'Mcmc.Algorithm.MC3.mc3Load'-mcmcContinue :: Algorithm a => Iterations -> Settings -> a -> IO a+mcmcContinue :: (Algorithm a) => Iterations -> Settings -> a -> IO a mcmcContinue dn s = mcmc s'   where     n' = Iterations $ fromIterations (sIterations s) + fromIterations dn

@@ -55,7 +55,7 @@ instance Contravariant MonitorParameterBatch where   contramap f (MonitorParameterBatch n m) = MonitorParameterBatch n (m . VB.map f) -mean :: Real a => VB.Vector a -> Double+mean :: (Real a) => VB.Vector a -> Double mean xs = realToFrac (VB.sum xs) / fromIntegral (VB.length xs) {-# SPECIALIZE mean :: VB.Vector Double -> Double #-} {-# SPECIALIZE mean :: VB.Vector Int -> Double #-}@@ -64,7 +64,7 @@ -- -- Print the mean with eight decimal places (half precision). monitorBatchMean ::-  Real a =>+  (Real a) =>   -- | Name.   String ->   MonitorParameterBatch a@@ -76,7 +76,7 @@ -- -- Print the mean with full precision. monitorBatchMeanF ::-  Real a =>+  (Real a) =>   -- | Name.   String ->   MonitorParameterBatch a@@ -88,7 +88,7 @@ -- -- Print the real float parameters such as 'Double' with scientific notation. monitorBatchMeanS ::-  Real a =>+  (Real a) =>   -- | Name.   String ->   MonitorParameterBatch a

@@ -46,5 +46,5 @@     ts = round dt  -- | Render a time stamp.-renderTime :: FormatTime t => t -> String+renderTime :: (FormatTime t) => t -> String renderTime = formatTime defaultTimeLocale "%B %-e, %Y, at %H:%M %P, %Z."

@@ -59,38 +59,38 @@ type PriorFunctionG a b = a -> Log b  -- | Flat prior function. Useful for testing and debugging.-noPrior :: RealFloat b => PriorFunctionG a b+noPrior :: (RealFloat b) => PriorFunctionG a b noPrior = const 1.0 {-# SPECIALIZE noPrior :: PriorFunction Double #-}  -- | Improper uniform prior; strictly greater than a given value.-greaterThan :: RealFloat a => LowerBoundary a -> PriorFunctionG a a+greaterThan :: (RealFloat a) => LowerBoundary a -> PriorFunctionG a a greaterThan a x   | x > a = 1.0   | otherwise = 0.0 {-# SPECIALIZE greaterThan :: Double -> PriorFunction Double #-}  -- | Improper uniform prior; strictly greater than zero.-positive :: RealFloat a => PriorFunctionG a a+positive :: (RealFloat a) => PriorFunctionG a a positive = greaterThan 0 {-# SPECIALIZE positive :: PriorFunction Double #-}  -- | Improper uniform prior; strictly less than a given value.-lessThan :: RealFloat a => UpperBoundary a -> PriorFunctionG a a+lessThan :: (RealFloat a) => UpperBoundary a -> PriorFunctionG a a lessThan a x   | x < a = 1.0   | otherwise = 0.0 {-# SPECIALIZE lessThan :: Double -> PriorFunction Double #-}  -- | Improper uniform prior; strictly less than zero.-negative :: RealFloat a => PriorFunctionG a a+negative :: (RealFloat a) => PriorFunctionG a a negative = lessThan 0.0 {-# SPECIALIZE negative :: PriorFunction Double #-}  -- | Exponential distributed prior. -- -- Call 'error' if the rate is zero or negative.-exponential :: RealFloat a => Rate a -> PriorFunctionG a a+exponential :: (RealFloat a) => Rate a -> PriorFunctionG a a exponential l x   | l <= 0.0 = error "exponential: Rate is zero or negative."   | x <= 0.0 = 0.0@@ -136,17 +136,17 @@  -- | Calculate mean and variance of the gamma distribution given the shape and -- the scale.-gammaShapeScaleToMeanVariance :: Num a => Shape a -> Scale a -> (Mean a, Variance a)+gammaShapeScaleToMeanVariance :: (Num a) => Shape a -> Scale a -> (Mean a, Variance a) gammaShapeScaleToMeanVariance k t = let m = k * t in (m, m * t) {-# SPECIALIZE gammaShapeScaleToMeanVariance :: Double -> Double -> (Double, Double) #-}  -- | Calculate shape and scale of the gamma distribution given the mean and -- the variance.-gammaMeanVarianceToShapeScale :: Fractional a => Mean a -> Variance a -> (Shape a, Scale a)+gammaMeanVarianceToShapeScale :: (Fractional a) => Mean a -> Variance a -> (Shape a, Scale a) gammaMeanVarianceToShapeScale m v = (m * m / v, v / m) {-# SPECIALIZE gammaMeanVarianceToShapeScale :: Double -> Double -> (Double, Double) #-} -mLnSqrt2Pi :: RealFloat a => a+mLnSqrt2Pi :: (RealFloat a) => a mLnSqrt2Pi = 0.9189385332046727417803297364056176398613974736377834128171 {-# INLINE mLnSqrt2Pi #-} @@ -158,7 +158,7 @@ -- \(\mu\) and \(\sigma\), but are not the same as \(\mu\) and \(\sigma\)! -- -- Call 'error' if the standard deviation is zero or negative.-logNormal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a+logNormal :: (RealFloat a) => Mean a -> StandardDeviation a -> PriorFunctionG a a logNormal m s x   | s <= 0.0 = error "logNormal: Standard deviation is zero or negative."   | x <= 0.0 = 0.0@@ -172,7 +172,7 @@ -- | Normal distributed prior. -- -- Call 'error' if the standard deviation is zero or negative.-normal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a+normal :: (RealFloat a) => Mean a -> StandardDeviation a -> PriorFunctionG a a normal m s x   | s <= 0 = error "normal: Standard deviation is zero or negative."   | otherwise = Exp $ (-xm * xm / (2 * s * s)) - denom@@ -184,7 +184,7 @@ -- | Uniform prior on [a, b]. -- -- Call 'error' if the lower boundary is greather than the upper boundary.-uniform :: RealFloat a => LowerBoundary a -> UpperBoundary a -> PriorFunctionG a a+uniform :: (RealFloat a) => LowerBoundary a -> UpperBoundary a -> PriorFunctionG a a uniform a b x   | a > b = error "uniform: Lower boundary is greater than upper boundary."   | x < a = 0.0@@ -205,11 +205,11 @@ -- -- Use with care because the elements are checked for positiveness, and this can -- take some time if the list is long and does not contain any zeroes.-product' :: RealFloat a => [Log a] -> Log a+product' :: (RealFloat a) => [Log a] -> Log a product' = fromMaybe 0 . prodM {-# SPECIALIZE product' :: [Log Double] -> Log Double #-}  -- The type could be generalized to any MonadPlus Integer-prodM :: RealFloat a => [Log a] -> Maybe (Log a)+prodM :: (RealFloat a) => [Log a] -> Maybe (Log a) prodM = foldM (\ !acc x -> (acc * x) <$ guard (acc /= 0.0)) 1.0 {-# SPECIALIZE prodM :: [Log Double] -> Maybe (Log Double) #-}

@@ -97,7 +97,7 @@ defaultHParams = HParams Nothing Nothing Nothing  hamiltonianPFunctionWithTuningParameters ::-  Traversable s =>+  (Traversable s) =>   Dimension ->   HStructure s ->   (s Double -> Target) ->@@ -159,7 +159,7 @@ -- -- May call 'error' during initialization. hamiltonian ::-  Traversable s =>+  (Traversable s) =>   HParams ->   HTuningConf ->   HStructure s ->

@@ -235,7 +235,7 @@  -- See Algorithm 4 in [4]. findReasonableEpsilon ::-  StatefulGen g m =>+  (StatefulGen g m) =>   Target ->   Masses ->   Positions ->@@ -325,7 +325,7 @@   where     err msg = error $ "hTuningFunctionWith: " <> msg -checkHStructureWith :: Foldable s => Masses -> HStructure s -> Maybe String+checkHStructureWith :: (Foldable s) => Masses -> HStructure s -> Maybe String checkHStructureWith ms (HStructure x toVec fromVec)   | toList (fromVec x xVec) /= toList x = eWith "'fromVectorWith x (toVector x) /= x' for sample state."   | L.size xVec /= nrows = eWith "Mass matrix and 'toVector x' have different sizes for sample state."@@ -337,7 +337,7 @@  -- Generate momenta for a new iteration. generateMomenta ::-  StatefulGen g m =>+  (StatefulGen g m) =>   Mu ->   Masses ->   g ->

@@ -171,7 +171,7 @@ defaultNParams = NParams Nothing Nothing  nutsPFunctionWithTuningParameters ::-  Traversable s =>+  (Traversable s) =>   Dimension ->   HStructure s ->   (s Double -> Target) ->@@ -264,7 +264,7 @@ -- -- May call 'error' during initialization. nuts ::-  Traversable s =>+  (Traversable s) =>   NParams ->   HTuningConf ->   HStructure s ->