mcmc 0.8.0.1 → 0.8.1.0

raw patch · 21 files changed

+556/−337 lines, 21 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• ChangeLog.md +8/−0
• bench/Bench.hs +1/−1
• bench/Poisson.hs +3/−4
• mcmc.cabal +1/−1
• src/Mcmc/Acceptance.hs +78/−44
• src/Mcmc/Algorithm.hs +2/−1
• src/Mcmc/Algorithm/MC3.hs +15/−15
• src/Mcmc/Algorithm/MHG.hs +27/−21
• src/Mcmc/Chain/Chain.hs +3/−4
• src/Mcmc/Chain/Save.hs +4/−4
• src/Mcmc/Cycle.hs +62/−33
• src/Mcmc/Internal/SpecFunctions.hs +2/−2
• src/Mcmc/MarginalLikelihood.hs +2/−2
• src/Mcmc/Mcmc.hs +54/−22
• src/Mcmc/Proposal.hs +54/−46
• src/Mcmc/Proposal/Generic.hs +17/−17
• src/Mcmc/Proposal/Hamiltonian/Common.hs +2/−1
• src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs +9/−19
• src/Mcmc/Proposal/Hamiltonian/Internal.hs +83/−48
• src/Mcmc/Proposal/Hamiltonian/Masses.hs +110/−33
• src/Mcmc/Proposal/Hamiltonian/Nuts.hs +19/−19

@@ -5,6 +5,14 @@ ## Unreleased changes  +## 0.8.1.0++-   Automatic intermediate tuning for HMC and NUTS.+-   Fix documentation for generic proposals.+-   Tooling updates.+-   Slight runtime improvements (strictness annotations).++ ## 0.8.0.1  -   Improve exception handling (also during execution of monitors; also improve

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

@@ -1,7 +1,6 @@ {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE TypeSynonymInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# OPTIONS_GHC -Wno-orphans #-} @@ -60,10 +59,10 @@ lh x = product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]  proposalAlpha :: Proposal I-proposalAlpha = (VB.element 0) @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune+proposalAlpha = VB.element 0 @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune  proposalBeta :: Proposal I-proposalBeta = (VB.element 1) @~ slideSymmetric 0.2 (PName "Beta") (pWeight 1) NoTune+proposalBeta = VB.element 1 @~ slideSymmetric 0.2 (PName "Beta") (pWeight 1) NoTune  proposals :: Cycle I proposals = cycleFromList [proposalAlpha, proposalBeta]@@ -103,7 +102,7 @@   void $ mcmc s a  toVec :: I -> VS.Vector Double-toVec xs = VS.generate 2 (\i -> xs VB.! i)+toVec xs = VS.generate 2 (xs VB.!)  fromVec :: I -> VS.Vector Double -> I fromVec _ xs = VB.mk2 (xs VS.! 0) (xs VS.! 1)

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

@@ -15,11 +15,13 @@   ( -- * Acceptance rates     AcceptanceRate,     AcceptanceCounts (..),-    Acceptance (fromAcceptance),+    AcceptanceRates (..),+    Acceptance,+    Acceptances (fromAcceptances),     emptyA,     pushAccept,     pushReject,-    pushAcceptanceCounts,+    ResetAcceptance (..),     resetA,     transformKeysA,     acceptanceRate,@@ -44,48 +46,74 @@  $(deriveJSON defaultOptions ''AcceptanceCounts) -addAccept :: AcceptanceCounts -> AcceptanceCounts-addAccept (AcceptanceCounts a r) = AcceptanceCounts (a + 1) r+-- | Proposals based on Hamiltonian dynamics use expected acceptance rates, not counts.+data AcceptanceRates = AcceptanceRates+  { totalAcceptanceRate :: !Double,+    nAcceptanceRates :: !Int+  }+  deriving (Show, Eq) -addReject :: AcceptanceCounts -> AcceptanceCounts-addReject (AcceptanceCounts a r) = AcceptanceCounts a (r + 1)+$(deriveJSON defaultOptions ''AcceptanceRates) -addAcceptanceCounts :: AcceptanceCounts -> AcceptanceCounts -> AcceptanceCounts-addAcceptanceCounts (AcceptanceCounts al rl) (AcceptanceCounts ar rr) =-  AcceptanceCounts (al + ar) (rl + rr)+-- | Stored actual acceptance counts and maybe expected acceptance rates.+data Acceptance = A AcceptanceCounts (Maybe AcceptanceRates)+  deriving (Show, Eq) +$(deriveJSON defaultOptions ''Acceptance)++addAccept :: Maybe AcceptanceRates -> Acceptance -> Acceptance+addAccept mr' (A (AcceptanceCounts a r) mr) = A (AcceptanceCounts (a + 1) r) (addAcceptanceRates mr' mr)++addReject :: Maybe AcceptanceRates -> Acceptance -> Acceptance+addReject mr' (A (AcceptanceCounts a r) mr) = A (AcceptanceCounts a (r + 1)) (addAcceptanceRates mr' mr)++addAcceptanceRates :: Maybe AcceptanceRates -> Maybe AcceptanceRates -> Maybe AcceptanceRates+addAcceptanceRates Nothing Nothing = Nothing+addAcceptanceRates (Just r) Nothing = Just r+addAcceptanceRates Nothing (Just r) = Just r+addAcceptanceRates (Just (AcceptanceRates al rl)) (Just (AcceptanceRates ar rr)) =+  Just $ AcceptanceRates (al + ar) (rl + rr)+ -- | For each key @k@, store the number of accepted and rejected proposals.-newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k AcceptanceCounts}+newtype Acceptances k = Acceptances {fromAcceptances :: M.Map k Acceptance}   deriving (Eq, Show) -instance ToJSONKey k => ToJSON (Acceptance k) where-  toJSON (Acceptance m) = toJSON m-  toEncoding (Acceptance m) = toEncoding m+instance ToJSONKey k => ToJSON (Acceptances k) where+  toJSON (Acceptances m) = toJSON m+  toEncoding (Acceptances m) = toEncoding m -instance (Ord k, FromJSONKey k) => FromJSON (Acceptance k) where-  parseJSON v = Acceptance <$> parseJSON v+instance (Ord k, FromJSONKey k) => FromJSON (Acceptances k) where+  parseJSON v = Acceptances <$> parseJSON v  -- | In the beginning there was the Word. -- -- Initialize an empty storage of accepted/rejected values.-emptyA :: Ord k => [k] -> Acceptance k-emptyA ks = Acceptance $ M.fromList [(k, AcceptanceCounts 0 0) | k <- ks]+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 => k -> Acceptance k -> Acceptance k-pushAccept k = Acceptance . M.adjust addAccept k . fromAcceptance+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 => k -> Acceptance k -> Acceptance k-pushReject k = Acceptance . M.adjust addReject k . fromAcceptance+pushReject :: Ord k => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k+pushReject mr k = Acceptances . M.adjust (addReject mr) k . fromAcceptances --- | For key @k@, add acceptance counts.-pushAcceptanceCounts :: Ord k => k -> AcceptanceCounts -> Acceptance k -> Acceptance k-pushAcceptanceCounts k c = Acceptance . M.adjust (addAcceptanceCounts c) k . fromAcceptance+-- | Reset acceptance specification.+data ResetAcceptance+  = -- | Reset actual acceptance counts and expected acceptance rates.+    ResetEverything+  | -- | Only reset expected acceptance rates.+    ResetExpectedRatesOnly  -- | Reset acceptance counts.-resetA :: Ord k => Acceptance k -> Acceptance k-resetA = emptyA . M.keys . fromAcceptance+resetA :: Ord k => ResetAcceptance -> Acceptances k -> Acceptances k+resetA ResetEverything = emptyA . M.keys . fromAcceptances+resetA ResetExpectedRatesOnly = Acceptances . M.map f . fromAcceptances+  where+    f (A cs _) = A cs Nothing  transformKeys :: (Ord k1, Ord k2) => [(k1, k2)] -> M.Map k1 v -> M.Map k2 v transformKeys ks m = foldl' insrt M.empty ks@@ -94,10 +122,11 @@  -- | Transform keys using the given lists. Keys not provided will not be present -- in the new 'Acceptance' variable.-transformKeysA :: (Ord k1, Ord k2) => [(k1, k2)] -> Acceptance k1 -> Acceptance k2-transformKeysA ks = Acceptance . transformKeys ks . fromAcceptance+transformKeysA :: (Ord k1, Ord k2) => [(k1, k2)] -> Acceptances k1 -> Acceptances k2+transformKeysA ks = Acceptances . transformKeys ks . fromAcceptances --- | Acceptance counts and rate for a specific proposal.+-- | Compute acceptance counts, and actual and expected acceptances rates for a+-- specific proposal. -- -- Return @Just (accepts, rejects, acceptance rate)@. --@@ -106,23 +135,28 @@ acceptanceRate ::   Ord k =>   k ->-  Acceptance k ->-  Maybe (Int, Int, AcceptanceRate)-acceptanceRate k a = case fromAcceptance a M.!? k of-  Just (AcceptanceCounts 0 0) -> Nothing-  Just (AcceptanceCounts as rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))+  Acceptances k ->+  -- | (nAccepts, nRejects, actualRate, expectedRate)+  (Int, Int, Maybe AcceptanceRate, Maybe AcceptanceRate)+acceptanceRate k a = case fromAcceptances a M.!? k of+  Just (A (AcceptanceCounts as rs) mrs) -> (as, rs, mar, mtr)+    where+      s = as + rs+      mar = if s <= 0 then Nothing else Just $ fromIntegral as / fromIntegral s+      mtr = case mrs of+        Nothing -> Nothing+        Just (AcceptanceRates xs n) -> Just $ xs / fromIntegral n   Nothing -> error "acceptanceRate: Key not found in map." --- | Acceptance rates for all proposals.+-- | Compute actual acceptance rates for all proposals. -- -- Set rate to 'Nothing' if no proposals have been accepted or rejected -- (division by zero).-acceptanceRates :: Acceptance k -> M.Map k (Maybe AcceptanceRate)-acceptanceRates =-  M.map-    ( \(AcceptanceCounts as rs) ->-        if as + rs == 0-          then Nothing-          else Just $ fromIntegral as / fromIntegral (as + rs)-    )-    . fromAcceptance+acceptanceRates :: Acceptances k -> M.Map k (Maybe AcceptanceRate)+acceptanceRates = M.map getRate . fromAcceptances+  where+    getRate (A (AcceptanceCounts as rs) _) =+      let s = as + rs+       in if s <= 0+            then Nothing+            else Just $ fromIntegral as / fromIntegral s

