mcmc 0.7.0.1 → 0.8.0.0

raw patch · 25 files changed

+255/−276 lines, 25 filesdep −data-defaultdep −double-conversiondep −pretty-showPVP ok

version bump matches the API change (PVP)

Dependencies removed: data-default, double-conversion, pretty-show

API changes (from Hackage documentation)

Files

• ChangeLog.md +15/−0
• bench/Poisson.hs +1/−1
• mcmc.cabal +3/−6
• src/Mcmc/Acceptance.hs +12/−6
• src/Mcmc/Algorithm.hs +2/−2
• src/Mcmc/Algorithm/MC3.hs +21/−8
• src/Mcmc/Algorithm/MHG.hs +17/−8
• src/Mcmc/Chain/Chain.hs +2/−3
• src/Mcmc/Chain/Save.hs +13/−14
• src/Mcmc/Chain/Trace.hs +4/−17
• src/Mcmc/Cycle.hs +9/−11
• src/Mcmc/Internal/Shuffle.hs +3/−9
• src/Mcmc/Internal/SpecFunctions.hs +16/−13
• src/Mcmc/Jacobian.hs +2/−6
• src/Mcmc/Likelihood.hs +2/−6
• src/Mcmc/MarginalLikelihood.hs +1/−2
• src/Mcmc/Mcmc.hs +18/−24
• src/Mcmc/Monitor.hs +18/−26
• src/Mcmc/Monitor/Log.hs +2/−2
• src/Mcmc/Monitor/Parameter.hs +7/−9
• src/Mcmc/Monitor/ParameterBatch.hs +9/−12
• src/Mcmc/Posterior.hs +3/−7
• src/Mcmc/Prior.hs +3/−16
• src/Mcmc/Proposal.hs +70/−67
• src/Mcmc/Proposal/Hamiltonian/Internal.hs +2/−1

@@ -5,6 +5,21 @@ ## Unreleased changes  +## 0.8.0.0++-   Use the new bytestring realfloat builder.+-   Fix a small memory leak during burn in with auto tuning.+-   Save chain after burn in.+-   New example (Poisson; see the respective [blog post](https://dschrempf.github.io/coding/2022-06-28-sample-from-a-posterior-using-markov-chain-monte-carlo-algorithms-and-haskell/)).+-   Simplify prior, likelihood, posterior, and Jacobian types.+-   Do not export `grabble`.+-   Remove chain index.+-   Remove `lengthT`.+-   Remove `ProposalOrder` default instance.+-   Use GHC 9.2.4 by default.+-   Reduce package dependencies as much as possible.++ ## 0.7.0.1  -   Use random-1.2. This means, that the random number generator type has changed

@@ -49,7 +49,7 @@     ys = [1976.0 .. 1985.0]     m = sum ys / fromIntegral (length ys) -f :: (RealFloat a, Typeable a) => Int -> Double -> IG a -> LikelihoodG a+f :: (RealFloat a, Typeable a) => Int -> Double -> LikelihoodFunctionG (IG a) a f ft yr xs = poisson l ft   where     a = xs VB.! 0

@@ -1,6 +1,6 @@ cabal-version:      3.0 name:               mcmc-version:            0.7.0.1+version:            0.8.0.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>@@ -92,22 +92,19 @@     , ad     , aeson     , async-    , base               >=4.7 && <5+    , base            >=4.7 && <5     , bytestring     , circular     , containers-    , covariance         >=0.2-    , data-default+    , covariance      >=0.2     , directory     , dirichlet-    , double-conversion     , hmatrix     , log-domain     , math-functions     , microlens     , mwc-random     , parallel-    , pretty-show     , primitive     , random     , splitmix

@@ -83,25 +83,31 @@ pushAcceptanceCounts :: Ord k => k -> AcceptanceCounts -> Acceptance k -> Acceptance k pushAcceptanceCounts k c = Acceptance . M.adjust (addAcceptanceCounts c) k . fromAcceptance --- | Reset acceptance storage.+-- | Reset acceptance counts. resetA :: Ord k => Acceptance k -> Acceptance k resetA = emptyA . M.keys . fromAcceptance -transformKeys :: (Ord k1, Ord k2) => [k1] -> [k2] -> M.Map k1 v -> M.Map k2 v-transformKeys ks1 ks2 m = foldl' insrt M.empty $ zip ks1 ks2+transformKeys :: (Ord k1, Ord k2) => [(k1, k2)] -> M.Map k1 v -> M.Map k2 v+transformKeys ks m = foldl' insrt M.empty ks   where     insrt m' (k1, k2) = M.insert k2 (m M.! k1) m'  -- | 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 ks1 ks2 = Acceptance . transformKeys ks1 ks2 . fromAcceptance+transformKeysA :: (Ord k1, Ord k2) => [(k1, k2)] -> Acceptance k1 -> Acceptance k2+transformKeysA ks = Acceptance . transformKeys ks . fromAcceptance  -- | Acceptance counts and rate for a specific proposal. --+-- Return @Just (accepts, rejects, acceptance rate)@.+-- -- Return 'Nothing' if no proposals have been accepted or rejected (division by -- zero).-acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, AcceptanceRate)+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))

@@ -28,8 +28,8 @@   -- | Current iteration.   aIteration :: a -> Int -  -- | Check if the current state is invalid. At the moment this should just-  -- check if the posterior probability is zero or NaN.+  -- | Check if the current state is invalid. A state is invalid if the+  -- posterior probability is positive or negative infinite or NaN.   aIsInvalidState :: a -> Bool    -- | Sample the next state.

