mcmc 0.6.2.3 → 0.6.2.4

raw patch · 16 files changed

+172/−99 lines, 16 filesdep ~covariancePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: covariance

API changes (from Hackage documentation)

Files

• ChangeLog.md +16/−0
• mcmc.cabal +2/−2
• src/Mcmc/Algorithm.hs +3/−2
• src/Mcmc/Algorithm/MC3.hs +7/−6
• src/Mcmc/Algorithm/MHG.hs +6/−6
• src/Mcmc/Cycle.hs +26/−10
• src/Mcmc/MarginalLikelihood.hs +7/−6
• src/Mcmc/Mcmc.hs +39/−26
• src/Mcmc/Proposal.hs +22/−8
• src/Mcmc/Proposal/Bactrian.hs +2/−2
• src/Mcmc/Proposal/Hamiltonian.hs +10/−10
• src/Mcmc/Proposal/Scale.hs +3/−3
• src/Mcmc/Proposal/Simplex.hs +2/−2
• src/Mcmc/Proposal/Slide.hs +4/−4
• src/Mcmc/Settings.hs +20/−9
• test/Mcmc/ProposalSpec.hs +3/−3

@@ -5,6 +5,22 @@ ## Unreleased changes  +## 0.6.2.4++-   Specify covariance version bounds. Use `covariance-0.2.0.0` (specifically+    state that sigma is rescaled with `rescaleSWith`).+++## 0.6.2.3++-   Allow burn in with fast proposals only (`BurnInWithCustomAutoTUning`).+    Sometimes it is advantageous to hold back slow proposals initially, especially+    when the state is so far off that it does not make sense to compute complex+    proposals.+-   Hamiltonian proposal: Use automatic differentiation specialized to `Double`+    (roughly 10 percent faster).++ ## 0.6.2.2  -   Remove dependency `monad-parallel`. Fix stackage build.

@@ -1,6 +1,6 @@ cabal-version:      3.0 name:               mcmc-version:            0.6.2.3+version:            0.6.2.4 synopsis:           Sample from a posterior using Markov chain Monte Carlo description:   Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme>@@ -90,7 +90,7 @@     , bytestring     , circular     , containers-    , covariance+    , covariance         >=0.2     , data-default     , deepseq     , directory

@@ -16,6 +16,7 @@  import qualified Data.ByteString.Lazy.Char8 as BL import Data.Time+import Mcmc.Cycle import Mcmc.Settings  -- | Class for algorithms used by MCMC samplers.@@ -31,7 +32,7 @@   aIsInValidState :: a -> Bool    -- | Sample the next state.-  aIterate :: ParallelizationMode -> a -> IO a+  aIterate :: IterationMode -> ParallelizationMode -> a -> IO a    -- | Auto tune all proposals over the last N iterations.   --@@ -43,7 +44,7 @@   aResetAcceptance :: a -> a    -- | Summarize the cycle.-  aSummarizeCycle :: a -> BL.ByteString+  aSummarizeCycle :: IterationMode -> a -> BL.ByteString    -- | Open required monitor files and setup corresponding file handles.   aOpenMonitors :: AnalysisName -> ExecutionMode -> a -> IO a