@@ -16,6 +16,7 @@  import qualified Data.ByteString.Lazy.Char8 as BL import Data.Time+import Mcmc.Acceptance import Mcmc.Cycle import Mcmc.Proposal import Mcmc.Settings@@ -42,7 +43,7 @@   aAutoTune :: TuningType -> Int -> a -> IO a    -- | Reset acceptance counts.-  aResetAcceptance :: a -> a+  aResetAcceptance :: ResetAcceptance -> a -> a    -- | Clean after burn in. In particular, this is used to reduce the length of   -- the trace, if required.

@@ -127,7 +127,7 @@     savedMC3Chains :: V.Vector (SavedChain a),     savedMC3ReciprocalTemperatures :: ReciprocalTemperatures,     savedMC3Iteration :: Int,-    savedMC3SwapAcceptance :: Acceptance Int,+    savedMC3SwapAcceptance :: Acceptances Int,     savedMC3Generator :: (Word64, Word64)   }   deriving (Eq, Show)@@ -172,7 +172,7 @@     -- | Current iteration.     mc3Iteration :: Int,     -- | Number of accepted and rejected swaps.-    mc3SwapAcceptance :: Acceptance Int,+    mc3SwapAcceptances :: Acceptances Int,     mc3Generator :: IOGenM StdGen   } @@ -252,7 +252,7 @@ initMHG prf lhf i beta a   | i < 0 = error "initMHG: Chain index negative."   -- Do nothing for the cold chain.-  | i == 0 = return $ MHG $ c+  | i == 0 = return $ MHG c   | otherwise = do       -- We have to push the current link in the trace, since it is not set by       -- 'setReciprocalTemperature'. The other links in the trace are still@@ -411,11 +411,11 @@       -- traceIO $ "Log priors (left, right, after swap): " <> show (ln prL') <> " " <> show (ln prR')       -- traceIO $ "Log likelihoods (left, right, before swap): " <> show (ln lhL) <> " " <> show (ln lhR)       -- traceIO $ "Log likelihood (left, right, after swap): " <> show (ln lhL') <> " " <> show (ln lhR')-      let !ac' = pushAccept i (mc3SwapAcceptance a)-      return $ a {mc3MHGChains = y, mc3SwapAcceptance = ac'}+      let !ac' = pushAccept Nothing i (mc3SwapAcceptances a)+      return $ a {mc3MHGChains = y, mc3SwapAcceptances = ac'}     else do-      let !ac' = pushReject i (mc3SwapAcceptance a)-      return $ a {mc3SwapAcceptance = ac'}+      let !ac' = pushReject Nothing i (mc3SwapAcceptances a)+      return $ a {mc3SwapAcceptances = ac'}   where     g = mc3Generator a @@ -483,7 +483,7 @@   mhgs' <- V.mapM (aAutoTune b l) $ mc3MHGChains a   -- 2. Auto tune temperatures.   let optimalRate = getOptimalRate PDimensionUnknown-      mCurrentRates = acceptanceRates $ mc3SwapAcceptance a+      mCurrentRates = acceptanceRates $ mc3SwapAcceptances a       -- We assume that the acceptance rate of state swaps between two chains is       -- roughly proportional to the ratio of the temperatures of the chains.       -- Hence, we focus on temperature ratios, actually reciprocal temperature@@ -512,15 +512,15 @@             (V.tail mhgs')   return $ a {mc3MHGChains = mhgs'', mc3ReciprocalTemperatures = bs'} -mc3ResetAcceptance :: ToJSON a => MC3 a -> MC3 a-mc3ResetAcceptance a = a'+mc3ResetAcceptance :: ToJSON a => ResetAcceptance -> MC3 a -> MC3 a+mc3ResetAcceptance x a = a'   where     -- 1. Reset acceptance of all chains.-    mhgs' = V.map aResetAcceptance (mc3MHGChains a)+    mhgs' = V.map (aResetAcceptance x) (mc3MHGChains a)     -- 2. Reset acceptance of swaps.-    ac' = resetA $ mc3SwapAcceptance a+    ac' = resetA x $ mc3SwapAcceptances a     ---    a' = a {mc3MHGChains = mhgs', mc3SwapAcceptance = ac'}+    a' = a {mc3MHGChains = mhgs', mc3SwapAcceptances = ac'}  mc3CleanAfterBurnIn :: ToJSON a => TraceLength -> MC3 a -> IO (MC3 a) mc3CleanAfterBurnIn tl a = do@@ -573,14 +573,14 @@     -- Acceptance rates may be 'Nothing' when no proposals have been undertaken.     -- The 'sequence' operations pull the 'Nothing's out of the inner     -- structures.-    as = sequence $ V.map (sequence . acceptanceRates . acceptance) cs+    as = sequence $ V.map (sequence . acceptanceRates . acceptances) cs     mVecAr = V.map (\mp -> sum mp / fromIntegral (length mp)) <$> as     mAr = (\vec -> V.sum vec / fromIntegral (V.length vec)) <$> mVecAr     bs = mc3ReciprocalTemperatures a     bsB = map (BB.toLazyByteString . BB.formatDouble (BB.standard 2)) $ U.toList bs     swapPeriod = fromSwapPeriod $ mc3SwapPeriod $ mc3Settings a     swapPeriodB = BB.toLazyByteString $ BB.intDec swapPeriod-    swapAcceptance = mc3SwapAcceptance a+    swapAcceptance = mc3SwapAcceptances a     n = fromNChains $ mc3NChains $ mc3Settings a     proposalHLine = BL.replicate (BL.length proposalHeader) '-' 

@@ -233,7 +233,7 @@ mhgPropose :: MHG a -> Proposal a -> IO (MHG a) mhgPropose (MHG c) p = do   -- 1. Sample new state.-  !(pres, mcs) <- liftIO $ s x g+  (!pres, !mcs) <- liftIO $ s x g   -- 2. Define new prior and likelihood calculation functions. Avoid actual   -- calculation of the values.   --@@ -242,15 +242,11 @@   -- https://stackoverflow.com/a/46603680/3536806.   let calcPrLh y = (pF y, lF y) `using` parTuple2 rdeepseq rdeepseq       accept y pr lh =-        let !ac' = case mcs of-              Nothing -> pushAccept p ac-              Just cs -> pushAcceptanceCounts p cs ac-         in pure $ MHG $ c {link = Link y pr lh, acceptance = ac'}+        let !ac' = pushAccept mcs p ac+         in pure $ MHG $ c {link = Link y pr lh, acceptances = ac'}       reject =-        let !ac' = case mcs of-              Nothing -> pushReject p ac-              Just cs -> pushAcceptanceCounts p cs ac-         in pure $ MHG $ c {acceptance = ac'}+        let !ac' = pushReject mcs p ac+         in pure $ MHG $ c {acceptances = ac'}   -- 3. Accept or reject.   --   -- 3a. When rejection is inevitable, avoid calculation of the prior, the@@ -274,7 +270,7 @@     (Link x pX lX) = link c     pF = priorFunction c     lF = likelihoodFunction c-    ac = acceptance c+    ac = acceptances c     g = generator c  mhgPush :: MHG a -> IO (MHG a)@@ -314,21 +310,31 @@     g = generator c  mhgAutoTune :: TuningType -> Int -> MHG a -> IO (MHG a)-mhgAutoTune b n (MHG c) = do-  mxs <--    if ccRequireTrace cc-      then Just . VB.map state <$> takeT n tr-      else pure Nothing-  return $ MHG $ c {cycle = autoTuneCycle b ac mxs cc}+mhgAutoTune tt n (MHG c)+  | isIntermediate =+      pure . MHG $+        if ccHasIntermediateTuners cc+          then -- Do not provide trace when tuning intermediately.+            c {cycle = autoTuneCycle tt ac Nothing cc}+          else -- Skip intermediate tuning completely when unnecessary.+            c+  | otherwise = do+      mxs <-+        -- Provide the trace if required.+        if ccRequireTrace cc+          then Just . VB.map state <$> takeT n tr+          else pure Nothing+      pure $ MHG c {cycle = autoTuneCycle tt ac mxs cc}   where-    ac = acceptance c+    isIntermediate = tt == IntermediateTuningFastProposalsOnly || tt == IntermediateTuningAllProposals+    ac = acceptances c     cc = cycle c     tr = trace c -mhgResetAcceptance :: MHG a -> MHG a-mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac}+mhgResetAcceptance :: ResetAcceptance -> MHG a -> MHG a+mhgResetAcceptance a (MHG c) = MHG $ c {acceptances = resetA a ac}   where-    ac = acceptance c+    ac = acceptances c  mhgCleanAfterBurnIn :: TraceLength -> MHG a -> IO (MHG a) mhgCleanAfterBurnIn tl (MHG c) = do@@ -346,7 +352,7 @@ mhgSummarizeCycle m (MHG c) = summarizeCycle m ac cc   where     cc = cycle c-    ac = acceptance c+    ac = acceptances c  mhgOpenMonitors ::   AnalysisName ->

@@ -92,10 +92,9 @@     -- | The 'Trace' of the Markov chain. In contrast to the link, the trace is     -- updated only after all proposals in the cycle have been executed.     trace :: Trace a,-    -- | For each 'Proposal', store the list of accepted (True) and rejected (False)-    -- proposals; for reasons of efficiency, the list is also stored in reverse-    -- order.-    acceptance :: Acceptance (Proposal a),+    -- | For each 'Proposal', store actual acceptance counts and for some+    -- proposals also expected acceptance rates.+    acceptances :: Acceptances (Proposal a),     -- | The random number generator.     generator :: IOGenM StdGen,     --

@@ -51,7 +51,7 @@   { savedLink :: Link a,     savedIteration :: Int,     savedTrace :: C.Stack VB.Vector (Link a),-    savedAcceptance :: Acceptance Int,+    savedAcceptances :: Acceptances Int,     savedSeed :: (Word64, Word64),     savedTuningParameters :: [Maybe (TuningParameter, AuxiliaryTuningParameters)]   }@@ -109,18 +109,18 @@                 "fromSave: Given likelihood:" <> show (lh $ state it) <> "."               ]        in error msg-  | length (fromAcceptance ac) /= length (ccProposals cc) =+  | length (fromAcceptances ac) /= length (ccProposals cc) =       let msg =             unlines               [ "fromSave: The number of proposals does not match.",-                "fromSave: Number of saved proposals:" <> show (length $ fromAcceptance ac) <> ".",+                "fromSave: Number of saved proposals:" <> show (length $ fromAcceptances ac) <> ".",                 "fromSave: Number of given proposals:" <> show (length $ ccProposals cc) <> "."               ]        in error msg   | otherwise = fromSavedChainUnsafe pr lh cc mn sv   where     it = savedLink sv-    ac = savedAcceptance sv+    ac = savedAcceptances sv  -- | See 'fromSavedChain' but do not perform sanity checks. Useful when -- restarting a run with changed prior function, likelihood function or