@@ -55,7 +55,6 @@ import Data.Aeson.TH import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL-import qualified Data.Double.Conversion.ByteString as BC import Data.List import qualified Data.Map.Strict as M import Data.Time@@ -80,6 +79,7 @@ import Mcmc.Settings import Numeric.Log hiding (sum) import System.Random.Stateful+import Text.Printf  -- | Total number of parallel chains. --@@ -251,14 +251,14 @@   IO (MHG a) initMHG prf lhf i beta a   | i < 0 = error "initMHG: Chain index negative."-  -- Only set the id for the cold chain.-  | i == 0 = return $ MHG $ c {chainId = Just 0}+  -- Do nothing for the cold chain.+  | 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       -- pointing to the link of the cold chain, but this has no effect.       t' <- pushT l t-      return $ MHG $ c {chainId = Just i, trace = t'}+      return $ MHG $ c {trace = t'}   where     a' = setReciprocalTemperature prf lhf beta a     c = fromMHG a'@@ -321,6 +321,8 @@  -- | Load an MC3 algorithm. --+-- Also create a backup of the save.+-- -- See 'Mcmc.Mcmc.mcmcContinue'. mc3Load ::   FromJSON a =>@@ -331,8 +333,11 @@   AnalysisName ->   IO (MC3 a) mc3Load pr lh cc mn nm = do-  savedMC3 <- eitherDecode . decompress <$> BL.readFile (mc3Fn nm)+  savedMC3 <- eitherDecode . decompress <$> BL.readFile fn   either error (fromSavedMC3 pr lh cc mn) savedMC3+  where+    -- fnBak = mc3Fn $ AnalysisName $ (fromAnalysisName nm ++ ".bak")+    fn = mc3Fn nm  -- I call the chains left and right, because it is easy to think about them as -- being left and right. Of course, the left chain may also have a larger index@@ -541,7 +546,7 @@         Nothing -> []         Just ar ->           [ "MC3: Average acceptance rate across all chains: "-              <> BL.fromStrict (BC.toFixed 2 ar)+              <> BB.toLazyByteString (BB.formatDouble (BB.standard 2) ar)               <> "."           ]       ++ [ "MC3: Reciprocal temperatures of the chains: " <> BL.intercalate ", " bsB <> ".",@@ -572,17 +577,25 @@     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 (BL.fromStrict . BC.toFixed 2) $ U.toList bs+    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     n = fromNChains $ mc3NChains $ mc3Settings a+    proposalHLine = BL.replicate (BL.length proposalHeader) '-'  -- No extra monitors are opened. mc3OpenMonitors :: ToJSON a => AnalysisName -> ExecutionMode -> MC3 a -> IO (MC3 a) mc3OpenMonitors nm em a = do-  mhgs' <- V.mapM (aOpenMonitors nm em) (mc3MHGChains a)+  mhgs' <- V.imapM mhgOpenMonitors (mc3MHGChains a)   return $ a {mc3MHGChains = mhgs'}+  where+    mhgOpenMonitors i (MHG c) = do+      m' <- mOpen pre suf em $ monitor c+      pure $ MHG c {monitor = m'}+      where+        pre = fromAnalysisName nm+        suf = printf "%02d" i  mc3ExecuteMonitors ::   ToJSON a =>

@@ -52,7 +52,6 @@ import Mcmc.Settings import Numeric.Log import System.Random.Stateful-import Text.Printf import Prelude hiding (cycle)  -- | The MHG algorithm.@@ -111,7 +110,7 @@   -- a monad.   tr <- replicateT tl l0   gm <- newIOGenM g-  return $ MHG $ Chain Nothing l0 0 tr ac gm 0 pr lh cc mn+  return $ MHG $ Chain l0 0 tr ac gm 0 pr lh cc mn   where     l0 = Link i0 (pr i0) (lh i0)     ac = emptyA $ ccProposals cc@@ -132,6 +131,8 @@  -- | Load an MHG algorithm. --+-- Also create a backup of the save.+-- -- See 'Mcmc.Mcmc.mcmcContinue'. mhgLoad ::   FromJSON a =>@@ -143,8 +144,10 @@   IO (MHG a) mhgLoad = mhgLoadWith fromSavedChain --- | See 'mhgLoad' but do not perform sanity checks.+-- | Like 'mhgLoad' but do not perform sanity checks. --+-- Also create a backup of the save.+-- -- Useful when restarting a run with changed prior function, likelihood function -- or proposals. Use with care! mhgLoadUnsafe ::@@ -168,9 +171,12 @@   AnalysisName ->   IO (MHG a) mhgLoadWith f pr lh cc mn nm = do-  savedChain <- eitherDecode . decompress <$> BL.readFile (mhgFn nm)+  savedChain <- eitherDecode . decompress <$> BL.readFile fn   chain <- either error (f pr lh cc mn) savedChain   return $ MHG chain+  where+    -- fnBak = mhgFn $ AnalysisName $ (fromAnalysisName nm ++ ".bak")+    fn = mhgFn nm  -- | MHG ratios are stored in log domain. type MHGRatio = Log Double@@ -342,14 +348,17 @@     cc = cycle c     ac = acceptance c -mhgOpenMonitors :: AnalysisName -> ExecutionMode -> MHG a -> IO (MHG a)+mhgOpenMonitors ::+  AnalysisName ->+  ExecutionMode ->+  MHG a ->+  IO (MHG a) mhgOpenMonitors nm em (MHG c) = do-  m' <- mOpen pre suf em m-  return $ MHG c {monitor = m'}+  m' <- mOpen pre "" em m+  pure $ MHG c {monitor = m'}   where     m = monitor c     pre = fromAnalysisName nm-    suf = maybe "" (printf "%02d") $ chainId c  mhgExecuteMonitors ::   Verbosity ->