@@ -426,10 +426,11 @@ -- However, we have to take care of the mutable traces. mc3Iterate ::   ToJSON a =>+  IterationMode ->   ParallelizationMode ->   MC3 a ->   IO (MC3 a)-mc3Iterate pm a = do+mc3Iterate m pm a = do   -- 1. Maybe propose swaps.   --   -- NOTE: Swaps have to be proposed first, because the traces are automatically@@ -446,10 +447,10 @@       else return a   -- 2. Iterate all chains and increment iteration.   mhgs <- case pm of-    Sequential -> V.mapM (aIterate pm) (mc3MHGChains a')+    Sequential -> V.mapM (aIterate m pm) (mc3MHGChains a')     Parallel ->       -- Go via a list, and use 'forkIO' ("Control.Concurrent.Async").-      V.fromList <$> mapConcurrently (aIterate pm) (V.toList $ mc3MHGChains a')+      V.fromList <$> mapConcurrently (aIterate m pm) (V.toList $ mc3MHGChains a')   let i = mc3Iteration a'   return $ a' {mc3MHGChains = mhgs, mc3Iteration = succ i} @@ -526,8 +527,8 @@ -- - 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 => MC3 a -> BL.ByteString-mc3SummarizeCycle a =+mc3SummarizeCycle :: ToJSON a => IterationMode -> MC3 a -> BL.ByteString+mc3SummarizeCycle m a =   BL.intercalate "\n" $     [ "MC3: Cycle of cold chain.",       coldMHGCycleSummary@@ -558,7 +559,7 @@       ++ [proposalHLine]   where     mhgs = mc3MHGChains a-    coldMHGCycleSummary = aSummarizeCycle $ V.head mhgs+    coldMHGCycleSummary = aSummarizeCycle m $ V.head mhgs     cs = V.map fromMHG mhgs     -- Acceptance rates may be 'Nothing' when no proposals have been undertaken.     -- The 'sequence' operations pull the 'Nothing's out of the inner

@@ -96,7 +96,7 @@       TraceMinimum n -> n     bi = case sBurnIn s of       BurnInWithAutoTuning _ n -> n-      BurnInWithCustomAutoTuning ns -> maximum ns+      BurnInWithCustomAutoTuning ns ms -> max (maximum $ 0 : ns) (maximum $ 0 : ms)       _ -> 0     traceLength = maximum $ minimumTraceLength : bi : batchMonitorSizes @@ -236,9 +236,9 @@ -- algorithm is just inherently sequential. Parallelization can be achieved by -- having parallel prior and/or likelihood functions, or by using algorithms -- running parallel chains such as 'MC3'.-mhgIterate :: ParallelizationMode -> MHG a -> IO (MHG a)-mhgIterate _ a = do-  ps <- prepareProposals cc g+mhgIterate :: IterationMode -> ParallelizationMode -> MHG a -> IO (MHG a)+mhgIterate m _ a = do+  ps <- prepareProposals m cc g   a' <- foldM mhgPropose a ps   mhgPush a'   where@@ -259,8 +259,8 @@   where     ac = acceptance c -mhgSummarizeCycle :: MHG a -> BL.ByteString-mhgSummarizeCycle (MHG c) = summarizeCycle ac cc+mhgSummarizeCycle :: IterationMode -> MHG a -> BL.ByteString+mhgSummarizeCycle m (MHG c) = summarizeCycle m ac cc   where     cc = cycle c     ac = acceptance c

@@ -18,6 +18,7 @@     Cycle (ccProposals),     cycleFromList,     setOrder,+    IterationMode (..),     prepareProposals,     autoTuneCycle, @@ -116,9 +117,13 @@ setOrder :: Order -> Cycle a -> Cycle a setOrder o c = c {ccOrder = o} +-- | Use all proposals, or use fast proposals only?+data IterationMode = AllProposals | FastProposals+  deriving (Eq)+ -- | Replicate 'Proposal's according to their weights and possibly shuffle them.-prepareProposals :: Cycle a -> GenIO -> IO [Proposal a]-prepareProposals (Cycle xs o) g = case o of+prepareProposals :: IterationMode -> Cycle a -> GenIO -> IO [Proposal a]+prepareProposals m (Cycle xs o) g = case o of   RandomO -> shuffle ps g   SequentialO -> return ps   RandomReversibleO -> do@@ -126,17 +131,28 @@     return $ psR ++ reverse psR   SequentialReversibleO -> return $ ps ++ reverse ps   where-    !ps = concat [replicate (fromPWeight $ prWeight p) p | p <- xs]+    !ps =+      concat+        [ replicate (fromPWeight $ prWeight p) p+          | p <- xs,+            case m of+              AllProposals -> True+              -- Only use proposal if it is fast.+              FastProposals -> prSpeed p == PFast+        ]  -- The number of proposals depends on the order.-getNProposalsPerCycle :: Cycle a -> Int-getNProposalsPerCycle (Cycle xs o) = case o of+getNProposalsPerCycle :: IterationMode -> Cycle a -> Int+getNProposalsPerCycle m (Cycle xs o) = case o of   RandomO -> once   SequentialO -> once   RandomReversibleO -> 2 * once   SequentialReversibleO -> 2 * once   where-    once = sum $ map (fromPWeight . prWeight) xs+    xs' = case m of+      AllProposals -> xs+      FastProposals -> filter (\x -> prSpeed x == PFast) xs+    once = sum $ map (fromPWeight . prWeight) xs'  -- See 'tuneWithTuningParameters' and 'Tuner'. tuneWithChainParameters :: AcceptanceRate -> VB.Vector a -> Proposal a -> Either String (Proposal a)@@ -167,8 +183,8 @@ proposalHLine = BL.replicate (BL.length proposalHeader) '-'  -- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance rates.-summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> BL.ByteString-summarizeCycle a c =+summarizeCycle :: IterationMode -> Acceptance (Proposal a) -> Cycle a -> BL.ByteString+summarizeCycle m a c =   BL.intercalate "\n" $     [ "Summary of proposal(s) in cycle.",       nProposalsFullStr,@@ -188,9 +204,9 @@       ++ [proposalHLine]   where     ps = ccProposals c-    nProposals = getNProposalsPerCycle c+    nProposals = getNProposalsPerCycle m c     nProposalsStr = BB.toLazyByteString $ BB.intDec nProposals     nProposalsFullStr = case nProposals of       1 -> nProposalsStr <> " proposal is performed per iteration."       _ -> nProposalsStr <> " proposals are performed per iterations."-    ar m = acceptanceRate m a+    ar pr = acceptanceRate pr a

@@ -21,11 +21,11 @@   ) where +import Control.Concurrent.Async hiding (link) import Control.Monad import Control.Monad.IO.Class-import Control.Concurrent.Async hiding (link)-import Control.Monad.Trans.Reader import Control.Monad.Trans.Class+import Control.Monad.Trans.Reader import Data.Aeson import Data.List hiding (cycle) import qualified Data.Map.Strict as M@@ -152,7 +152,7 @@       ac = acceptance ch''       mAr = sequence $ acceptanceRates ac   logDebugB "sampleAtPoint: Summarize cycle."-  logDebugB $ summarizeCycle ac $ cycle ch''+  logDebugB $ summarizeCycle AllProposals ac $ cycle ch''   case mAr of     Nothing -> logWarnB "Some acceptance rates are unavailable. The tuning period may be too small."     Just ar -> do@@ -279,10 +279,11 @@    -- Parallel execution of both path integrals.   r <- ask-  (lhssForward, lhssBackward) <- lift $+  (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 (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)   logInfoEndTime    logDebugB "tiWrapper: Calculate mean log likelihoods."

@@ -26,6 +26,7 @@ import Control.Monad.IO.Class import Control.Monad.Trans.Reader import Mcmc.Algorithm+import Mcmc.Cycle import Mcmc.Environment import Mcmc.Logger import Mcmc.Settings@@ -88,17 +89,17 @@   mStdLog <- liftIO $ mask_ $ aExecuteMonitors vb t0 iTotal a   forM_ mStdLog (logOutB "   ") -mcmcIterate :: Algorithm a => Int -> a -> MCMC a-mcmcIterate n a+mcmcIterate :: Algorithm a => IterationMode -> Int -> a -> MCMC a+mcmcIterate m n a   | n < 0 = error "mcmcIterate: Number of iterations is negative."   | n == 0 = return a   | otherwise = do       e <- ask       p <- sParallelizationMode . settings <$> ask       -- NOTE: User interrupt is handled during iterations.-      a' <- liftIO $ catch (aIterate p a) (mcmcExceptionHandler e a)+      a' <- liftIO $ catch (aIterate m p a) (mcmcExceptionHandler e a)       mcmcExecuteMonitors a'-      mcmcIterate (n - 1) a'+      mcmcIterate m (n - 1) a'  mcmcNewRun :: Algorithm a => a -> MCMC a mcmcNewRun a = do@@ -108,12 +109,11 @@   logInfoB $ aStdMonitorHeader a   mcmcExecuteMonitors a   when (aIsInValidState a) (logWarnB "The initial state is invalid!")-  logInfoB $ aSummarizeCycle a   a' <- mcmcBurnIn a   let i = fromIterations $ sIterations s   logInfoS $ "Run chain for " ++ show i ++ " iterations."   logInfoB $ aStdMonitorHeader a'-  mcmcIterate i a'+  mcmcIterate AllProposals i a'  mcmcContinueRun :: Algorithm a => a -> MCMC a mcmcContinueRun a = do@@ -129,28 +129,31 @@   logInfoS $ "Current iteration: " ++ show iCurrent ++ "."   when (iCurrent < iBurnIn) $ error "mcmcContinueRun: Can not continue burn in."   let di = iTotal - iCurrent-  logInfoB $ aSummarizeCycle a+  logInfoB $ aSummarizeCycle AllProposals a   logInfoS $ "Run chain for " ++ show di ++ " iterations."   logInfoB $ aStdMonitorHeader a-  mcmcIterate di a+  mcmcIterate AllProposals di a  mcmcBurnIn :: Algorithm a => a -> MCMC a mcmcBurnIn a = do   s <- reader settings   case sBurnIn s of     NoBurnIn -> do+      logInfoB $ aSummarizeCycle AllProposals a       logInfoS "No burn in."       return a     BurnInWithoutAutoTuning n -> do+      logInfoB $ aSummarizeCycle AllProposals a       logInfoS $ "Burn in for " <> show n <> " iterations."       logInfoS "Auto tuning is disabled."       logInfoB $ aStdMonitorHeader a-      a' <- mcmcIterate n a-      logInfoB $ aSummarizeCycle a'+      a' <- mcmcIterate AllProposals n a+      logInfoB $ aSummarizeCycle AllProposals a'       a'' <- mcmcResetAcceptance a'       logInfoB "Burn in finished."       return a''     BurnInWithAutoTuning n t -> do+      logInfoB $ aSummarizeCycle AllProposals a       logInfoS $ "Burn in for " ++ show n ++ " iterations."       logInfoS $ "Auto tuning is enabled with a period of " ++ show t ++ "."       logInfoB $ aStdMonitorHeader a@@ -158,16 +161,26 @@           -- Don't add another auto tune period if r == 0, because then we auto           -- tune without acceptance counts and get NaNs.           xs = replicate m t <> [r | r > 0]-      a' <- mcmcBurnInWithAutoTuning xs a+      a' <- mcmcBurnInWithAutoTuning AllProposals xs a       logInfoB "Burn in finished."       return a'-    BurnInWithCustomAutoTuning xs -> do-      logInfoS $ "Burn in for " ++ show (sum xs) ++ " iterations."-      logInfoS $ "Custom auto tuning is enabled with periods " ++ show xs ++ "."+    BurnInWithCustomAutoTuning xs ys -> do+      logInfoS $ "Burn in for " ++ show (sum xs + sum ys) ++ " iterations."+      a' <-+        if null xs+          then do+            logInfoB $ aSummarizeCycle AllProposals a+            pure a+          else do+            logInfoB $ aSummarizeCycle FastProposals a+            logInfoS $ "Fast custom auto tuning with periods " ++ show xs ++ "."+            logInfoB $ aStdMonitorHeader a+            mcmcBurnInWithAutoTuning FastProposals xs a+      logInfoS $ "Full custom auto tuning with periods " ++ show ys ++ "."       logInfoB $ aStdMonitorHeader a-      a' <- mcmcBurnInWithAutoTuning xs a+      a'' <- mcmcBurnInWithAutoTuning AllProposals ys a'       logInfoB "Burn in finished."-      return a'+      return a''  -- Auto tune the proposals. mcmcAutotune :: Algorithm a => Int -> a -> MCMC a@@ -175,23 +188,23 @@   logDebugB "Auto tune."   liftIO $ aAutoTune n a -mcmcBurnInWithAutoTuning :: Algorithm a => [Int] -> a -> MCMC a-mcmcBurnInWithAutoTuning [] _ = error "mcmcBurnInWithAutoTuning: Empty lisst."-mcmcBurnInWithAutoTuning [x] a = do+mcmcBurnInWithAutoTuning :: Algorithm a => IterationMode -> [Int] -> a -> MCMC a+mcmcBurnInWithAutoTuning _ [] _ = error "mcmcBurnInWithAutoTuning: Empty list."+mcmcBurnInWithAutoTuning m [x] a = do   -- Last round.-  a' <- mcmcIterate x a+  a' <- mcmcIterate m x a   a'' <- mcmcAutotune x a'-  logInfoB $ aSummarizeCycle a''+  logInfoB $ aSummarizeCycle m a''   logInfoS $ "Acceptance rates calculated over the last " <> show x <> " iterations."   mcmcResetAcceptance a''-mcmcBurnInWithAutoTuning (x : xs) a = do-  a' <- mcmcIterate x a+mcmcBurnInWithAutoTuning m (x : xs) a = do+  a' <- mcmcIterate m x a   a'' <- mcmcAutotune x a'-  logDebugB $ aSummarizeCycle a''+  logDebugB $ aSummarizeCycle m a''   logDebugS $ "Acceptance rates calculated over the last " <> show x <> " iterations."   logDebugB $ aStdMonitorHeader a''   a''' <- mcmcResetAcceptance a''-  mcmcBurnInWithAutoTuning xs a'''+  mcmcBurnInWithAutoTuning m xs a'''  mcmcInitialize :: Algorithm a => a -> MCMC a mcmcInitialize a = do@@ -221,7 +234,7 @@ mcmcClose :: Algorithm a => a -> MCMC a mcmcClose a = do   logDebugB "Closing MCMC run."-  logInfoB $ aSummarizeCycle a+  logInfoB $ aSummarizeCycle AllProposals a   logInfoS $ aName a ++ " algorithm finished."   mcmcSave a   logInfoEndTime

@@ -20,6 +20,7 @@     PWeight (fromPWeight),     pWeight,     PDimension (..),+    PSpeed (..),     Proposal (..),     KernelRatio,     Jacobian,@@ -90,6 +91,9 @@ -- -- The number of affected, independent parameters. --+-- The dimension is used to calculate the optimal acceptance rate, and does not+-- have to be exact.+-- -- Usually, the optimal acceptance rate of low dimensional proposals is higher -- than for high dimensional ones. However, this is not always true (see below). --@@ -131,6 +135,15 @@   | -- | Provide dimension ('Int') and desired acceptance rate ('Double').     PSpecial Int Double +-- | Rough indication whether a proposal is fast or slow.+--+-- Useful during burn in. Slow proposals are not executed during fast auto+-- tuning periods.+--+-- See 'Mcmc.Settings.BurnInSettings'.+data PSpeed = PFast | PSlow+  deriving (Eq)+ -- | A 'Proposal' is an instruction about how the Markov chain will traverse the -- state space @a@. Essentially, it is a probability mass or probability density -- conditioned on the current state (i.e., a Markov kernel).@@ -146,8 +159,7 @@     prName :: PName,     -- | Description of the proposal type and parameters.     prDescription :: PDescription,-    -- | Dimension of the proposal. The dimension is used to calculate the-    -- optimal acceptance rate, and does not have to be exact.+    prSpeed :: PSpeed,     prDimension :: PDimension,     -- | The weight determines how often a 'Proposal' is executed per iteration of     -- the Markov chain.@@ -212,8 +224,8 @@ -- For further reference, please see the [example -- @Pair@](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Pair/Pair.hs). liftProposalWith :: JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b-liftProposalWith jf l (Proposal n r d w s t) =-  Proposal n r d w (liftProposalSimpleWith jf l s) (liftTunerWith jf l <$> t)+liftProposalWith jf l (Proposal n r d p w s t) =+  Proposal n r d p w (liftProposalSimpleWith jf l s) (liftTunerWith jf l <$> t)  -- | Simple proposal without tuning information. --@@ -256,6 +268,8 @@   PDescription ->   -- | Function creating a simple proposal for a given tuning parameter.   (TuningParameter -> ProposalSimple a) ->+  -- | Speed.+  PSpeed ->   -- | Dimension.   PDimension ->   -- | Name.@@ -265,15 +279,15 @@   -- | Activate tuning?   Tune ->   Proposal a-createProposal r f d n w Tune =-  Proposal n r d w (f 1.0) (Just tuner)+createProposal r f s d n w Tune =+  Proposal n r s d w (f 1.0) (Just tuner)   where     fT = defaultTuningFunctionWith d     fTs = noAuxiliaryTuningFunction     g t _ = Right $ f t     tuner = Tuner 1.0 fT VU.empty fTs g-createProposal r f d n w NoTune =-  Proposal n r d w (f 1.0) Nothing+createProposal r f s d n w NoTune =+  Proposal n r s d w (f 1.0) Nothing  -- | Required information to tune 'Proposal's. data Tuner a = Tuner

@@ -90,7 +90,7 @@   PWeight ->   Tune ->   Proposal Double-slideBactrian m s = createProposal description (bactrianAdditiveSimple m s) (PDimension 1)+slideBactrian m s = createProposal description (bactrianAdditiveSimple m s) PFast (PDimension 1)   where     description = PDescription $ "Slide Bactrian; spike: " ++ show m ++ ", sd: " ++ show s @@ -137,6 +137,6 @@   PWeight ->   Tune ->   Proposal Double-scaleBactrian m s = createProposal description (bactrianMultSimple m s) (PDimension 1)+scaleBactrian m s = createProposal description (bactrianMultSimple m s) PFast (PDimension 1)   where     description = PDescription $ "Scale Bactrian; spike: " ++ show m <> ", sd: " <> show s

@@ -212,7 +212,10 @@  -- | Specifications of the Hamilton Monte Carlo proposal. data HSettings a = HSettings-  { -- | Extract values to be manipulated by the Hamiltonian proposal from the+  { -- | The sample state is used for error checks and to calculate the dimension+    -- of the proposal.+    hSample :: a,+    -- | Extract values to be manipulated by the Hamiltonian proposal from the     -- state.     hToVector :: a -> Values,     -- | Put those values back into the state.@@ -225,8 +228,8 @@     hTune :: HTune   } -checkHSettings :: Eq a => a -> HSettings a -> Maybe String-checkHSettings x (HSettings toVec fromVec _ _ masses l eps _)+checkHSettings :: Eq a => HSettings a -> Maybe String+checkHSettings (HSettings x toVec fromVec _ _ masses l eps _)   | any (<= 0) diagonalMasses = Just "checkHSettings: Some diagonal entries of the mass matrix are zero or negative."   | nrows /= ncols = Just "checkHSettings: Mass matrix is not square."   | fromVec x xVec /= x = Just "checkHSettings: 'fromVectorWith x (toVector x) /= x' for sample state."@@ -456,7 +459,7 @@           kernelR = prPhi' / prPhi        in return (fromVec x theta', kernelR, 1.0)   where-    (HSettings toVec fromVec gradient mVal masses l e _) = st+    (HSettings _ toVec fromVec gradient mVal masses l e _) = st     theta = toVec x     lL = maximum [1 :: Int, floor $ (0.8 :: Double) * fromIntegral l]     lR = maximum [lL, ceiling $ (1.2 :: Double) * fromIntegral l]@@ -569,23 +572,20 @@ -- | Hamiltonian Monte Carlo proposal. hamiltonian ::   Eq a =>-  -- | The sample state is used for error checks and to calculate the dimension-  -- of the proposal.-  a ->   HSettings a ->   PName ->   PWeight ->   Proposal a-hamiltonian x s n w = case checkHSettings x s of+hamiltonian s n w = case checkHSettings s of   Just err -> error err   Nothing ->     let desc = PDescription "Hamiltonian Monte Carlo (HMC)"         toVec = hToVector s-        dim = (L.size $ toVec x)+        dim = (L.size $ toVec $ hSample s)         pDim = PSpecial dim 0.65         ts = massesToTuningParameters (hMasses s)         ps = hamiltonianSimple s-        hamiltonianWith = Proposal n desc pDim w ps+        hamiltonianWith = Proposal n desc PSlow pDim w ps         tSet@(HTune tlf tms) = hTune s         tFun = case tlf of           HNoTuneLeapfrog -> noTuningFunction

@@ -52,7 +52,7 @@   PWeight ->   Tune ->   Proposal Double-scale k th = createProposal description (scaleSimple k th) (PDimension 1)+scale k th = createProposal description (scaleSimple k th) PFast (PDimension 1)   where     description = PDescription $ "Scale; shape: " ++ show k ++ ", scale: " ++ show th @@ -66,7 +66,7 @@   PWeight ->   Tune ->   Proposal Double-scaleUnbiased k = createProposal description (scaleSimple k (1 / k)) (PDimension 1)+scaleUnbiased k = createProposal description (scaleSimple k (1 / k)) PFast (PDimension 1)   where     description = PDescription $ "Scale unbiased; shape: " ++ show k @@ -96,6 +96,6 @@   PWeight ->   Tune ->   Proposal (Double, Double)-scaleContrarily k th = createProposal description (scaleContrarilySimple k th) (PDimension 2)+scaleContrarily k th = createProposal description (scaleContrarilySimple k th) PFast (PDimension 2)   where     description = PDescription $ "Scale contrariliy; shape: " ++ show k ++ ", scale: " ++ show th

@@ -140,7 +140,7 @@ -- For high dimensional simplices, this proposal may have low acceptance rates. -- In this case, please see the coordinate wise 'beta' proposal. dirichlet :: PDimension -> PName -> PWeight -> Tune -> Proposal Simplex-dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple+dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple PFast  -- The tuning parameter is the inverted mean of the shape values. --@@ -206,6 +206,6 @@ -- This proposal has been assigned a dimension of 2. See the discussion at -- 'PDimension'. beta :: Dimension -> PName -> PWeight -> Tune -> Proposal Simplex-beta i = createProposal description (betaSimple i) (PDimension 2)+beta i = createProposal description (betaSimple i) PFast (PDimension 2)   where     description = PDescription $ "Beta; coordinate: " ++ show i

@@ -40,7 +40,7 @@   PWeight ->   Tune ->   Proposal Double-slide m s = createProposal description (slideSimple m s) (PDimension 1)+slide m s = createProposal description (slideSimple m s) PFast (PDimension 1)   where     description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s @@ -60,7 +60,7 @@   PWeight ->   Tune ->   Proposal Double-slideSymmetric s = createProposal description (slideSymmetricSimple s) (PDimension 1)+slideSymmetric s = createProposal description (slideSymmetricSimple s) PFast (PDimension 1)   where     description = PDescription $ "Slide symmetric; sd: " ++ show s @@ -80,7 +80,7 @@   PWeight ->   Tune ->   Proposal Double-slideUniformSymmetric d = createProposal description (slideUniformSimple d) (PDimension 1)+slideUniformSymmetric d = createProposal description (slideUniformSimple d) PFast (PDimension 1)   where     description = PDescription $ "Slide uniform symmetric; delta: " ++ show d @@ -108,6 +108,6 @@   PWeight ->   Tune ->   Proposal (Double, Double)-slideContrarily m s = createProposal description (slideContrarilySimple m s) (PDimension 2)+slideContrarily m s = createProposal description (slideContrarilySimple m s) PFast (PDimension 2)   where     description = PDescription $ "Slide contrarily; mean: " ++ show m ++ ", sd: " ++ show s

@@ -69,15 +69,24 @@   | -- | Burn in for a given number of iterations. Enable auto tuning with a     -- given period.     BurnInWithAutoTuning Int Int-  | -- | Burn in with the given list of auto tuning periods.+  | -- | Burn in with the given list of fast and full auto tuning periods.     ---    -- For example, @BurnInWithCustomAutoTuning [100,200]@ performs 300-    -- iterations with two auto tuning steps. One after 100 iterations, the-    -- second one after 200 more iterations.+    -- The list of fast auto tuning periods may be empty. All periods have to be+    -- strictly positive.     --+    -- See also 'Mcmc.Proposals.PSpeed'.+    --+    -- For example, @BurnInWithCustomAutoTuning [50] [100,200]@ performs+    -- 1a. 50 iterations without any slow proposals such as Hamiltonian proposals;+    -- 1b. Auto tuning;+    -- 2a. 100 iterations with all proposals;+    -- 2b Auto tuning;+    -- 3a. 200 iterations with all proposals;+    -- 3b. Auto tuning.+    --     -- Usually it is useful to auto tune more frequently in the beginning of the     -- MCMC run.-    BurnInWithCustomAutoTuning [Int]+    BurnInWithCustomAutoTuning [Int] [Int]   deriving (Eq, Read, Show)  $(deriveJSON defaultOptions ''BurnInSettings)@@ -89,22 +98,24 @@   bsInt x <> " iterations; no auto tune." burnInPrettyPrint (BurnInWithAutoTuning x y) =   bsInt x <> " iterations; auto tune with a period of " <> bsInt y <> "."-burnInPrettyPrint (BurnInWithCustomAutoTuning xs) =-  bsInt (sum xs) <> " iterations; custom auto tune periods."+burnInPrettyPrint (BurnInWithCustomAutoTuning xs ys) =+  bsInt (sum xs) <> " fast," <> bsInt (sum ys) <> " slow iterations; custom auto tune periods."  -- Check if the burn in settings are valid. burnInValid :: BurnInSettings -> Bool burnInValid NoBurnIn = True burnInValid (BurnInWithoutAutoTuning n) = n > 0 burnInValid (BurnInWithAutoTuning n t) = n > 0 && t > 0-burnInValid (BurnInWithCustomAutoTuning xs) = not (null xs) && all (> 0) xs+-- The list of fast auto tuning periods may be empty, the list of full auto+-- tuning periods must be non-empty. All periods have to be strictly positive.+burnInValid (BurnInWithCustomAutoTuning xs ys) = all (> 0) xs && not (null ys) && all (> 0) ys  -- | Get the number of burn in iterations. burnInIterations :: BurnInSettings -> Int burnInIterations NoBurnIn = 0 burnInIterations (BurnInWithoutAutoTuning n) = n burnInIterations (BurnInWithAutoTuning n _) = n-burnInIterations (BurnInWithCustomAutoTuning xs) = sum xs+burnInIterations (BurnInWithCustomAutoTuning xs ys) = sum xs + sum ys  -- | Number of normal iterations after burn in. --

@@ -35,10 +35,10 @@     it "returns the correct number of proposals in a cycle" $       do         g <- create-        l1 <- length <$> prepareProposals c g+        l1 <- length <$> prepareProposals AllProposals c g         l1 `shouldBe` 4-        l2 <- length <$> prepareProposals (setOrder RandomReversibleO c) g+        l2 <- length <$> prepareProposals AllProposals (setOrder RandomReversibleO c) g         l2 `shouldBe` 8-        o3 <- prepareProposals (setOrder SequentialReversibleO c) g+        o3 <- prepareProposals AllProposals (setOrder SequentialReversibleO c) g         length o3 `shouldBe` 8         o3 == [p1, p2, p2, p2, p2, p2, p2, p1] `shouldBe` True