@@ -15,7 +15,7 @@ module Mcmc.Cycle   ( -- * Cycles     Order (..),-    Cycle (ccProposals, ccRequireTrace),+    Cycle (ccProposals, ccRequireTrace, ccHasIntermediateTuners),     cycleFromList,     setOrder,     IterationMode (..),@@ -27,10 +27,12 @@   ) where +import Control.Applicative import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL import Data.List import qualified Data.Map.Strict as M+import Data.Maybe import qualified Data.Vector as VB import Mcmc.Acceptance import Mcmc.Internal.Shuffle@@ -91,7 +93,10 @@   { ccProposals :: [Proposal a],     ccOrder :: Order,     -- | Does the cycle require the trace when auto tuning? See 'tRequireTrace'.-    ccRequireTrace :: Bool+    ccRequireTrace :: Bool,+    -- | Does the cycle include proposals that can be tuned every iterations?+    -- See 'tSuitableForIntermediateTuning'.+    ccHasIntermediateTuners :: Bool   }  -- | Create a 'Cycle' from a list of 'Proposal's; use 'RandomO', but see 'setOrder'.@@ -100,7 +105,7 @@   error "cycleFromList: Received an empty list but cannot create an empty Cycle." cycleFromList xs =   if length uniqueXs == length xs-    then Cycle xs RandomO (any needsTrace xs)+    then Cycle xs RandomO (any needsTrace xs) (any isIntermediate xs)     else error $ "\n" ++ msg ++ "cycleFromList: Proposals are not unique."   where     uniqueXs = nub xs@@ -109,9 +114,8 @@     removedDescriptions = map (show . prDescription) removedXs     removedMsgs = zipWith (\n d -> n ++ " " ++ d) removedNames removedDescriptions     msg = unlines removedMsgs-    needsTrace p = case prTuner p of-      Nothing -> False-      Just t -> tRequireTrace t+    needsTrace p = maybe False tRequireTrace (prTuner p)+    isIntermediate p = maybe False tSuitableForIntermediateTuning (prTuner p)  -- | Set the order of 'Proposal's in a 'Cycle'. setOrder :: Order -> Cycle a -> Cycle a@@ -123,9 +127,13 @@  -- | Replicate 'Proposal's according to their weights and possibly shuffle them. prepareProposals :: StatefulGen g m => IterationMode -> Cycle a -> g -> m [Proposal a]-prepareProposals m (Cycle xs o _) g =+prepareProposals m (Cycle xs o _ _) g =   if null ps-    then error "prepareProposals: No proposals found."+    then+      let msg = case m of+            FastProposals -> "no fast proposals found"+            AllProposals -> "no proposals found"+       in error $ "prepareProposals: " <> msg     else case o of       RandomO -> shuffle ps g       SequentialO -> return ps@@ -146,7 +154,7 @@  -- The number of proposals depends on the order. getNProposalsPerCycle :: IterationMode -> Cycle a -> Int-getNProposalsPerCycle m (Cycle xs o _) = case o of+getNProposalsPerCycle m (Cycle xs o _ _) = case o of   RandomO -> once   SequentialO -> once   RandomReversibleO -> 2 * once@@ -160,40 +168,61 @@ -- See 'tuneWithTuningParameters' and 'Tuner'. tuneWithChainParameters ::   TuningType ->-  AcceptanceRate ->+  Maybe AcceptanceRate ->   Maybe (VB.Vector a) ->   Proposal a ->   Either String (Proposal a)-tuneWithChainParameters b ar mxs p = case prTuner p of+tuneWithChainParameters tt mar mxs p = case prTuner p of   Nothing -> Right p-  Just (Tuner t ts rt fT _) -> case (rt, mxs) of-    (True, Nothing) -> error "tuneWithChainParameters: trace required"-    _ ->-      let (t', ts') = fT b d ar mxs (t, ts)-       in tuneWithTuningParameters t' ts' p-      where-        d = prDimension p+  Just (Tuner t ts rt it fT _) -> case (tt, it, prSpeed p) of+    (IntermediateTuningFastProposalsOnly, True, PFast) -> tuneIntermediate+    (IntermediateTuningAllProposals, True, _) -> tuneIntermediate+    (NormalTuningFastProposalsOnly, _, PFast) -> tuneNormally+    (NormalTuningAllProposals, _, _) -> tuneNormally+    (LastTuningFastProposalsOnly, _, _) -> tuneNormally+    (LastTuningAllProposals, _, _) -> tuneNormally+    _ -> Right p+    where+      hasTrace = isJust mxs+      err m = Left $ "tuneWithChainParameters: " <> m+      tuneIntermediate =+        if hasTrace+          then err "intermediate tuning but trace provided"+          else tune+      tuneNormally =+        if rt && not hasTrace+          then err "trace required"+          else tune+      tune =+        let (t', ts') = fT tt (prDimension p) mar mxs (t, ts)+         in tuneWithTuningParameters t' ts' p +-- (_, False, Just _) ->+-- (True, _, Just _) ->+-- (False, True, Nothing) ->+-- _ ->+ -- | Calculate acceptance rates and auto tunes the 'Proposal's in the 'Cycle'. -- -- Do not change 'Proposal's that are not tuneable.-autoTuneCycle :: TuningType -> Acceptance (Proposal a) -> Maybe (VB.Vector a) -> Cycle a -> Cycle a-autoTuneCycle b a mxs c = case (ccRequireTrace c, mxs) of-  (False, Just _) -> error "autoTuneCycle: trace not required"-  (True, Nothing) -> error "autoTuneCycle: trace required"-  _ ->-    if sort (M.keys ar) == sort ps-      then c {ccProposals = map tuneF ps}-      else error "autoTuneCycle: Proposals in map and cycle do not match."-    where-      ar = acceptanceRates a-      ps = ccProposals c-      tuneF p = case ar M.!? p of-        Just (Just x) -> either error id $ tuneWithChainParameters b x mxs p-        _ -> p+autoTuneCycle :: TuningType -> Acceptances (Proposal a) -> Maybe (VB.Vector a) -> Cycle a -> Cycle a+autoTuneCycle tt a mxs c+  | isJust mxs && not (ccRequireTrace c) = err "trace provided but not required"+  | otherwise =+      if sort (M.keys $ fromAcceptances a) == sort ps+        then c {ccProposals = map tuneF ps}+        else err "proposals in map and cycle do not match"+  where+    err msg = error $ "autoTuneCycle: " <> msg+    ps = ccProposals c+    tuneF p =+      let (_, _, mar, mtr) = acceptanceRate p a+          -- Favor the expected rate, if available.+          mr = mtr <|> mar+       in either error id $ tuneWithChainParameters tt mr mxs p  -- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance rates.-summarizeCycle :: IterationMode -> Acceptance (Proposal a) -> Cycle a -> BL.ByteString+summarizeCycle :: IterationMode -> Acceptances (Proposal a) -> Cycle a -> BL.ByteString summarizeCycle m a c =   BL.intercalate "\n" $     [ "Summary of proposal(s) in cycle.",

@@ -35,7 +35,7 @@ -- 'Numeric.SpecFunctions.logGamma'. logGammaG :: (Typeable a, RealFloat a) => a -> a logGammaG z-  | typeOf z == typeOf (0 :: Double) = unsafeCoerce logGamma z+  | typeOf z == typeRep (Proxy :: Proxy Double) = unsafeCoerce logGamma z   | otherwise = logGammaNonDouble z {-# SPECIALIZE logGammaG :: Double -> Double #-} @@ -214,7 +214,7 @@ -- 'Numeric.SpecFunctions.logFactorial'. logFactorialG :: forall a b. (Integral a, RealFloat b, Typeable b) => a -> b logFactorialG n-  | typeOf (undefined :: b) == typeOf (0 :: Double) = unsafeCoerce $ logFactorial n+  | typeRep (Proxy :: Proxy b) == typeRep (Proxy :: Proxy Double) = unsafeCoerce $ logFactorial n   | otherwise = logFactorialNonDouble n {-# SPECIALIZE logFactorialG :: Int -> Double #-} 

@@ -147,12 +147,12 @@ sampleAtPoint x ss lhf a = do   a'' <- liftIO $ mcmc ss' a'   let ch'' = fromMHG a''-      ac = acceptance ch''+      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. The tuning period may be too small."+    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."

@@ -25,11 +25,13 @@ import Control.Monad import Control.Monad.IO.Class import Control.Monad.Trans.Reader+import Data.Functor+import Mcmc.Acceptance (ResetAcceptance (ResetEverything, ResetExpectedRatesOnly)) import Mcmc.Algorithm import Mcmc.Cycle import Mcmc.Environment import Mcmc.Logger-import Mcmc.Proposal (TuningType (LastTuningStep, NormalTuningStep))+import Mcmc.Proposal import Mcmc.Settings import System.IO import Prelude hiding (cycle)@@ -52,7 +54,7 @@ mcmcResetAcceptance :: Algorithm a => a -> MCMC a mcmcResetAcceptance a = do   logDebugB "Reset acceptance rates."-  return $ aResetAcceptance a+  return $ aResetAcceptance ResetEverything a  mcmcExceptionHandler :: Algorithm a => Environment Settings -> a -> AsyncException -> IO b mcmcExceptionHandler e a err = do@@ -77,8 +79,16 @@   mStdLog <- liftIO $ aExecuteMonitors vb t0 iTotal a   forM_ mStdLog (logOutB "   ") -mcmcIterate :: Algorithm a => IterationMode -> Int -> a -> MCMC a-mcmcIterate m n a+-- When intermediate tuning is activated, specific proposals get tuned every+-- iterations.+data IntermediateTuningSpec+  = IntermediateTuningFastProposalsOnlyOn+  | IntermediateTuningAllProposalsOn+  | 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@@ -87,7 +97,20 @@       -- NOTE: Handle interrupts during iterations, before writing monitors,       -- using the old algorithm state @a@.       let handlerOld = mcmcExceptionHandler e a-          actionIterate = aIterate m p 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'@.@@ -100,8 +123,8 @@       -- recover from partly written monitor files.       let handlerNew = mcmcExceptionHandler e a'           actionWrite = runReaderT (mcmcExecuteMonitors a') e-      liftIO $ (uninterruptibleMask_ actionWrite) `catch` handlerNew-      mcmcIterate m (n - 1) a'+      liftIO $ uninterruptibleMask_ actionWrite `catch` handlerNew+      mcmcIterate t m (n - 1) a'  mcmcNewRun :: Algorithm a => a -> MCMC a mcmcNewRun a = do@@ -112,15 +135,15 @@   mcmcExecuteMonitors a   when (aIsInvalidState a) (logWarnB "The initial state is invalid!")   a' <- mcmcBurnIn a-  logInfoS $ "Cleaning chain after burn in."+  logInfoS "Cleaning chain after burn in."   let tl = sTraceLength s   a'' <- liftIO $ aCleanAfterBurnIn tl a'-  logInfoS $ "Saving chain after burn in."+  logInfoS "Saving chain after burn in."   mcmcSave a''   let i = fromIterations $ sIterations s   logInfoS $ "Running chain for " ++ show i ++ " iterations."   logInfoB $ aStdMonitorHeader a''-  mcmcIterate AllProposals i a''+  mcmcIterate IntermediateTuningOff AllProposals i a''  mcmcContinueRun :: Algorithm a => a -> MCMC a mcmcContinueRun a = do@@ -139,7 +162,7 @@   logInfoB $ aSummarizeCycle AllProposals a   logInfoS $ "Running chain for " ++ show di ++ " iterations."   logInfoB $ aStdMonitorHeader a-  mcmcIterate AllProposals di a+  mcmcIterate IntermediateTuningOff AllProposals di a  mcmcBurnIn :: Algorithm a => a -> MCMC a mcmcBurnIn a = do@@ -154,7 +177,7 @@       logInfoS $ "Burning in for " <> show n <> " iterations."       logInfoS "Auto tuning is disabled."       logInfoB $ aStdMonitorHeader a-      a' <- mcmcIterate AllProposals n a+      a' <- mcmcIterate IntermediateTuningOff AllProposals n a       logInfoB $ aSummarizeCycle AllProposals a'       a'' <- mcmcResetAcceptance a'       logInfoB "Burn in finished."@@ -191,25 +214,34 @@  -- Auto tune the proposals. mcmcAutotune :: Algorithm a => TuningType -> Int -> a -> MCMC a-mcmcAutotune NormalTuningStep n a = do-  logDebugB "Intermediate auto tune."-  liftIO $ aAutoTune NormalTuningStep n a-mcmcAutotune LastTuningStep n a = do-  logDebugB "Last auto tune."-  liftIO $ aAutoTune LastTuningStep n a+mcmcAutotune t n a = do+  case t of+    NormalTuningFastProposalsOnly -> logDebugB "Normal auto tune; fast proposals only."+    IntermediateTuningFastProposalsOnly -> pure ()+    LastTuningFastProposalsOnly -> logDebugB "Last auto tune; fast proposals only."+    NormalTuningAllProposals -> logDebugB "Normal auto tune; all proposals."+    IntermediateTuningAllProposals -> pure ()+    LastTuningAllProposals -> logDebugB "Last auto tune; all proposals."+  liftIO $ aAutoTune t n a  mcmcBurnInWithAutoTuning :: Algorithm a => IterationMode -> [Int] -> a -> MCMC a mcmcBurnInWithAutoTuning _ [] _ = error "mcmcBurnInWithAutoTuning: Empty list." mcmcBurnInWithAutoTuning m [x] a = do   -- Last round.-  a' <- mcmcIterate m x a-  a'' <- mcmcAutotune LastTuningStep x a'+  let (tti, ttl) = case m of+        FastProposals -> (IntermediateTuningFastProposalsOnlyOn, LastTuningFastProposalsOnly)+        AllProposals -> (IntermediateTuningAllProposalsOn, LastTuningAllProposals)+  a' <- mcmcIterate tti m x a+  a'' <- mcmcAutotune ttl x a'   logInfoB $ aSummarizeCycle m a''   logInfoS $ "Acceptance rates calculated over the last " <> show x <> " iterations."   mcmcResetAcceptance a'' mcmcBurnInWithAutoTuning m (x : xs) a = do-  a' <- mcmcIterate m x a-  a'' <- mcmcAutotune NormalTuningStep x a'+  let (tti, ttn) = case m of+        FastProposals -> (IntermediateTuningFastProposalsOnlyOn, NormalTuningFastProposalsOnly)+        AllProposals -> (IntermediateTuningAllProposalsOn, NormalTuningAllProposals)+  a' <- mcmcIterate tti m x a+  a'' <- mcmcAutotune ttn x a'   logDebugB $ aSummarizeCycle m a''   logDebugS $ "Acceptance rates calculated over the last " <> show x <> " iterations."   logDebugB $ aStdMonitorHeader a''

@@ -38,8 +38,6 @@     TuningFunction,     AuxiliaryTuningParameters,     tuningFunction,-    tuningFunctionWithAux,-    tuningFunctionOnlyAux,     tuningParameterMin,     tuningParameterMax,     tuneWithTuningParameters,@@ -207,16 +205,14 @@     Propose !a !KernelRatio !Jacobian   deriving (Show, Eq) --- TODO @Dominik (high, feature): Proposals should be aware of: Is this Burn in, or not?- -- | Simple proposal function without tuning information. -- -- Instruction about randomly moving from the current state to a new state, -- given some source of randomness. ----- Maybe report acceptance counts internal to the proposal (e.g., used by+-- Maybe report acceptance rates internal to the proposal (e.g., used by -- proposals based on Hamiltonian dynamics).-type PFunction a = a -> IOGenM StdGen -> IO (PResult a, Maybe AcceptanceCounts)+type PFunction a = a -> IOGenM StdGen -> IO (PResult a, Maybe AcceptanceRates)  -- | Create a proposal with a single tuning parameter. --@@ -243,7 +239,7 @@   where     fT = tuningFunction     g t _ = Right $ f t-    tuner = Tuner 1.0 VU.empty False fT g+    tuner = Tuner 1.0 VU.empty False False fT g createProposal r f s d n w NoTune =   Proposal n r s d w (f 1.0) Nothing @@ -253,6 +249,8 @@     tAuxiliaryTuningParameters :: AuxiliaryTuningParameters,     -- | Does the tuner require the trace over the last tuning period?     tRequireTrace :: Bool,+    -- | Can the tuner be used for intermediate tuning (see 'TuningType')?+    tSuitableForIntermediateTuning :: Bool,     tTuningFunction :: TuningFunction a,     -- | Given the tuning parameter, and the auxiliary tuning parameters, get     -- the tuned propose function.@@ -275,16 +273,35 @@ -- expected acceptance rate; and vice versa. type TuningParameter = Double --- | The last tuning step may be special.-data TuningType = NormalTuningStep | LastTuningStep+-- | Tuning type. To distinguish between fast and slow proposals, see+-- 'Mcmc.Cycle.IterationMode'.+data TuningType+  = -- | Normal tuning step with fast proposals only.+    NormalTuningFastProposalsOnly+  | -- | Intermediate tuning step executed after each iteration with fast+    -- proposals only. Only suitable for proposals which can calculate expected+    -- acceptance rates.+    IntermediateTuningFastProposalsOnly+  | -- | The last tuning step with fast proposals only may be special.+    LastTuningFastProposalsOnly+  | -- | Normal tuning step of all proposals.+    NormalTuningAllProposals+  | -- | Intermediate tuning step of all proposals.+    IntermediateTuningAllProposals+  | -- | The last tuning step with all proposals.+    LastTuningAllProposals+  deriving (Eq)  -- | Compute new tuning parameters. type TuningFunction a =   TuningType ->   PDimension ->-  -- | Acceptance rate of last tuning period.-  AcceptanceRate ->-  -- | Trace of last tuning period. Only available when requested by proposal.+  -- | Acceptance rate of last tuning period. May not always be available+  -- because proposals may be skipped.+  Maybe AcceptanceRate ->+  -- | Trace of last tuning period. Not available for intermediate tuning' steps+  -- (see 'TuningType'), and only available for other tuning types when+  -- requested by proposal.   Maybe (VB.Vector a) ->   (TuningParameter, AuxiliaryTuningParameters) ->   (TuningParameter, AuxiliaryTuningParameters)@@ -301,27 +318,15 @@  -- | Default tuning function. ----- The default tuning function only uses the acceptance rate. In particular, it--- does not handle auxiliary tuning parameters and ignores the actual samples--- attained during the last tuning period.+-- The default tuning function only uses the actual acceptance rate. In+-- particular, it does not handle auxiliary tuning parameters, ignores+-- intermediate tuning steps, and ignores the actual samples attained during the+-- last tuning period. tuningFunction :: TuningFunction a-tuningFunction _ d r _ (!t, !ts) = bimap (tuningFunctionSimple d r) id (t, ts)---- | Also tune auxiliary tuning parameters.-tuningFunctionWithAux ::-  -- | Auxiliary tuning function.-  (TuningType -> VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters) ->-  TuningFunction a-tuningFunctionWithAux _ _ _ _ Nothing _ = error "tuningFunctionWithAux: empty trace"-tuningFunctionWithAux f b d r (Just xs) (!t, !ts) = bimap (tuningFunctionSimple d r) (f b xs) (t, ts)---- | Only tune auxiliary tuning parameters.-tuningFunctionOnlyAux ::-  -- | Auxiliary tuning function.-  (TuningType -> VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters) ->-  TuningFunction a-tuningFunctionOnlyAux _ _ _ _ Nothing _ = error "tuningFunctionOnlyAux: empty trace"-tuningFunctionOnlyAux f b _ _ (Just xs) (!t, !ts) = bimap id (f b xs) (t, ts)+tuningFunction IntermediateTuningFastProposalsOnly _ _ _ t = t+tuningFunction IntermediateTuningAllProposals _ _ _ t = t+tuningFunction _ _ Nothing _ t = t+tuningFunction _ d (Just r) _ (!t, !ts) = first (tuningFunctionSimple d r) (t, ts)  -- IDEA: Per proposal type tuning parameter boundaries. For example, a sliding -- proposal with a large tuning parameter is not a problem. But then, if the@@ -330,7 +335,7 @@  -- | Minimal tuning parameter; subject to change. tuningParameterMin :: TuningParameter-tuningParameterMin = 1e-5+tuningParameterMin = 1e-6  -- | Maximal tuning parameter; subject to change. tuningParameterMax :: TuningParameter@@ -358,14 +363,14 @@   Either String (Proposal a) tuneWithTuningParameters t ts p = case prTuner p of   Nothing -> Left "tuneWithTuningParameters: Proposal is not tunable."-  Just (Tuner _ _ nt fT g) ->+  Just (Tuner _ _ reqTr inTn fT g) ->     -- Ensure that the tuning parameter is strictly positive and well bounded.     let t' = max tuningParameterMin t         t'' = min tuningParameterMax t'         psE = g t'' ts      in case psE of           Left err -> Left $ "tune: " <> err-          Right ps -> Right $ p {prFunction = ps, prTuner = Just $ Tuner t'' ts nt fT g}+          Right ps -> Right $ p {prFunction = ps, prTuner = Just $ Tuner t'' ts reqTr inTn fT g}  -- | See 'PDimension'. getOptimalRate :: PDimension -> Double@@ -431,7 +436,7 @@  -- Lift tuner from one data type to another. liftTunerWith :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b-liftTunerWith jf l (Tuner p ps nt fP g) = Tuner p ps nt fP' g'+liftTunerWith jf l (Tuner p ps reqTr inTn fP g) = Tuner p ps reqTr inTn fP' g'   where     fP' b d r = fP b d r . fmap (VB.map (^. l))     g' x xs = liftPFunctionWith jf l <$> g x xs@@ -455,14 +460,15 @@   BL.ByteString ->   BL.ByteString ->   BL.ByteString-renderRow name ptype weight nAccept nReject acceptRate optimalRate tuneParam manualAdjustment = nm <> pt <> wt <> na <> nr <> ra <> ro <> tp <> mt+renderRow name ptype weight nAccept nReject acceptRateActual optimalRate tuneParam manualAdjustment =+  nm <> pt <> wt <> na <> nr <> ra <> ro <> tp <> mt   where     nm = alignLeft 30 name     pt = alignLeft 50 ptype     wt = alignRight 8 weight     na = alignRight 14 nAccept     nr = alignRight 14 nReject-    ra = alignRight 14 acceptRate+    ra = alignRight 14 acceptRateActual     ro = alignRight 14 optimalRate     tp = alignRight 20 tuneParam     mt = alignRight 30 manualAdjustment@@ -476,7 +482,7 @@     "Weight"     "Accepted"     "Rejected"-    "Rate"+    "Actual rate"     "Optimal rate"     "Tuning parameter"     "Consider manual adjustment"@@ -488,7 +494,7 @@   PWeight ->   Maybe TuningParameter ->   PDimension ->-  Maybe (Int, Int, Double) ->+  (Int, Int, Maybe Double, Maybe Double) ->   BL.ByteString summarizeProposal name description weight tuningParameter dimension ar =   renderRow@@ -497,18 +503,19 @@     weightStr     nAccept     nReject-    acceptRate+    acceptRateActual+    -- acceptRateExpected     optimalRate     tuneParamStr     manualAdjustmentStr   where     fN n = BB.formatDouble (BB.standard n)     weightStr = BB.toLazyByteString $ BB.intDec $ fromPWeight weight-    nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) ar-    nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) ar-    acceptRate = BB.toLazyByteString $ maybe "" (fN 2 . (^. _3)) ar+    nAccept = BB.toLazyByteString $ BB.intDec $ ar ^. _1+    nReject = BB.toLazyByteString $ BB.intDec $ ar ^. _2+    acceptRateActual = BB.toLazyByteString $ maybe "" (fN 2) (ar ^. _3)     optimalRate = BB.toLazyByteString $ fN 2 $ getOptimalRate dimension-    tuneParamStr = BB.toLazyByteString $ maybe "" (fN 4) tuningParameter+    tuneParamStr = BB.toLazyByteString $ maybe "" (fN 6) tuningParameter     checkRate rate       | rate < rateMin = Just "rate too low"       | rate > rateMax = Just "rate too high"@@ -518,7 +525,8 @@       | tp >= (0.9 * tuningParameterMax) = Just "tuning parameter too high"       | otherwise = Nothing     tps = checkTuningParam =<< tuningParameter-    ars = (checkRate . (^. _3)) =<< ar+    -- Use actual acceptance rate.+    ars = checkRate =<< (ar ^. _3)     manualAdjustmentStr =       let        in case (ars, tps) of

@@ -19,26 +19,25 @@ import Numeric.Log import Statistics.Distribution --- | Generic function to create proposals for continuous parameters (e.g.,--- 'Double').+-- | Generic function to create proposals using a continuous auxiliary variable+-- of type 'Double'. -- -- The procedure is as follows: Let \(\mathbb{X}\) be the state space and \(x\) -- be the current state. -- -- 1. Let \(D\) be a continuous probability distribution on \(\mathbb{D}\);---    sample an auxiliary variable \(epsilon \sim D\).+--    sample an auxiliary variable \(u \sim D\). ----- 2. Suppose \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). PFunction a---    new state \(x' = x \odot \epsilon\).+-- 2. Let \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). Propose a+--    new state \(x' = x \odot u\). ----- 3. If the proposal is unbiased, the Metropolis-Hastings-Green ratio can---    directly be calculated using the posterior function.+-- If the proposal is unbiased, the Metropolis-Hastings-Green ratio can directly+-- be calculated using the posterior function. ----- 4. However, if the proposal is biased: Suppose \(g : \mathbb{D} \to---    \mathbb{D}\) inverses the auxiliary variable \(\epsilon\) such that \(x =---    x' \odot g(\epsilon)\). Calculate the Metropolis-Hastings-Green ratio---    using the posterior function, \(g\), \(D\), \(\epsilon\), and possibly a---    Jacobian function.+-- However, if the proposal is biased: Let \(g : \mathbb{D} \to \mathbb{D}\);+-- \(g\) inverses the auxiliary variable \(u\) such that \(x = x' \odot g(u)\).+-- Calculate the Metropolis-Hastings-Green ratio using the posterior function,+-- \(g\), \(D\), \(u\), and possibly a Jacobian function. genericContinuous ::   (ContDistr d, ContGen d) =>   -- | Probability distribution@@ -46,12 +45,12 @@   -- | Forward operator \(\odot\).   --   -- For example, for a multiplicative proposal on one variable the forward-  -- operator is @(*)@, so that @x * u = y@.+  -- operator is @(*)@, so that \(x' = x * u\).   (a -> Double -> a) ->   -- | Inverse operator \(g\) of the auxiliary variable.   --   -- For example, 'recip' for a multiplicative proposal on one variable, since-  -- @y * (recip u) = x * u * (recip u) = x@.+  -- \(x' * u^{-1} = x * u * u^{-1} = x\).   --   -- Required for biased proposals.   Maybe (Double -> Double) ->@@ -86,7 +85,8 @@   pure (Propose (x `f` u) r j, Nothing) {-# INLINEABLE genericContinuous #-} --- | Generic function to create proposals for discrete parameters (e.g., 'Int').+-- | Generic function to create proposals using a discrete auxiliary variable of+-- type 'Int'. -- -- See 'genericContinuous'. genericDiscrete ::@@ -95,11 +95,11 @@   d ->   -- | Forward operator.   ---  -- For example, (+), so that x + dx = x'.+  -- For example, (+), so that \(x + dx = x'\).   (a -> Int -> a) ->   -- | Inverse operator \(g\) of the auxiliary variable.   ---  -- For example, 'negate', so that x' + (negate dx) = x.+  -- For example, 'negate', so that \(x' - dx = x + dx - dx = x\).   --   -- Only required for biased proposals.   Maybe (Int -> Int) ->

@@ -59,7 +59,8 @@ import qualified Numeric.LinearAlgebra as L  -- NOTE: Implementing the Riemannian adaptation (state-dependent mass matrix).--- seems a little bit of an overkill.+-- seems a little bit of an overkill. See also+-- https://discourse.mc-stan.org/t/riemann-manifold-hmc-in-stan/19466/5.  -- | The Hamiltonian proposal acts on a vector of floating point values referred -- to as positions.

@@ -108,14 +108,6 @@   hParamsI <- fromAuxiliaryTuningParameters d ts   pure $ hamiltonianPFunction hParamsI hstruct targetWith --- TODO @Dominik (high, issue): Acceptance counts. How to combine with values--- reported here and from the NUTS sampler.---- TODO @Dominik (high, feature): The expected acceptance counts should not be--- calculated after burn in. Rather, the actual acceptance counts should be--- reported. For this to work, the proposal needs to know if it is in "burn in--- phase" or not.- -- The inverted covariance matrix and the log determinant of the covariance -- matrix are calculated by 'hamiltonianPFunction'. hamiltonianPFunction ::@@ -130,11 +122,13 @@   -- of epsilon. I still think it should vary because otherwise, there may be   -- dragons due to periodicity.   let lM = la / eRan-      lL = maximum [1 :: Int, floor $ 0.9 * lM]-      lR = maximum [lL, ceiling $ 1.1 * lM]+      lL = max (1 :: Int) (floor $ 0.9 * lM)+      lR = max lL (ceiling $ 1.1 * lM)   lRan <- uniformRM (lL, lR) g   case leapfrog (targetWith x) msI lRan eRan q p of-    Nothing -> pure (ForceReject, Just $ AcceptanceCounts 0 100)+    -- NOTE: I am not sure if it is correct to set the expected acceptance rate+    -- to 0 when the leapfrog integrator fails.+    Nothing -> pure (ForceReject, Just $ AcceptanceRates 0 1)     -- Check if next state is accepted here, because the Jacobian is included in     -- the target function. If not: pure (x, 0.0, 1.0).     Just (q', p', prQ, prQ') -> do@@ -149,13 +143,9 @@       -- chain back to the previous state. However, we are only interested in       -- the positions, and are not even storing the momenta.       let pr = if accept then ForceAccept (fromVec x q') else ForceReject-          ar = exp $ ln r-          getCounts s = max 0 $ min 100 $ round $ s * 100-          ac =-            if ar >= 0-              then let cs = getCounts ar in AcceptanceCounts cs (100 - cs)-              else error $ "hamiltonianPFunction: Acceptance rate negative."-      pure (pr, Just ac)+          -- Limit expected acceptance rate between 0 and 1.+          ar = max 0 $ min 1 (exp $ ln r)+      pure (pr, Just $ AcceptanceRates ar 1)   where     (HParamsI e la ms _ _ msI mus) = hparamsi     (HStructure _ toVec fromVec) = hstruct@@ -204,7 +194,7 @@       tuner = do         tfun <- hTuningFunctionWith dim toVec htconf         let pfun = hamiltonianPFunctionWithTuningParameters dim hstruct targetWith-        pure $ Tuner 1.0 ts True tfun pfun+        pure $ Tuner 1.0 ts True True tfun pfun    in case checkHStructureWith (hpsMasses hParamsI) hstruct of         Just err -> error err         Nothing -> hamiltonianWith tuner

@@ -55,6 +55,7 @@ import Control.Monad import Control.Monad.ST import Data.Foldable+import Data.Maybe import qualified Data.Vector.Storable as VS import qualified Data.Vector.Unboxed as VU import Mcmc.Proposal@@ -92,31 +93,48 @@   }   deriving (Show) --- The default tuning parameters in [4] are:+-- The default tuning parameters in [4] which have been tweaked for tuning the+-- proposals after every iteration are: -- --   mu = log $ 10 * eps --   ga = 0.05 --   t0 = 10 --   ka = 0.75 ----- However, these default tuning parameters will be off, because the authors--- suggesting these values tune the proposal after every single iteration.+-- For reference, I used the following default parameters with longer auto+-- tuning intervals. ----- The following values are tweaked for our case, where tuning does not happen--- after each iteration. Of course, we could tune the leapfrog parameters after--- each generation. Even the mass parameters could be tuned each iteration when--- the masses are estimated from more past iterations spanning many tuning--- intervals.+--   mu = log $ 10 * eps+--   ga = 0.1+--   t0 = 3+--   ka = 0.5 ----- NOTE: In theory, these we could expose these internal tuning parameters to--- the user.+-- Another good resource:+-- https://mc-stan.org/docs/2_29/reference-manual/hmc-algorithm-parameters.html.+--+-- NOTE: In theory, we could expose these internal tuning parameters to the+-- user. tParamsFixedWith :: LeapfrogScalingFactor -> TParamsFixed tParamsFixedWith eps = TParamsFixed eps mu ga t0 ka   where+    -- "Mu is a freely chosen point that the iterators are shrunk towards". I am+    -- not exactly sure what this means. The parameter does not seem to have+    -- much of an effect.     mu = log $ 10 * eps-    ga = 0.1-    t0 = 3-    ka = 0.5+    -- Gamma "controls the amount of shrinkage towards mu". The larger gamma is,+    -- the less variant epsilon is.+    --+    -- I changed this parameter from 0.05 to get better results in test runs.+    ga = 0.15+    -- "Free parameter that stabilizes the initial iterations". The larger t0+    -- is, the stabler epsilon is in the first iterations.+    t0 = 10+    -- "Setting the parameter ka < 1 allows us to give higher weight to more+    -- recent iterates and to more quickly forget the iterates produced during+    -- the early warmup stages."+    --+    -- I changed this parameter from 0.75 to get better results in test runs.+    ka = 0.75  -- All internal parameters. data HParamsI = HParamsI@@ -235,7 +253,7 @@           a :: Double           a = if rI > 0.5 then 1 else (-1)           go e r =-            if r ** a > 2 ** (negate a)+            if r ** a > 2 ** negate a               then case leapfrog t msI 1 e q p of                 Nothing -> e                 Just (_, p'', _, prQ'') ->@@ -257,38 +275,55 @@   (a -> Positions) ->   HTuningConf ->   Maybe (TuningFunction a)-hTuningFunctionWith n toVec (HTuningConf lc mc) = case (lc, mc) of-  (HNoTuneLeapfrog, HNoTuneMasses) -> Nothing-  (_, _) -> Just $-    \tt pdim ar mxs (_, !ts) ->-      case mxs of-        Nothing -> error "hTuningFunctionWith: empty trace"-        Just xs ->-          let (HParamsI eps la ms tpv tpf msI mus) =-                -- NOTE: Use error here, because a dimension mismatch is a serious bug.-                either error id $ fromAuxiliaryTuningParameters n ts-              (TParamsVar epsMean h m) = tpv-              (TParamsFixed eps0 mu ga t0 ka) = tpf-              (ms', msI') = case mc of-                HNoTuneMasses -> (ms, msI)-                HTuneDiagonalMassesOnly -> tuneDiagonalMassesOnly toVec xs (ms, msI)-                HTuneAllMasses -> tuneAllMasses toVec xs (ms, msI)-              (eps'', epsMean'', h'') = case lc of-                HNoTuneLeapfrog -> (eps, epsMean, h)-                HTuneLeapfrog ->-                  let delta = getOptimalRate pdim-                      c = recip $ m + t0-                      h' = (1.0 - c) * h + c * (delta - ar)-                      logEps' = mu - (sqrt m / ga) * h'-                      eps' = exp logEps'-                      mMKa = m ** (negate ka)-                      epsMean' = exp $ mMKa * logEps' + (1 - mMKa) * log epsMean-                   in (eps', epsMean', h')-              eps''' = case tt of-                NormalTuningStep -> eps''-                LastTuningStep -> epsMean''-              tpv' = TParamsVar epsMean'' h'' (m + 1.0)-           in (eps''' / eps0, toAuxiliaryTuningParameters $ HParamsI eps''' la ms' tpv' tpf msI' mus)+hTuningFunctionWith _ _ (HTuningConf HNoTuneLeapfrog HNoTuneMasses) = Nothing+hTuningFunctionWith n toVec (HTuningConf lc mc) = Just $ \tt pdim mar mxs (_, !ts) ->+  case tt of+    IntermediateTuningFastProposalsOnly -> err "fast intermediate tuning step but slow proposal"+    NormalTuningFastProposalsOnly -> err "fast normal tuning step but slow proposal"+    _ ->+      let (HParamsI eps la ms tpv tpf msI mus) =+            -- NOTE: Use error here, because a dimension mismatch is a serious bug.+            either error id $ fromAuxiliaryTuningParameters n ts+          (TParamsVar epsMean h m) = tpv+          (TParamsFixed eps0 mu ga t0 ka) = tpf+          m' = SmoothingParameter $ round m+          (ms', msI') = case tt of+            IntermediateTuningAllProposals -> (ms, msI)+            _ ->+              let xs = fromMaybe (err "empty trace") mxs+               in case mc of+                    HNoTuneMasses -> (ms, msI)+                    HTuneDiagonalMassesOnly -> tuneDiagonalMassesOnly m' toVec xs (ms, msI)+                    HTuneAllMasses -> tuneAllMasses m' toVec xs (ms, msI)+          (eps'', epsMean'', h'') = case tt of+            LastTuningFastProposalsOnly -> (eps, epsMean, h)+            _ -> case lc of+              HNoTuneLeapfrog -> (eps, epsMean, h)+              HTuneLeapfrog ->+                let ar = fromMaybe (err "no acceptance rate") mar+                    delta = getOptimalRate pdim+                    -- Algorithm 6; explained in Section 3.2.+                    --+                    -- Another good resource is the Tensorflow API+                    -- documentation:+                    -- https://www.tensorflow.org/probability/api_docs/python/tfp/mcmc/DualAveragingStepSizeAdaptation.+                    --+                    -- See also Nesterov (2007) Primal-dual subgradient methods+                    -- for convex problems, Mathematical Programming.+                    c = recip $ m + t0+                    h' = (1.0 - c) * h + c * (delta - ar)+                    eps' = exp $ mu - (sqrt m / ga) * h'+                    mMKa = m ** negate ka+                    -- Original formula is:+                    -- epsMean' = exp $ mMKa * logEps' + (1 - mMKa) * log epsMean+                    -- Which is the same as:+                    epsMean' = (eps' ** mMKa) * (epsMean ** (1 - mMKa))+                    epsF = if tt == LastTuningAllProposals then epsMean' else eps'+                 in (epsF, epsMean', h')+          tpv' = TParamsVar epsMean'' h'' (m + 1.0)+       in (eps'' / eps0, toAuxiliaryTuningParameters $ HParamsI eps'' la ms' tpv' tpf msI' mus)+  where+    err msg = error $ "hTuningFunctionWith: " <> msg  checkHStructureWith :: Foldable s => Masses -> HStructure s -> Maybe String checkHStructureWith ms (HStructure x toVec fromVec)@@ -361,7 +396,7 @@           else Nothing   -- L-1 full steps for positions and momenta. This gives the positions q_{L-1},   -- and the momenta p_{L-1/2}.-  (qLM1, pLM1Half) <- go (l - 1) $ Just $ (q, pHalf)+  (qLM1, pLM1Half) <- go (l - 1) $ Just (q, pHalf)   -- The last full step of the positions.   let qL = leapfrogStepPositions msI eps qLM1 pLM1Half   -- The last half step of the momenta.@@ -379,7 +414,7 @@           let qs' = leapfrogStepPositions msI eps qs ps               (x, ps') = leapfrogStepMomenta eps tF qs' p            in if x > 0.0-                then go (n - 1) $ Just $ (qs', ps')+                then go (n - 1) $ Just (qs', ps')                 else Nothing  leapfrogStepMomenta ::

@@ -21,6 +21,7 @@     massesToVector,      -- * Tuning+    SmoothingParameter (..),     tuneDiagonalMassesOnly,     tuneAllMasses,   )@@ -32,6 +33,7 @@ import qualified Data.Vector.Unboxed as VU import Mcmc.Proposal.Hamiltonian.Common import qualified Numeric.LinearAlgebra as L+import Numeric.Natural import qualified Statistics.Covariance as S import qualified Statistics.Function as S import qualified Statistics.Sample as S@@ -108,7 +110,7 @@ -- crucial. The Hamiltonian algorithms also work when the masses are off. cleanMatrix :: L.Matrix Double -> L.Matrix Double cleanMatrix xs =-  (L.diag $ L.cmap cleanDiag xsDiag) + L.cmap cleanOffDiag xsOffDiag+  L.diag (L.cmap cleanDiag xsDiag) + L.cmap cleanOffDiag xsOffDiag   where     xsDiag = L.takeDiag xs     cleanDiag x@@ -129,9 +131,9 @@   | n == 0 || m == 0 = error "getMassesI: Matrix empty."   | n /= m = error "getMassesI: Matrix not square."   | sign /= 1.0 = error "getMassesI: Determinant of matrix is negative."-  | otherwise = toGMatrix $ cleanMatrix $ xsI+  | otherwise = toGMatrix $ cleanMatrix xsI   where-    xs' = L.unSym $ xs+    xs' = L.unSym xs     n = L.rows xs'     m = L.cols xs'     (xsI, (_, sign)) = L.invlndet xs'@@ -179,20 +181,68 @@ getSampleSize :: VS.Vector Double -> Int getSampleSize = VS.length . VS.uniq . S.gsort -getNewMassDiagonalWithRescue :: Int -> Double -> Double -> Double-getNewMassDiagonalWithRescue sampleSize massOld massEstimate+rescueAfter :: Double -> Double+rescueAfter y+  | massMin > abs y = signum y * massMin+  | abs y > massMax = signum y * massMax+  | otherwise = y++rescueBefore :: (Double -> Double -> Double) -> Double -> Double -> Double+rescueBefore f old new+  -- Be permissive with NaN and infinite values.+  | isNaN new = old+  | isInfinite new = old+  | otherwise = f old new++rescue :: (Double -> Double -> Double) -> Double -> Double -> Double+rescue f old new = rescueAfter $ rescueBefore f old new++-- -- I do not know where I got this function from, but it works pretty well!+-- interpolate :: Double -> Double -> Double+-- interpolate old new = interSqrt ** 2+--   where+--     interSqrt = recip 3 * (sqrt old + 2 * sqrt new)++-- | This parameter plays the same role as @m@ in the dual averaging algorithm.+-- In the beginning, when @m@ is zero, the newly estimated masses will take full+-- precedence over the current masses. After some auto tuning steps, when @m@ is+-- larger, the newly estimated masses influence the current masses only+-- slightly.+newtype SmoothingParameter = SmoothingParameter Natural++-- Another interpolation function I came up with. It is pretty cool, because+-- (similar to above) it interpolates the square roots, which is what we want.+interpolate' :: SmoothingParameter -> Double -> Double -> Double+interpolate' (SmoothingParameter m) oldSquared newSquared = finSign * (fin ** 2)+  where+    oldSign = signum oldSquared+    newSign = signum newSquared+    sqrt' = sqrt . abs+    old = oldSign * sqrt' oldSquared+    new = newSign * sqrt' newSquared+    -- The new mass will be the second last boundary. That is, the larger the+    -- number of bins is, the more informative the new mass is compared to the+    -- old mass.+    nbins = max (100 - fromIntegral m) 2+    fin = recip nbins * (old + (nbins - 1) * new)+    finSign = signum fin++getNewMassDiagonalWithSampleSize :: SmoothingParameter -> Int -> Double -> Double -> Double+getNewMassDiagonalWithSampleSize m sampleSize massOld massEstimate   | sampleSize < samplesDiagonalMin = massOld-  -- Be permissive with NaN and negative diagonal masses. Diagonal masses are-  -- variances which are strictly positive.-  | isNaN massEstimate = massOld+  -- Diagonal masses are variances which are strictly positive.   | massEstimate <= 0 = massOld-  | massMin > massNew = massMin-  | massNew > massMax = massMax-  | otherwise = massNew-  where-    massNewSqrt = recip 3 * (sqrt massOld + 2 * sqrt massEstimate)-    massNew = massNewSqrt ** 2+  | otherwise = rescue (interpolate' m) massOld massEstimate +getNewMassDiagonal :: SmoothingParameter -> Double -> Double -> Double+getNewMassDiagonal m massOld massEstimate+  -- Diagonal masses are variances which are strictly positive.+  | massEstimate <= 0 = massOld+  | otherwise = rescue (interpolate' m) massOld massEstimate++getNewMassOffDiagonal :: SmoothingParameter -> Double -> Double -> Double+getNewMassOffDiagonal m = rescue (interpolate' m)+ -- The Cholesky decomposition, which is performed when sampling new momenta with -- 'generateMomenta', requires a positive definite covariance matrix. The -- Graphical Lasso algorithm finds positive definite covariance matrices, but@@ -214,7 +264,7 @@     m = L.cols a     b = L.unSym $ L.sym a     (_, s, v) = L.svd b-    h = (L.tr v) L.<> (L.diag s L.<> v)+    h = L.tr v L.<> (L.diag s L.<> v)     a2 = L.scale 0.5 (b + h)     a3 = L.unSym $ L.sym a2     isPositiveDefinite = isJust . L.mbChol . L.trustSym@@ -231,6 +281,7 @@            in go x' (k + 1)  tuneDiagonalMassesOnly ::+  SmoothingParameter ->   -- Conversion from value to vector.   (a -> Positions) ->   -- Value vector.@@ -242,7 +293,7 @@ -- NOTE: Here, we lose time because we convert the states to vectors again, -- something that has already been done. But then, auto tuning is not a runtime -- determining factor.-tuneDiagonalMassesOnly toVec xs (ms, msI)+tuneDiagonalMassesOnly m toVec xs (ms, msI)   -- If not enough data is available, do not tune.   | VB.length xs < samplesDiagonalMin = (ms, msI)   | dimState /= dimMs = error "tuneDiagonalMassesOnly: Dimension mismatch."@@ -271,16 +322,40 @@     msDiagonalEstimate = VS.fromList $ map (recip . S.variance) $ L.toColumns xs''     msDiagonalNew =       VS.zipWith3-        getNewMassDiagonalWithRescue+        (getNewMassDiagonalWithSampleSize m)         sampleSizes         msDiagonalOld         msDiagonalEstimate  -- This value was carefully tuned using the example "hamiltonian". defaultGraphicalLassoPenalty :: Double-defaultGraphicalLassoPenalty = 0.4+defaultGraphicalLassoPenalty = 0.3 +interpolateElementWise :: SmoothingParameter -> L.Matrix Double -> L.Matrix Double -> L.Matrix Double+interpolateElementWise m old new+  | mO /= mN = err "different number of rows"+  | nO /= nN = err "different number of columns"+  | mO /= nO = err "not square"+  | mO < 1 = err "empty matrix"+  | otherwise = L.build (mO, nO) f+  where+    mO = L.rows old+    nO = L.cols old+    mN = L.rows new+    nN = L.cols new+    err msg = error $ "interpolateElementWise: " <> msg+    f iD jD =+      let -- This sucks a bit, because we need a function (e -> e -> e), and+          -- since the return type is Double, the indices are also Doubble.+          i = round iD+          j = round jD+          g = if i == j then getNewMassDiagonal m else getNewMassOffDiagonal m+          eO = old `L.atIndex` (i, j)+          eN = new `L.atIndex` (i, j)+       in g eO eN+ tuneAllMasses ::+  SmoothingParameter ->   -- Conversion from value to vector.   (a -> Positions) ->   -- Value vector.@@ -292,16 +367,16 @@ -- NOTE: Here, we lose time because we convert the states to vectors again, -- something that has already been done. But then, auto tuning is not a runtime -- determining factor.-tuneAllMasses toVec xs (ms, msI)+tuneAllMasses m toVec xs (ms, msI)   -- If not enough data is available, do not tune.   | VB.length xs < samplesDiagonalMin = (ms, msI)   -- If not enough data is available, only the diagonal masses are tuned.   | VB.length xs < samplesAllMinWith dimMs = fallbackDiagonal   | L.rank xs'' /= dimState = fallbackDiagonal   | dimState /= dimMs = error "tuneAllMasses: Dimension mismatch."-  | otherwise = (msPD', msPDI')+  | otherwise = (msNew, msINew)   where-    fallbackDiagonal = tuneDiagonalMassesOnly toVec xs (ms, msI)+    fallbackDiagonal = tuneDiagonalMassesOnly m toVec xs (ms, msI)     -- xs: Each element contains all parameters of one iteration.     -- xs': Each element is a vector containing all parameters changed by the     -- proposal of one iteration.@@ -314,16 +389,18 @@     dimState = VS.length $ VB.head xs'     dimMs = L.rows $ L.unSym ms     (_, ss, xsNormalized) = S.scale xs''-    sigmaNormalized =-      L.unSym $-        either error fst $-          S.graphicalLasso defaultGraphicalLassoPenalty xsNormalized-    -- Sigma is the inverted mass matrix.-    sigma = S.rescaleSWith ss sigmaNormalized-    msI' = toGMatrix sigma-    -- The masses should be positive definite, but sometimes they happen to be-    -- not because of numerical errors.-    ms' = L.sym $ cleanMatrix $ L.inv sigma+    -- The first value is the covariance matrix sigma, which the inverted mass+    -- matrix (precision matrix). However, we interpolate the new mass matrix+    -- using the old one and the new estimate, so we have to recalculate the+    -- covariance matrix anyways.+    (_, precNormalized) =+      either error id $+        S.graphicalLasso defaultGraphicalLassoPenalty xsNormalized+    ms' = S.rescalePWith ss (L.unSym precNormalized)+    -- Clean NaNs, infinities; ensure positive definiteness. The masses should+    -- be positive definite, but sometimes they happen to be not because of+    -- numerical errors.+    msNewDirty = interpolateElementWise m (L.unSym ms) ms'     -- Positive definite matrices are symmetric.-    msPD' = L.trustSym $ findClosestPositiveDefiniteMatrix $ L.unSym ms'-    msPDI' = if L.unSym ms' == L.unSym msPD' then msI' else getMassesI msPD'+    msNew = L.trustSym $ findClosestPositiveDefiniteMatrix $ cleanMatrix msNewDirty+    msINew = getMassesI msNew

@@ -79,7 +79,7 @@ -- Second function in Algorithm 3 and Algorithm 6, respectively in [4]. buildTreeWith ::   -- The exponent of the total energy of the starting state is used to-  -- calcaulate the expected acceptance rate 'Alpha'.+  -- calculate the expected acceptance rate 'Alpha'.   Log Double ->   MassesI ->   Target ->@@ -108,8 +108,8 @@                 n' = if u <= expEPot' * expEKin' then 1 else 0                 errorIsSmall = u < deltaMax * expETot'                 alpha' = expETot' / expETot0-                alpha = min 1.0 alpha'-+                -- Limit expected acceptance rate between 0 and 1.+                alpha = max 0 $ min 1 alpha'   -- Recursive case. This is complicated because the algorithm is written for an   -- imperative language, and because we have two stacked monads.   | otherwise = do@@ -138,10 +138,10 @@             Nothing -> pure Nothing             Just (qm'', pm'', qp'', pp'', q''', n''', a''', na''') -> do               b <- uniformRM (0, 1) g :: IO Double-              let q'''' = if b < fromIntegral n''' / (fromIntegral $ n' + n''') then q''' else q'+              let n'''' = n' + n'''+                  q'''' = if b < fromIntegral n''' / fromIntegral n'''' then q''' else q'                   a'''' = a' + a'''                   na'''' = na' + na'''-                  n'''' = n' + n'''                   -- Important: Check for U-turn. This formula differs from the                   -- formula using indicator functions in Algorithm 3. However,                   -- check Equation (4).@@ -152,7 +152,7 @@   where     buildTree = buildTreeWith expETot0 msI tfun g --- | Paramters of the NUTS proposal.+-- | Parameters of the NUTS proposal. -- -- Includes tuning parameters and tuning configuration. data NParams = NParams@@ -184,8 +184,8 @@  data IsNew   = Old-  | OldWith {_acceptanceCountsOld :: AcceptanceCounts}-  | NewWith {_acceptanceCountsNew :: AcceptanceCounts}+  | OldWith {_acceptanceRatesOld :: AcceptanceRates}+  | NewWith {_acceptanceRatesNew :: AcceptanceRates}  -- First function in Algorithm 3. nutsPFunction ::@@ -230,28 +230,28 @@           Nothing -> pure (y, isNew)           Just (qm'', pm'', qp'', pp'', y'', n'', a, na) -> do             let r = fromIntegral n'' / fromIntegral n :: Double-                ar = (exp $ ln a) / fromIntegral na :: Double-                getCounts s = max 0 $ min 100 $ round $ s * 100-                ac =-                  if ar >= 0-                    then let cs = getCounts ar in AcceptanceCounts cs (100 - cs)-                    else error $ "nutsPFunction: Acceptance rate negative."+                -- Individual expected acceptance rates are limited in+                -- 'buildTreeWith'.+                ar = max 0 $ exp $ ln a+                ars = AcceptanceRates ar na             isAccept <-               if r > 1.0                 then pure True                 else do                   b <- uniformRM (0, 1) g                   pure $ b < r-            let (y''', isNew') = if isAccept then (y'', NewWith ac) else (y, OldWith ac)+            let (y''', isNew') = if isAccept then (y'', NewWith ars) else (y, OldWith ars)                 isUTurn = let dq = (qp'' - qm'') in (dq * pm'' < 0) || (dq * pp'' < 0)             if isUTurn               then pure (y''', isNew')               else go qm'' pm'' qp'' pp'' (j + 1) y''' (n + n'') isNew'   (x', isNew) <- go q p q p 0 q 1 Old   pure $ case isNew of-    Old -> (ForceReject, Just $ AcceptanceCounts 0 100)-    OldWith ac -> (ForceReject, Just $ ac)-    NewWith ac -> (ForceAccept $ fromVec x', Just $ ac)+    -- NOTE: I am not sure if it is correct to set the expected acceptance rate+    -- to 0 when the no u-turn leapfrog integrator fails.+    Old -> (ForceReject, Just $ AcceptanceRates 0 1)+    OldWith ac -> (ForceReject, Just ac)+    NewWith ac -> (ForceAccept $ fromVec x', Just ac)   where     (HParamsI e _ ms _ _ msI mus) = hparamsi     (HStructure _ toVec fromVecWith) = hstruct@@ -299,7 +299,7 @@       tuner = do         tfun <- hTuningFunctionWith dim toVec htconf         let pfun = nutsPFunctionWithTuningParameters dim hstruct targetWith-        pure $ Tuner 1.0 ts True tfun pfun+        pure $ Tuner 1.0 ts True True tfun pfun    in case checkHStructureWith (hpsMasses hParamsI) hstruct of         Just err -> error err         Nothing -> nutsWith tuner