@@ -82,13 +82,12 @@ data Chain a = Chain   { -- Variables; saved. -    -- | Chain index; useful if more chains are run.-    chainId :: Maybe Int,     -- | The current 'Link' of the chain combines the current state and the     -- current likelihood. The link is updated after a proposal has been     -- executed.     link :: Link a,-    -- | The current iteration or completed number of cycles.+    -- | The current iteration. In contrast to the link, the iteration is+    -- updated only after all proposals in the cycle have been executed.     iteration :: Int,     -- | 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.

@@ -1,6 +1,12 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE TemplateHaskell #-} +-- |+-- Module      :  Mcmc.Chain.Save+-- Description :  Save and load a Markov chain+-- Copyright   :  2021 Dominik Schrempf+-- License     :  GPL-3.0-or-later+-- -- Maintainer  :  dominik.schrempf@gmail.com -- Stability   :  unstable -- Portability :  portable@@ -10,12 +16,6 @@ -- Save and load chains. It is easy to save and restore the current state and -- likelihood (or the trace), but it is not feasible to store all the proposals -- and so on, so they have to be provided again when continuing a run.---- |--- Module      :  Mcmc.Chain.Save--- Description :  Save and load a Markov chain--- Copyright   :  2021 Dominik Schrempf--- License     :  GPL-3.0-or-later module Mcmc.Chain.Save   ( SavedChain (..),     toSavedChain,@@ -48,8 +48,7 @@ -- -- See 'toSavedChain'. data SavedChain a = SavedChain-  { savedId :: Maybe Int,-    savedLink :: Link a,+  { savedLink :: Link a,     savedIteration :: Int,     savedTrace :: C.Stack VB.Vector (Link a),     savedAcceptance :: Acceptance Int,@@ -64,13 +63,13 @@ toSavedChain ::   Chain a ->   IO (SavedChain a)-toSavedChain (Chain ci it i tr ac g _ _ _ cc _) = do+toSavedChain (Chain it i tr ac g _ _ _ cc _) = do   g' <- saveGen g   tr' <- freezeT tr-  return $ SavedChain ci it i tr' ac' g' ts+  return $ SavedChain it i tr' ac' g' ts   where     ps = ccProposals cc-    ac' = transformKeysA ps [0 ..] ac+    ac' = transformKeysA (zip ps [0 ..]) ac     ts =       [ (\t -> (tTuningParameter t, tAuxiliaryTuningParameters t)) <$> mt         | mt <- map prTuner ps@@ -133,12 +132,12 @@   Monitor a ->   SavedChain a ->   IO (Chain a)-fromSavedChainUnsafe pr lh cc mn (SavedChain ci it i tr ac' g' ts) = do+fromSavedChainUnsafe pr lh cc mn (SavedChain it i tr ac' g' ts) = do   g <- loadGen g'   tr' <- thawT tr-  return $ Chain ci it i tr' ac g i pr lh cc' mn+  return $ Chain it i tr' ac g i pr lh cc' mn   where-    ac = transformKeysA [0 ..] (ccProposals cc) ac'+    ac = transformKeysA (zip [0 ..] (ccProposals cc)) ac'     tunePs mt p = case mt of       Nothing -> p       Just (x, xs) -> either (error . (<> err)) id $ tuneWithTuningParameters x xs p

@@ -13,7 +13,6 @@   ( Trace,     replicateT,     fromVectorT,-    lengthT,     pushT,     headT,     takeT,@@ -48,10 +47,6 @@ fromVectorT :: VB.Vector (Link a) -> IO (Trace a) fromVectorT xs = Trace <$> C.fromVector xs --- | Get the length of the trace.-lengthT :: Trace a -> Int-lengthT = C.size . fromTrace- -- | Push a 'Link' on the 'Trace'. pushT :: Link a -> Trace a -> IO (Trace a) pushT x t = do@@ -59,27 +54,19 @@   return $ Trace s' {-# INLINEABLE pushT #-} --- | Get the most recent link of the trace.------ See 'C.get'.+-- | Get the most recent link of the trace (see 'C.get'). headT :: Trace a -> IO (Link a) headT = C.get . fromTrace {-# INLINEABLE headT #-} --- | Get the k most recent links of the trace.------ See 'C.take'.+-- | Get the k most recent links of the trace (see 'C.take'). takeT :: Int -> Trace a -> IO (VB.Vector (Link a)) takeT k = C.take k . fromTrace --- | Freeze the mutable trace for storage.------ See 'C.freeze'.+-- | Freeze the mutable trace for storage (see 'C.freeze'). freezeT :: Trace a -> IO (C.Stack VB.Vector (Link a)) freezeT = C.freeze . fromTrace --- | Thaw a circular stack.------ See 'See.thaw'.+-- | Thaw a circular stack (see 'C.thaw'). thawT :: C.Stack VB.Vector (Link a) -> IO (Trace a) thawT t = Trace <$> C.thaw t

@@ -23,14 +23,12 @@     autoTuneCycle,      -- * Output-    proposalHLine,     summarizeCycle,   ) where  import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL-import Data.Default import Data.List import qualified Data.Map.Strict as M import qualified Data.Vector as VB@@ -63,8 +61,6 @@     SequentialReversibleO   deriving (Eq, Show) -instance Default Order where def = RandomO- -- Describe the order. describeOrder :: Order -> BL.ByteString describeOrder RandomO = "The proposals are executed in random order."@@ -98,13 +94,13 @@     ccRequireTrace :: Bool   } --- | Create a 'Cycle' from a list of 'Proposal's.+-- | Create a 'Cycle' from a list of 'Proposal's; use 'RandomO', but see 'setOrder'. cycleFromList :: [Proposal a] -> Cycle a cycleFromList [] =   error "cycleFromList: Received an empty list but cannot create an empty Cycle." cycleFromList xs =   if length uniqueXs == length xs-    then Cycle xs def (any needsTrace xs)+    then Cycle xs RandomO (any needsTrace xs)     else error $ "\n" ++ msg ++ "cycleFromList: Proposals are not unique."   where     uniqueXs = nub xs@@ -162,7 +158,12 @@     once = sum $ map (fromPWeight . prWeight) xs'  -- See 'tuneWithTuningParameters' and 'Tuner'.-tuneWithChainParameters :: TuningType -> AcceptanceRate -> Maybe (VB.Vector a) -> Proposal a -> Either String (Proposal a)+tuneWithChainParameters ::+  TuningType ->+  AcceptanceRate ->+  Maybe (VB.Vector a) ->+  Proposal a ->+  Either String (Proposal a) tuneWithChainParameters b ar mxs p = case prTuner p of   Nothing -> Right p   Just (Tuner t ts rt fT _) -> case (rt, mxs) of@@ -191,10 +192,6 @@         Just (Just x) -> either error id $ tuneWithChainParameters b x mxs p         _ -> p --- | Horizontal line of proposal summaries.-proposalHLine :: BL.ByteString-proposalHLine = BL.replicate (BL.length proposalHeader) '-'- -- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance rates. summarizeCycle :: IterationMode -> Acceptance (Proposal a) -> Cycle a -> BL.ByteString summarizeCycle m a c =@@ -223,3 +220,4 @@       1 -> nProposalsStr <> " proposal is performed per iteration."       _ -> nProposalsStr <> " proposals are performed per iterations."     ar pr = acceptanceRate pr a+    proposalHLine = BL.replicate (BL.length proposalHeader) '-'

@@ -13,7 +13,6 @@ -- From https://wiki.haskell.org/Random_shuffle. module Mcmc.Internal.Shuffle   ( shuffle,-    grabble,   ) where @@ -23,19 +22,14 @@ import qualified Data.Vector.Mutable as M import System.Random.Stateful --- Fisher-Yates shuffle. See also--- 'System.Random.MWC.Distributions.uniformPermutation' which is a little--- cleaner, in my opinion. However, I would like to move away from MWC so I--- leave the custom implementation for now.------ NOTE: I do not need to move away from MWC anymore true because I can use--- SplitMix with it.+-- TODO: Remove shuffle, use System.Random.MWC.Distributions.uniformShuffle and+-- vectors.  -- | Shuffle a vector. shuffle :: StatefulGen g m => [a] -> g -> m [a] shuffle xs = grabble xs (length xs) --- | @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@+-- @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@ -- elements from @xs@, without replacement, and that @m@ times. grabble :: StatefulGen g m => [a] -> Int -> g -> m [a] grabble xs m gen = do

@@ -39,7 +39,7 @@   | otherwise = logGammaNonDouble z {-# SPECIALIZE logGammaG :: Double -> Double #-} --- | See 'Numeric.SpecFunctions.logGamma'.+-- See 'Numeric.SpecFunctions.logGamma'. logGammaNonDouble :: RealFloat a => a -> a logGammaNonDouble z   | z <= 0 = 1 / 0@@ -53,10 +53,11 @@  lgamma1_15G :: RealFloat a => a -> a -> a lgamma1_15G zm1 zm2 =-  r * y + r-    * ( evaluatePolynomial zm1 tableLogGamma_1_15PG-          / evaluatePolynomial zm1 tableLogGamma_1_15QG-      )+  r * y+    + r+      * ( evaluatePolynomial zm1 tableLogGamma_1_15PG+            / evaluatePolynomial zm1 tableLogGamma_1_15QG+        )   where     r = zm1 * zm2     y = 0.52815341949462890625@@ -89,10 +90,11 @@  lgamma15_2G :: RealFloat a => a -> a -> a lgamma15_2G zm1 zm2 =-  r * y + r-    * ( evaluatePolynomial (-zm2) tableLogGamma_15_2PG-          / evaluatePolynomial (-zm2) tableLogGamma_15_2QG-      )+  r * y+    + r+      * ( evaluatePolynomial (-zm2) tableLogGamma_15_2PG+            / evaluatePolynomial (-zm2) tableLogGamma_15_2QG+        )   where     r = zm1 * zm2     y = 0.452017307281494140625@@ -133,10 +135,11 @@  lgamma2_3G :: RealFloat a => a -> a lgamma2_3G z =-  r * y + r-    * ( evaluatePolynomial zm2 tableLogGamma_2_3PG-          / evaluatePolynomial zm2 tableLogGamma_2_3QG-      )+  r * y+    + r+      * ( evaluatePolynomial zm2 tableLogGamma_2_3PG+            / evaluatePolynomial zm2 tableLogGamma_2_3QG+        )   where     r = zm2 * (z + 1)     zm2 = z - 2

@@ -11,7 +11,6 @@ -- Creation date: Tue May 31 09:37:54 2022. module Mcmc.Jacobian   ( Jacobian,-    JacobianG,     JacobianFunction,     JacobianFunctionG,   )@@ -22,11 +21,8 @@ -- | Absolute value of the determinant of the Jacobian matrix. type Jacobian = Log Double --- | Generalized Jacobian.-type JacobianG a = Log a- -- | Function calculating the 'Jacobian'.-type JacobianFunction a = JacobianFunctionG a Double+type JacobianFunction a = a -> Log Double  -- | Function calculating the 'Jacobian'.-type JacobianFunctionG a b = a -> JacobianG b+type JacobianFunctionG a b = a -> Log b

@@ -11,7 +11,6 @@ -- Creation date: Wed Mar  3 11:39:04 2021. module Mcmc.Likelihood   ( Likelihood,-    LikelihoodG,     LikelihoodFunction,     LikelihoodFunctionG,     noLikelihood,@@ -23,14 +22,11 @@ -- | Likelihood values are stored in log domain. type Likelihood = Log Double --- | Generalized likelihood.-type LikelihoodG a = Log a- -- | Likelihood function.-type LikelihoodFunction a = LikelihoodFunctionG a Double+type LikelihoodFunction a = a -> Log Double  -- | Generalized likelihood function.-type LikelihoodFunctionG a b = a -> LikelihoodG b+type LikelihoodFunctionG a b = a -> Log b  -- | Flat likelihood function. Useful for testing and debugging. noLikelihood :: RealFloat b => LikelihoodFunctionG a b

@@ -48,7 +48,6 @@ import System.Directory import System.Random.Stateful import Text.Printf-import Text.Show.Pretty import Prelude hiding (cycle)  -- | Marginal likelihood values are stored in log domain.@@ -393,7 +392,7 @@         logInfoStartingTime         logInfoB "Estimate marginal likelihood."         logDebugB "marginalLikelihood: The marginal likelihood settings are:"-        logDebugS $ ppShow s+        logDebugS $ show s         val <- case mlAlgorithm s of           ThermodynamicIntegration -> tiWrapper s prf lhf cc mn i0 g           SteppingStoneSampling -> sssWrapper s prf lhf cc mn i0 g

@@ -63,21 +63,11 @@   putStrLn "Closing output files."   _ <- aCloseMonitors a   closeEnvironment e-  putStrLn "Saving settings."-  let s = settings e-  settingsSave s-  putStrLn "Saving compressed MCMC analysis."-  putStrLn "For long traces, or complex objects, this may take a while."-  let nm = sAnalysisName s-  aSave nm a-  putStrLn "Markov chain saved. Analysis can be continued."+  runReaderT (mcmcSave a) e   putStrLn "Graceful termination successful."   putStrLn "Rethrowing error."   throw err --- XXX: Exception handling. Is it enough to mask execution of monitors and catch--- UserInterrupt during iterations?- mcmcExecuteMonitors :: Algorithm a => a -> MCMC () mcmcExecuteMonitors a = do   e <- ask@@ -96,25 +86,29 @@   | otherwise = do       e <- ask       p <- sParallelizationMode . settings <$> ask-      -- NOTE: User interrupt is handled during iterations.+      -- NOTE: User interrupt is only handled during iterations.       a' <- liftIO $ catch (aIterate m p a) (mcmcExceptionHandler e a)+      -- NOTE: We may want to mask execution of monitors and catch exceptions+      -- even then?       mcmcExecuteMonitors a'       mcmcIterate m (n - 1) a'  mcmcNewRun :: Algorithm a => a -> MCMC a mcmcNewRun a = do   s <- reader settings-  logInfoB "Start new MCMC sampler."+  logInfoB "Starting new MCMC sampler."   logInfoB "Initial state."   logInfoB $ aStdMonitorHeader a   mcmcExecuteMonitors a   when (aIsInvalidState a) (logWarnB "The initial state is invalid!")   a' <- mcmcBurnIn a-  logInfoS $ "Clean 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."+  mcmcSave a''   let i = fromIterations $ sIterations s-  logInfoS $ "Run chain for " ++ show i ++ " iterations."+  logInfoS $ "Running chain for " ++ show i ++ " iterations."   logInfoB $ aStdMonitorHeader a''   mcmcIterate AllProposals i a'' @@ -133,7 +127,7 @@   when (iCurrent < iBurnIn) $ error "mcmcContinueRun: Can not continue burn in."   let di = iTotal - iCurrent   logInfoB $ aSummarizeCycle AllProposals a-  logInfoS $ "Run chain for " ++ show di ++ " iterations."+  logInfoS $ "Running chain for " ++ show di ++ " iterations."   logInfoB $ aStdMonitorHeader a   mcmcIterate AllProposals di a @@ -147,7 +141,7 @@       return a     BurnInWithoutAutoTuning n -> do       logInfoB $ aSummarizeCycle AllProposals a-      logInfoS $ "Burn in for " <> show n <> " iterations."+      logInfoS $ "Burning in for " <> show n <> " iterations."       logInfoS "Auto tuning is disabled."       logInfoB $ aStdMonitorHeader a       a' <- mcmcIterate AllProposals n a@@ -157,7 +151,7 @@       return a''     BurnInWithAutoTuning n t -> do       logInfoB $ aSummarizeCycle AllProposals a-      logInfoS $ "Burn in for " ++ show n ++ " iterations."+      logInfoS $ "Burning in for " ++ show n ++ " iterations."       logInfoS $ "Auto tuning is enabled with a period of " ++ show t ++ "."       logInfoB $ aStdMonitorHeader a       let (m, r) = n `divMod` t@@ -168,7 +162,7 @@       logInfoB "Burn in finished."       return a'     BurnInWithCustomAutoTuning xs ys -> do-      logInfoS $ "Burn in for " ++ show (sum xs + sum ys) ++ " iterations."+      logInfoS $ "Burning in for " ++ show (sum xs + sum ys) ++ " iterations."       a' <-         if null xs           then do@@ -188,7 +182,7 @@ -- Auto tune the proposals. mcmcAutotune :: Algorithm a => TuningType -> Int -> a -> MCMC a mcmcAutotune NormalTuningStep n a = do-  logDebugB "Auto tune."+  logDebugB "Intermediate auto tune."   liftIO $ aAutoTune NormalTuningStep n a mcmcAutotune LastTuningStep n a = do   logDebugB "Last auto tune."@@ -226,15 +220,15 @@ mcmcSave a = do   s <- reader settings   case sSaveMode s of-    NoSave -> logInfoB "Do not save the MCMC analysis."+    NoSave -> logInfoB "NoSave set; Do not save the MCMC analysis."     Save -> do-      logInfoB "Save settings."+      logInfoB "Saving settings."       liftIO $ settingsSave s       let nm = sAnalysisName s-      logInfoB "Save compressed MCMC analysis."+      logInfoB "Saving compressed MCMC analysis."       logInfoB "For long traces, or complex objects, this may take a while."       liftIO $ aSave nm a-      logInfoB "Markov chain saved."+      logInfoB "Markov chain saved. Analysis can be continued."  -- Report and finish up. mcmcClose :: Algorithm a => a -> MCMC a

@@ -62,24 +62,23 @@     mStdOut :: MonitorStdOut a,     -- | Monitors writing to files.     mFiles :: [MonitorFile a],-    -- | Monitors calculating batch means and-    -- writing to files.+    -- | Monitors calculating summary statistics over the last batch of+    -- iterations and writing to files.     mBatches :: [MonitorBatch a]   } --- | Monitor period.+-- | Execute monitor every given number of iterations. type Period = Int --- | Do not monitor parameters.------ Monitor prior and likelihood with given period.+-- | Only monitor prior, likelihood and posterior functions with given period.+-- Do not monitor parameters. simpleMonitor :: Period -> Monitor a simpleMonitor p   | p < 1 = error "simpleMonitor: Monitor period must be 1 or larger."   | otherwise =       Monitor (MonitorStdOut [] p) [] [] --- | Monitor to standard output; constructed with 'monitorStdOut'.+-- | Monitor to standard output; construct with 'monitorStdOut'. data MonitorStdOut a = MonitorStdOut   { msParams :: [MonitorParameter a],     msPeriod :: Period@@ -107,7 +106,7 @@  -- | Header of monitor to standard output. msHeader :: MonitorStdOut a -> BL.ByteString-msHeader m = BL.intercalate "\n" [row, sep]+msHeader m = row <> "\n" <> sep   where     row =       msRenderRow $@@ -130,10 +129,6 @@   let dt = ct `diffUTCTime` st       -- NOTE: Don't evaluate this when i == ss.       timePerIter = dt / fromIntegral (i - ss)-      -- -- Always 0; doesn't make much sense.-      -- tpi = if (i - ss) < 10-      --       then ""-      --       else renderDurationS timePerIter       eta =         if (i - ss) < 10           then ""@@ -154,9 +149,6 @@   IO (Maybe BL.ByteString) msExec i it ss st j m   | i `mod` msPeriod m /= 0 = return Nothing-  -- -- | i `mod` (msPeriod m * 100) == 0 = do-  -- --   l <- msDataLine i it ss st j m-  -- --   return $ Just $ msHeader m <> "\n" <> l   | otherwise = Just <$> msDataLine i it ss st j m  -- | Monitor to a file; constructed with 'monitorFile'.@@ -216,18 +208,18 @@       Just h ->         BL.hPutStrLn h $           mfRenderRow $-            BL.pack (show i) :-            renderLog p :-            renderLog l :-            renderLog (p * l) :-              [BB.toLazyByteString $ mpFunc mp x | mp <- mfParams m]+            BL.pack (show i)+              : renderLog p+              : renderLog l+              : renderLog (p * l)+              : [BB.toLazyByteString $ mpFunc mp x | mp <- mfParams m]  mfClose :: MonitorFile a -> IO () mfClose m = case mfHandle m of   Just h -> hClose h   Nothing -> error $ "mfClose: File was not opened for monitor " <> mfName m <> "." --- | Batch size.+-- | Size of the batch used to calculate summary statistics. type BatchSize = Int  -- | Batch monitor to a file.@@ -303,11 +295,11 @@             mlos = mean los         BL.hPutStrLn h $           mfRenderRow $-            BL.pack (show i) :-            renderLog mlps :-            renderLog mlls :-            renderLog mlos :-              [BB.toLazyByteString $ mbpFunc mbp (VB.map state xs) | mbp <- mbParams m]+            BL.pack (show i)+              : renderLog mlps+              : renderLog mlls+              : renderLog mlos+              : [BB.toLazyByteString $ mbpFunc mbp (VB.map state xs) | mbp <- mbParams m]  mbClose :: MonitorBatch a -> IO () mbClose m = case mbHandle m of

@@ -14,11 +14,11 @@   ) where +import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL-import qualified Data.Double.Conversion.ByteString as BC import Numeric.Log  -- | Print a log value. renderLog :: Log Double -> BL.ByteString-renderLog = BL.fromStrict . BC.toFixed 8 . ln+renderLog = BB.toLazyByteString . BB.formatDouble (BB.standard 8) . ln {-# INLINEABLE renderLog #-}

@@ -18,12 +18,11 @@     monitorInt,     monitorDouble,     monitorDoubleF,-    monitorDoubleE,+    monitorDoubleS,   ) where  import qualified Data.ByteString.Builder as BB-import qualified Data.Double.Conversion.ByteString as BC import Data.Functor.Contravariant  -- XXX: 'MonitorParameter' has a drawback. Extracting and monitoring multiple@@ -60,19 +59,18 @@   -- | Name.   String ->   MonitorParameter Double-monitorDouble n = MonitorParameter n (BB.byteString . BC.toFixed 8)+monitorDouble n = MonitorParameter n (BB.formatDouble $ BB.standard 8) --- | Monitor 'Double' with full precision computing the shortest string of--- digits that correctly represent the number.+-- | Monitor 'Double' with full precision. monitorDoubleF ::   -- | Name.   String ->   MonitorParameter Double-monitorDoubleF n = MonitorParameter n (BB.byteString . BC.toShortest)+monitorDoubleF n = MonitorParameter n BB.doubleDec --- | Monitor 'Double' in exponential format and half precision.-monitorDoubleE ::+-- | Monitor 'Double' in scientific format.+monitorDoubleS ::   -- | Name.   String ->   MonitorParameter Double-monitorDoubleE n = MonitorParameter n (BB.byteString . BC.toExponential 8)+monitorDoubleS n = MonitorParameter n (BB.formatDouble BB.scientific)

@@ -22,12 +22,11 @@     (>$<),     monitorBatchMean,     monitorBatchMeanF,-    monitorBatchMeanE,+    monitorBatchMeanS,   ) where  import qualified Data.ByteString.Builder as BB-import qualified Data.Double.Conversion.ByteString as BC import Data.Functor.Contravariant import qualified Data.Vector as VB @@ -69,32 +68,30 @@   -- | Name.   String ->   MonitorParameterBatch a-monitorBatchMean n = MonitorParameterBatch n (BB.byteString . BC.toFixed 8 . mean)+monitorBatchMean n = MonitorParameterBatch n (BB.formatDouble (BB.standard 8) . mean) {-# SPECIALIZE monitorBatchMean :: String -> MonitorParameterBatch Int #-} {-# SPECIALIZE monitorBatchMean :: String -> MonitorParameterBatch Double #-}  -- | Batch mean monitor. ----- Print the mean with full precision computing the shortest string of digits--- that correctly represent the number.+-- Print the mean with full precision. monitorBatchMeanF ::   Real a =>   -- | Name.   String ->   MonitorParameterBatch a-monitorBatchMeanF n = MonitorParameterBatch n (BB.byteString . BC.toShortest . mean)+monitorBatchMeanF n = MonitorParameterBatch n (BB.doubleDec . mean) {-# SPECIALIZE monitorBatchMeanF :: String -> MonitorParameterBatch Int #-} {-# SPECIALIZE monitorBatchMeanF :: String -> MonitorParameterBatch Double #-}  -- | Batch mean monitor. ----- Print the real float parameters such as 'Double' with scientific notation and--- eight decimal places.-monitorBatchMeanE ::+-- Print the real float parameters such as 'Double' with scientific notation.+monitorBatchMeanS ::   Real a =>   -- | Name.   String ->   MonitorParameterBatch a-monitorBatchMeanE n = MonitorParameterBatch n (BB.byteString . BC.toExponential 8 . mean)-{-# SPECIALIZE monitorBatchMeanE :: String -> MonitorParameterBatch Int #-}-{-# SPECIALIZE monitorBatchMeanE :: String -> MonitorParameterBatch Double #-}+monitorBatchMeanS n = MonitorParameterBatch n (BB.formatDouble BB.scientific . mean)+{-# SPECIALIZE monitorBatchMeanS :: String -> MonitorParameterBatch Int #-}+{-# SPECIALIZE monitorBatchMeanS :: String -> MonitorParameterBatch Double #-}

@@ -11,7 +11,6 @@ -- Creation date: Fri May 28 12:26:35 2021. module Mcmc.Posterior   ( Posterior,-    PosteriorG,     PosteriorFunction,     PosteriorFunctionG,   )@@ -20,13 +19,10 @@ import Numeric.Log  -- | Posterior values are stored in log domain.-type Posterior = PosteriorG Double---- | Generalized posterior.-type PosteriorG a = Log a+type Posterior = Log Double  -- | Posterior function.-type PosteriorFunction a = PosteriorFunctionG a Double+type PosteriorFunction a = a -> Log Double  -- | Generalized posterior function.-type PosteriorFunctionG a b = a -> PosteriorG b+type PosteriorFunctionG a b = a -> Log b

@@ -13,7 +13,6 @@ -- Creation date: Thu Jul 23 13:26:14 2020. module Mcmc.Prior   ( Prior,-    PriorG,     PriorFunction,     PriorFunctionG, @@ -50,26 +49,14 @@ import Mcmc.Statistics.Types import Numeric.Log --- TODO (high): Think about using a "structure" variable.------ For example,------ type PriorFunctionG s a = s a -> PriorG a------ Many of the prior functions would need to use the Identity Functor. This may--- be slow.- -- | Prior values are stored in log domain.-type Prior = PriorG Double---- | Generalized prior.-type PriorG a = Log a+type Prior = Log Double  -- | Prior function.-type PriorFunction a = PriorFunctionG a Double+type PriorFunction a = a -> Log Double  -- | Generalized prior function.-type PriorFunctionG a b = a -> PriorG b+type PriorFunctionG a b = a -> Log b  -- | Flat prior function. Useful for testing and debugging. noPrior :: RealFloat b => PriorFunctionG a b

@@ -1,3 +1,4 @@+{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DerivingVia #-} {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE RankNTypes #-}@@ -26,9 +27,6 @@     PResult (..),     Jacobian,     JacobianFunction,-    (@~),-    liftProposal,-    liftProposalWith,     PFunction,     createProposal, @@ -47,6 +45,11 @@     tuneWithTuningParameters,     getOptimalRate, +    -- * Lift proposals+    (@~),+    liftProposal,+    liftProposalWith,+     -- * Output     proposalHeader,     summarizeProposal,@@ -56,7 +59,6 @@ import Data.Bifunctor import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL-import qualified Data.Double.Conversion.ByteString as BC import Data.Function import qualified Data.Vector as VB import qualified Data.Vector.Unboxed as VU@@ -205,37 +207,6 @@     Propose !a !KernelRatio !Jacobian   deriving (Show, Eq) --- | Lift a proposal from one data type to another.------ Assume the Jacobian is 1.0.------ For example:------ @--- scaleFirstEntryOfTuple = _1 @~ scale--- @------ See also 'liftProposalWith'.-infixl 7 @~--(@~) :: Lens' b a -> Proposal a -> Proposal b-(@~) = liftProposal---- | See '(@~)'.-liftProposal :: Lens' b a -> Proposal a -> Proposal b-liftProposal = liftProposalWith (const 1.0)---- | Lift a proposal from one data type to another.------ A function to calculate the Jacobian has to be provided. If the Jabobian is--- 1.0, use 'liftProposal'.------ 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 p w s t) =-  Proposal n r d p w (liftPFunctionWith jf l s) (liftTunerWith jf l <$> t)- -- TODO @Dominik (high, feature): Proposals should be aware of: Is this Burn in, or not?  -- | Simple proposal function without tuning information.@@ -247,23 +218,6 @@ -- proposals based on Hamiltonian dynamics). type PFunction a = a -> IOGenM StdGen -> IO (PResult a, Maybe AcceptanceCounts) --- Lift a proposal function from one data type to another.-liftPFunctionWith :: JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b-liftPFunctionWith jf l s = s'-  where-    s' y g = do-      (pr, ac) <- s (y ^. l) g-      let pr' = case pr of-            ForceAccept x' -> ForceAccept $ set l x' y-            ForceReject -> ForceReject-            Propose x' r j ->-              let y' = set l x' y-                  jxy = jf y-                  jyx = jf y'-                  j' = j * jyx / jxy-               in Propose y' r j'-      pure (pr', ac)- -- | Create a proposal with a single tuning parameter. -- -- Proposals with auxiliary tuning parameters have to be created manually. See@@ -311,13 +265,6 @@       Either String (PFunction a)   } --- 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'-  where-    fP' b d r = fP b d r . fmap (VB.map (^. l))-    g' x xs = liftPFunctionWith jf l <$> g x xs- -- | Tune proposal? data Tune = Tune | NoTune   deriving (Show, Eq)@@ -358,23 +305,23 @@ -- does not handle auxiliary tuning parameters and ignores the actual samples -- attained during the last tuning period. tuningFunction :: TuningFunction a-tuningFunction _ d r _ = bimap (tuningFunctionSimple d r) id+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) = bimap (tuningFunctionSimple d r) (f b xs)+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) = bimap id (f b xs)+tuningFunctionOnlyAux _ _ _ _ Nothing _ = error "tuningFunctionOnlyAux: empty trace"+tuningFunctionOnlyAux f b _ _ (Just xs) (!t, !ts) = bimap id (f b xs) (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@@ -434,6 +381,61 @@ getOptimalRate PDimensionUnknown = 0.234 getOptimalRate (PSpecial _ r) = r +-- | Lift a proposal from one data type to another.+--+-- Assume the Jacobian is 1.0.+--+-- For example:+--+-- @+-- scaleFirstEntryOfTuple = _1 @~ scale+-- @+--+-- See also 'liftProposalWith'.+infixl 7 @~++(@~) :: Lens' b a -> Proposal a -> Proposal b+(@~) = liftProposal++-- | See '(@~)'.+liftProposal :: Lens' b a -> Proposal a -> Proposal b+liftProposal = liftProposalWith (const 1.0)++-- | Lift a proposal from one data type to another.+--+-- A function to calculate the Jacobian has to be provided. If the Jabobian is+-- 1.0, use 'liftProposal'.+--+-- 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 p w s t) =+  Proposal n r d p w (liftPFunctionWith jf l s) (liftTunerWith jf l <$> t)++-- Lift a proposal function from one data type to another.+liftPFunctionWith :: JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b+liftPFunctionWith jf l s = s'+  where+    s' y g = do+      (pr, ac) <- s (y ^. l) g+      let pr' = case pr of+            ForceAccept x' -> ForceAccept $ set l x' y+            ForceReject -> ForceReject+            Propose x' r j ->+              let y' = set l x' y+                  jxy = jf y+                  jyx = jf y'+                  j' = j * jyx / jxy+               in Propose y' r j'+      pure (pr', ac)++-- 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'+  where+    fP' b d r = fP b d r . fmap (VB.map (^. l))+    g' x xs = liftPFunctionWith jf l <$> g x xs+ -- Warn if acceptance rate is lower. rateMin :: Double rateMin = 0.1@@ -500,12 +502,13 @@     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 = BL.fromStrict $ maybe "" (BC.toFixed 2 . (^. _3)) ar-    optimalRate = BL.fromStrict $ BC.toFixed 2 $ getOptimalRate dimension-    tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 4) tuningParameter+    acceptRate = BB.toLazyByteString $ maybe "" (fN 2 . (^. _3)) ar+    optimalRate = BB.toLazyByteString $ fN 2 $ getOptimalRate dimension+    tuneParamStr = BB.toLazyByteString $ maybe "" (fN 4) tuningParameter     checkRate rate       | rate < rateMin = Just "rate too low"       | rate > rateMax = Just "rate too high"

@@ -17,6 +17,7 @@ -- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler: --   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of --   Machine Learning Research.+{-# LANGUAGE BangPatterns #-}  -- | -- Module      :  Mcmc.Proposal.Hamiltonian.Internal@@ -259,7 +260,7 @@ hTuningFunctionWith n toVec (HTuningConf lc mc) = case (lc, mc) of   (HNoTuneLeapfrog, HNoTuneMasses) -> Nothing   (_, _) -> Just $-    \tt pdim ar mxs (_, ts) ->+    \tt pdim ar mxs (_, !ts) ->       case mxs of         Nothing -> error "hTuningFunctionWith: empty trace"         Just xs ->