mcmc 0.3.0 → 0.4.0.0

raw patch · 39 files changed

+2653/−1630 lines, 39 filesdep +circulardep +deepseqdep +monad-paralleldep −QuickCheckdep −hspec-discoverdep −math-functionsdep ~aesondep ~bytestringdep ~containersPVP ok

version bump matches the API change (PVP)

Dependencies added: circular, deepseq, monad-parallel, pretty-show

Dependencies removed: QuickCheck, hspec-discover, math-functions

Dependency ranges changed: aeson, bytestring, containers, criterion, data-default, directory, dirichlet, double-conversion, hspec, log-domain, microlens, mwc-random, primitive, statistics, time, transformers, vector, zlib

API changes (from Hackage documentation)

Files

• ChangeLog.md +17/−0
• README.md +22/−6
• bench/Bench.hs +72/−3
• bench/Normal.hs +69/−33
• bench/Poisson.hs +15/−21
• mcmc.cabal +107/−105
• src/Mcmc.hs +100/−73
• src/Mcmc/Algorithm.hs +62/−0
• src/Mcmc/Algorithm/MC3.hs +546/−0
• src/Mcmc/Algorithm/Metropolis.hs +254/−0
• src/Mcmc/Chain/Chain.hs +95/−0
• src/Mcmc/Chain/Link.hs +47/−0
• src/Mcmc/Chain/Save.hs +108/−0
• src/Mcmc/Chain/Trace.hs +75/−0
• src/Mcmc/Environment.hs +49/−0
• src/Mcmc/Internal/ByteString.hs +9/−3
• src/Mcmc/Internal/Random.hs +48/−0
• src/Mcmc/Internal/Shuffle.hs +14/−28
• src/Mcmc/Item.hs +0/−48
• src/Mcmc/Mcmc.hs +234/−222
• src/Mcmc/Metropolis.hs +0/−196
• src/Mcmc/Monitor.hs +73/−92
• src/Mcmc/Monitor/Parameter.hs +0/−14
• src/Mcmc/Monitor/ParameterBatch.hs +28/−42
• src/Mcmc/Monitor/Time.hs +7/−1
• src/Mcmc/Prior.hs +42/−28
• src/Mcmc/Proposal.hs +229/−115
• src/Mcmc/Proposal/Bactrian.hs +9/−9
• src/Mcmc/Proposal/Generic.hs +1/−1
• src/Mcmc/Proposal/Scale.hs +9/−22
• src/Mcmc/Proposal/Simplex.hs +13/−7
• src/Mcmc/Proposal/Slide.hs +13/−13
• src/Mcmc/Save.hs +0/−180
• src/Mcmc/Settings.hs +236/−0
• src/Mcmc/Status.hs +0/−218
• src/Mcmc/Trace.hs +0/−63
• src/Mcmc/Verbosity.hs +0/−40
• test/Mcmc/ProposalSpec.hs +6/−6
• test/Mcmc/SaveSpec.hs +44/−41

@@ -5,6 +5,23 @@ ## Unreleased changes  +## 0.4.0.0++-   Greatly improve documentation.+-   Major design change: Introduction of the `Algorithm` type class; algorithms+    are data types. See `MHG`.+-   Metropolic-coupled Markov chain Monte Carlo algorithm (parallel chains).+-   Optimal acceptance rate depends on dimension of proposal.+-   Use a circular trace with constant memory usage (big change).+-   Therefore, batch monitors use vectors now.+-   Always save chain with complete trace (but with sensible length).+-   Determine necessary trace length at initialization.+-   Rename `Item` to `Link`.+-   Rename `Status` to `Chain` and separate `Settings` and `Environment` from the+    `Chain`.+-   Many bug fixes.++ ## 0.3.0  -   New shorter example/test for dating trees.

@@ -1,14 +1,17 @@ -# Markov chain Monte Carlo+# Markov chain Monte Carlo sampler  <p align="center"><img src="https://travis-ci.org/dschrempf/mcmc.svg?branch=master"/></p> -Sample from a posterior using Markov chain Monte Carlo methods.+Sample from a posterior using Markov chain Monte Carlo (MCMC) algorithms. -At the moment, the library is tailored to the Metropolis-Hastings algorithm-since it covers most use cases. More algorithms will be implemented soon.+At the moment, the following algorithms are available: +-   Metropolis-Hastings-Green <sup><a id="fnr.1" class="footref" href="#fn.1">1</a></sup>;+-   Metropolis-coupled Markov chain Monte Carlo (also known as parallel+    tempering) <sup><a id="fnr.2" class="footref" href="#fn.2">2</a></sup> <sup>, </sup><sup><a id="fnr.3" class="footref" href="#fn.3">3</a></sup>. + ## Documentation  The [source code](https://hackage.haskell.org/package/mcmc) contains detailed documentation about general concepts as well@@ -17,8 +20,8 @@  ## Examples -Have a look at the [example MCMC analyses](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples). They can be built with [Stack](https://docs.haskellstack.org/en/stable/README/) and are-attached to this repository.+[Example MCMC analyses](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples) can be built with [Stack](https://docs.haskellstack.org/en/stable/README/) and are attached to this+repository.      git clone https://github.com/dschrempf/mcmc.git     cd mcmc@@ -28,3 +31,16 @@      stack exec archery ++# Footnotes++<sup><a id="fn.1" href="#fnr.1">1</a></sup> Geyer, C. J., Introduction to Markov chain Monte Carlo, In Handbook of+Markov Chain Monte Carlo (pp. 45) (2011). CRC press.++<sup><a id="fn.2" href="#fnr.2">2</a></sup> Geyer, C. J., Markov chain monte carlo maximum likelihood, Computing+Science and Statistics, Proceedings of the 23rd Symposium on the Interface,+(), (1991).++<sup><a id="fn.3" href="#fnr.3">3</a></sup> Altekar, G., Dwarkadas, S., Huelsenbeck, J. P., & Ronquist, F., Parallel+metropolis coupled markov chain monte carlo for bayesian phylogenetic inference,+Bioinformatics, 20(3), 407–415 (2004).

@@ -22,7 +22,76 @@ main = do   g <- create   defaultMain-    [ bench "Normal" $ nfIO (normalBench g),-      bench "NormalBactrian" $ nfIO (normalBactrianBench g),-      bench "Poisson" $ nfIO (poissonBench g)+    [ bgroup+        "Normal"+        [ bench "Slide" $ nfIO (normalSlideBench g),+          bench "Bactrian" $ nfIO (normalBactrianBench g),+          bench "LargeCycle" $ nfIO (normalLargeCycleBench g)+        ],+      bench "Poisson" $ nfIO (poissonBench g),+      bgroup+        "MC3"+        [ bench "MC3 2" $ nfIO (normalMC3 g 2),+          bench "MC3 3" $ nfIO (normalMC3 g 3),+          bench "MC3 4" $ nfIO (normalMC3 g 4),+          bench "MC3 5" $ nfIO (normalMC3 g 5),+          bench "MC3 10" $ nfIO (normalMC3 g 10)+        ]     ]++-- benchmarking Normal/Slide+-- time                 42.31 ms   (41.88 ms .. 42.60 ms)+--                      1.000 R²   (0.999 R² .. 1.000 R²)+-- mean                 42.75 ms   (42.52 ms .. 43.29 ms)+-- std dev              661.4 μs   (347.4 μs .. 1.074 ms)++-- benchmarking Normal/Bactrian+-- time                 45.51 ms   (45.30 ms .. 45.92 ms)+--                      1.000 R²   (0.999 R² .. 1.000 R²)+-- mean                 45.41 ms   (45.31 ms .. 45.61 ms)+-- std dev              276.0 μs   (141.5 μs .. 460.8 μs)++-- benchmarking Normal/LargeCycle+-- time                 68.82 ms   (67.18 ms .. 70.81 ms)+--                      0.999 R²   (0.997 R² .. 1.000 R²)+-- mean                 67.68 ms   (67.26 ms .. 68.59 ms)+-- std dev              1.074 ms   (618.2 μs .. 1.602 ms)++-- benchmarking Poisson+-- time                 72.94 ms   (63.07 ms .. 87.73 ms)+--                      0.953 R²   (0.920 R² .. 1.000 R²)+-- mean                 64.76 ms   (62.84 ms .. 71.62 ms)+-- std dev              5.785 ms   (783.9 μs .. 10.08 ms)+-- variance introduced by outliers: 26% (moderately inflated)++-- benchmarking MC3/MC3 2+-- time                 13.08 ms   (12.73 ms .. 13.44 ms)+--                      0.993 R²   (0.986 R² .. 0.997 R²)+-- mean                 13.41 ms   (13.16 ms .. 13.72 ms)+-- std dev              682.7 μs   (520.3 μs .. 874.0 μs)+-- variance introduced by outliers: 22% (moderately inflated)++-- benchmarking MC3/MC3 3+-- time                 19.19 ms   (18.86 ms .. 19.59 ms)+--                      0.998 R²   (0.996 R² .. 1.000 R²)+-- mean                 19.28 ms   (19.11 ms .. 19.51 ms)+-- std dev              454.1 μs   (339.1 μs .. 608.4 μs)++-- benchmarking MC3/MC3 4+-- time                 25.01 ms   (24.21 ms .. 25.66 ms)+--                      0.997 R²   (0.996 R² .. 0.999 R²)+-- mean                 24.21 ms   (23.99 ms .. 24.55 ms)+-- std dev              606.1 μs   (414.7 μs .. 738.5 μs)++-- benchmarking MC3/MC3 5+-- time                 28.39 ms   (26.99 ms .. 29.34 ms)+--                      0.995 R²   (0.990 R² .. 0.999 R²)+-- mean                 31.13 ms   (30.09 ms .. 33.49 ms)+-- std dev              3.009 ms   (984.0 μs .. 4.319 ms)+-- variance introduced by outliers: 40% (moderately inflated)++-- benchmarking MC3/MC3 10+-- time                 57.25 ms   (56.98 ms .. 57.61 ms)+--                      1.000 R²   (1.000 R² .. 1.000 R²)+-- mean                 57.46 ms   (57.34 ms .. 57.56 ms)+-- std dev              192.4 μs   (140.1 μs .. 284.0 μs)

@@ -1,6 +1,6 @@ -- | -- Module      :  Normal--- Description :  Benchmark Metropolis-Hastings algorithm+-- Description :  Benchmark Metropolis-Hastings-Green algorithm -- Copyright   :  (c) Dominik Schrempf 2020 -- License     :  GPL-3.0-or-later --@@ -10,34 +10,28 @@ -- -- Creation date: Wed May  6 00:10:11 2020. module Normal-  ( normalBench,+  ( normalSlideBench,     normalBactrianBench,+    normalLargeCycleBench,+    normalMC3,   ) where  import Control.Monad import Mcmc-import Numeric.Log as L-import Statistics.Distribution hiding-  ( mean,-    stdDev,-  )-import Statistics.Distribution.Normal import System.Random.MWC  trueMean :: Double trueMean = 5 -trueStdDev :: Double-trueStdDev = 4+stdDev :: Double+stdDev = 4 -lh :: Double -> Log Double-lh = Exp . logDensity (normalDistr trueMean trueStdDev)+lh :: LikelihoodFunction Double+lh = normal trueMean stdDev -proposals :: Cycle Double-proposals =-  fromList-    [slideSymmetric 1.0 (PName "Medium") (PWeight 1) Tune]+cc :: Cycle Double+cc = cycleFromList [slideSymmetric 1.0 (PName "Medium") (PWeight 1) Tune]  mons :: [MonitorParameter Double] mons = [monitorDouble "mu"]@@ -48,26 +42,68 @@ mon :: Monitor Double mon = Monitor monStd [] [] -nBurn :: Maybe Int-nBurn = Just 2000--nAutoTune :: Maybe Int-nAutoTune = Just 200+normalSlideBench :: GenIO -> IO ()+normalSlideBench g = do+  let s =+        Settings+          (AnalysisName "Normal")+          (BurnInWithAutoTuning 2000 500)+          (Iterations 20000)+          Overwrite+          Sequential+          NoSave+          Quiet+  a <- mhg noPrior lh cc mon 0 g+  void $ mcmc s a -nIter :: Int-nIter = 20000+ccLarge :: Cycle Double+ccLarge =+  cycleFromList+    [slideSymmetric 1.0 (PName $ "Medium " ++ show i) (PWeight 1) Tune | i <- [0 .. 100 :: Int]] -normalBench :: GenIO -> IO ()-normalBench g = do-  let s = quiet $ status "Normal" (const 1) lh proposals mon 0 nBurn nAutoTune nIter g-  void $ mh s+-- Should have the same run time as 'normalSlide'.+normalLargeCycleBench :: GenIO -> IO ()+normalLargeCycleBench g = do+  let s =+        Settings+          (AnalysisName "Normal")+          (BurnInWithAutoTuning 20 5)+          (Iterations 200)+          Overwrite+          Sequential+          NoSave+          Quiet+  a <- mhg noPrior lh ccLarge mon 0 g+  void $ mcmc s a -proposalsBactrian :: Cycle Double-proposalsBactrian =-  fromList-    [slideBactrian 0.5 1.0 (PName "Bactrian") (PWeight 1) Tune]+ccBactrian :: Cycle Double+ccBactrian = cycleFromList [slideBactrian 0.5 1.0 (PName "Bactrian") (PWeight 1) Tune]  normalBactrianBench :: GenIO -> IO () normalBactrianBench g = do-  let s = quiet $ status "NormalBactrian" (const 1) lh proposalsBactrian mon 0 nBurn nAutoTune nIter g-  void $ mh s+  let s =+        Settings+          (AnalysisName "NormalBactrian")+          (BurnInWithAutoTuning 2000 200)+          (Iterations 20000)+          Overwrite+          Sequential+          NoSave+          Quiet+  a <- mhg noPrior lh ccBactrian mon 0 g+  void $ mcmc s a++normalMC3 :: GenIO -> Int -> IO ()+normalMC3 g n = do+  let mcmcS =+        Settings+          (AnalysisName "MC3")+          (BurnInWithAutoTuning 200 20)+          (Iterations 2000)+          Overwrite+          Sequential+          NoSave+          Quiet+      mc3S = MC3Settings (NChains n) (SwapPeriod 2) (NSwaps 1)+  a <- mc3 mc3S noPrior lh cc mon 0 g+  void $ mcmc mcmcS a

@@ -20,11 +20,6 @@ import Lens.Micro import Mcmc import Numeric.Log hiding (sum)-import Statistics.Distribution hiding-  ( mean,-    stdDev,-  )-import Statistics.Distribution.Poisson import System.Random.MWC  type I = (Double, Double)@@ -39,13 +34,12 @@     m = sum ys / fromIntegral (length ys)  f :: Int -> Double -> I -> Log Double-f ft yr (a, b) = Exp $ logProbability (poisson l) (fromIntegral ft)+f ft yr (a, b) = poisson l (fromIntegral ft)   where     l = exp $ a + b * yr -lh :: I -> Log Double-lh x =-  product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]+lh :: LikelihoodFunction I+lh x = product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]  proposalAlpha :: Proposal I proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (PWeight 1) NoTune@@ -54,7 +48,7 @@ proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (PWeight 1) NoTune  proposals :: Cycle I-proposals = fromList [proposalAlpha, proposalBeta]+proposals = cycleFromList [proposalAlpha, proposalBeta]  initial :: I initial = (0, 0)@@ -71,16 +65,16 @@ mon :: Monitor I mon = Monitor monStd [] [] -nBurn :: Maybe Int-nBurn = Just 2000--nAutoTune :: Maybe Int-nAutoTune = Just 200--nIter :: Int-nIter = 10000- poissonBench :: GenIO -> IO () poissonBench g = do-  let s = quiet $ status "Poisson" (const 1) lh proposals mon initial nBurn nAutoTune nIter g-  void $ mh s+  let s =+        Settings+          (AnalysisName "Poisson")+          (BurnInWithAutoTuning 2000 200)+          (Iterations 10000)+          Overwrite+          Sequential+          NoSave+          Quiet+  a <- mhg noPrior lh proposals mon initial g+  void $ mcmc s a

@@ -1,117 +1,119 @@-cabal-version:  2.2-name:           mcmc-version:        0.3.0-synopsis:       Sample from a posterior using Markov chain Monte Carlo-description:    Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme>-category:       Math, Statistics-homepage:       https://github.com/dschrempf/mcmc#readme-bug-reports:    https://github.com/dschrempf/mcmc/issues-author:         Dominik Schrempf-maintainer:     dominik.schrempf@gmail.com-copyright:      Dominik Schrempf (2020)-license:        GPL-3.0-or-later-build-type:     Simple+cabal-version:      2.2+name:               mcmc+version:            0.4.0.0+license:            GPL-3.0-or-later+copyright:          Dominik Schrempf (2020)+maintainer:         dominik.schrempf@gmail.com+author:             Dominik Schrempf+homepage:           https://github.com/dschrempf/mcmc#readme+bug-reports:        https://github.com/dschrempf/mcmc/issues+synopsis:           Sample from a posterior using Markov chain Monte Carlo+description:+    Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme> +category:           Math, Statistics+build-type:         Simple extra-source-files:     README.md     ChangeLog.md  source-repository head-  type: git-  location: https://github.com/dschrempf/mcmc+    type:     git+    location: https://github.com/dschrempf/mcmc  library-  exposed-modules:-      Mcmc-      Mcmc.Item-      Mcmc.Mcmc-      Mcmc.Metropolis-      Mcmc.Monitor-      Mcmc.Monitor.Log-      Mcmc.Monitor.Parameter-      Mcmc.Monitor.ParameterBatch-      Mcmc.Monitor.Time-      Mcmc.Prior-      Mcmc.Proposal-      Mcmc.Proposal.Bactrian-      Mcmc.Proposal.Generic-      Mcmc.Proposal.Scale-      Mcmc.Proposal.Slide-      Mcmc.Proposal.Simplex-      Mcmc.Save-      Mcmc.Status-      Mcmc.Trace-      Mcmc.Verbosity-  other-modules:-      Mcmc.Internal.ByteString-      Mcmc.Internal.Shuffle-      Paths_mcmc-  autogen-modules:-      Paths_mcmc-  hs-source-dirs:-      src-  ghc-options: -Wall-  build-depends:-      aeson-    , base >=4.7 && <5-    , bytestring-    , containers-    , data-default-    , directory-    , dirichlet-    , double-conversion-    , log-domain-    , math-functions-    , microlens-    , mwc-random-    , primitive-    , statistics-    , time-    , transformers-    , vector-    , zlib-  default-language: Haskell2010+    exposed-modules:+        Mcmc+        Mcmc.Algorithm+        Mcmc.Algorithm.Metropolis+        Mcmc.Algorithm.MC3+        Mcmc.Chain.Chain+        Mcmc.Chain.Link+        Mcmc.Chain.Save+        Mcmc.Chain.Trace+        Mcmc.Environment+        Mcmc.Mcmc+        Mcmc.Monitor+        Mcmc.Monitor.Log+        Mcmc.Monitor.Parameter+        Mcmc.Monitor.ParameterBatch+        Mcmc.Monitor.Time+        Mcmc.Prior+        Mcmc.Proposal+        Mcmc.Proposal.Bactrian+        Mcmc.Proposal.Generic+        Mcmc.Proposal.Scale+        Mcmc.Proposal.Slide+        Mcmc.Proposal.Simplex+        Mcmc.Settings +    hs-source-dirs:   src+    other-modules:+        Mcmc.Internal.ByteString+        Mcmc.Internal.Random+        Mcmc.Internal.Shuffle+        Paths_mcmc++    autogen-modules:  Paths_mcmc+    default-language: Haskell2010+    ghc-options:      -Wall -Wunused-packages+    build-depends:+        aeson >=1.5.4.1 && <1.6,+        base >=4.7 && <5,+        bytestring >=0.10.10.0 && <0.11,+        circular >=0.3.1.1 && <0.4,+        containers >=0.6.2.1 && <0.7,+        data-default >=0.7.1.1 && <0.8,+        deepseq >=1.4.4.0 && <1.5,+        directory >=1.3.6.0 && <1.4,+        dirichlet >=0.1.0.0 && <0.2,+        double-conversion >=2.0.2.0 && <2.1,+        log-domain ==0.13.*,+        microlens >=0.4.11.2 && <0.5,+        mwc-random >=0.14.0.0 && <0.15,+        monad-parallel >=0.7.2.3 && <0.8,+        pretty-show ==1.10.*,+        primitive >=0.7.1.0 && <0.8,+        statistics >=0.15.2.0 && <0.16,+        time >=1.9.3 && <1.10,+        transformers >=0.5.6.2 && <0.6,+        vector >=0.12.1.2 && <0.13,+        zlib >=0.6.2.2 && <0.7+ test-suite mcmc-test-  type: exitcode-stdio-1.0-  main-is: Spec.hs-  other-modules:-      Mcmc.ProposalSpec-      Mcmc.SaveSpec-      Paths_mcmc-  hs-source-dirs:-      test-  ghc-options: -Wall-  build-depends:-      QuickCheck-    , base >=4.7 && <5-    , directory-    , hspec-    , hspec-discover-    , log-domain-    , mcmc-    , mwc-random-    , statistics-    , vector-  default-language: Haskell2010+    type:             exitcode-stdio-1.0+    main-is:          Spec.hs+    hs-source-dirs:   test+    other-modules:+        Mcmc.ProposalSpec+        Mcmc.SaveSpec+        Paths_mcmc +    default-language: Haskell2010+    ghc-options:      -Wall -Wunused-packages+    build-depends:+        base >=4.7 && <5,+        hspec >=2.7.4 && <2.8,+        log-domain ==0.13.*,+        mcmc -any,+        mwc-random >=0.14.0.0 && <0.15,+        statistics >=0.15.2.0 && <0.16+ benchmark mcmc-bench-  type: exitcode-stdio-1.0-  main-is: Bench.hs-  other-modules:-      Normal-      Poisson-      Paths_mcmc-  hs-source-dirs:-      bench-  ghc-options: -Wall-  build-depends:-      base >=4.7 && <5-    , criterion-    , log-domain-    , mcmc-    , microlens-    , mwc-random-    , statistics-    , vector-  default-language: Haskell2010+    type:             exitcode-stdio-1.0+    main-is:          Bench.hs+    hs-source-dirs:   bench+    other-modules:+        Normal+        Poisson+        Paths_mcmc++    default-language: Haskell2010+    ghc-options:      -Wall -Wunused-packages+    build-depends:+        base >=4.7 && <5,+        criterion >=1.5.7.0 && <1.6,+        log-domain ==0.13.*,+        mcmc -any,+        microlens >=0.4.11.2 && <0.5,+        mwc-random >=0.14.0.0 && <0.15

@@ -2,7 +2,7 @@  -- | -- Module      :  Mcmc--- Description :  Markov chain Monte Carlo algorithms, batteries included+-- Description :  Markov chain Monte Carlo samplers, batteries included -- Copyright   :  (c) Dominik Schrempf 2020 -- License     :  GPL-3.0-or-later --@@ -12,37 +12,64 @@ -- -- Creation date: Tue May  5 18:01:15 2020. ----- A short introduction to update mechanisms using the Metropolis-Hastings--- algorithm (see Geyer, C. J., 2011; Introduction to Markov Chain Monte Carlo. In--- Handbook of Markov Chain Monte Carlo (pp. 45), Chapman \& Hall/CRC).+-- For an introduction to Markov chain Monte Carlo (MCMC) samplers and update+-- mechanisms using the Metropolis-Hastings-Green algorithm, please see Geyer,+-- C. J., (2011), Introduction to Markov Chain Monte Carlo, In Handbook of+-- Markov Chain Monte Carlo (pp. 45), CRC press. ----- For examples, please see--- [mcmc-examples](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples).+-- This library focusses on classical Markov chain Monte Carlo algorithms such+-- as the Metropolis-Hastings-Green [1] algorithm, or population methods+-- involving parallel chains such as the Metropolic-coupled Markov chain Monte+-- Carlo [2] algorithm. In particular, sequential Monte Carlo [3] algorithms+-- following a moving posterior distribution are not provided. --+-- An MCMC sampler can be run with 'mcmc', for example using the+-- Metropolis-Hastings-Green algorithm 'mhg'.+--+-- Usually, it is best to start with an example:+--+-- - Basic inference of the [accuracy of an+--   archer](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples/Archery/Archery.hs)+--+-- - [More involved+--   examples](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples/Archery/Archery.hs)+-- -- __The import of this module alone should cover most use cases.__+--+-- @[1]@ Geyer, C. J. (2011), Introduction to markov chain monte carlo, In+-- Handbook of Markov Chain Monte Carlo (pp. 45), CRC press.+--+-- @[2]@ Geyer, C. J. (1991), Markov chain monte carlo maximum likelihood,+-- Computing Science and Statistics, Proceedings of the 23rd Symposium on the+-- Interface.+--+-- @[3]@ Sequential monte carlo methods in practice (2001), Editors: Arnaud+-- Doucet, Nando de Freitas, and Neil Gordon, Springer New York. module Mcmc   ( -- * Proposals -    -- | A 'Proposal' is an instruction about how to advance a given Markov chain so-    -- that it possibly reaches a new state. That is, 'Proposal's specify how the-    -- chain traverses the state space. As far as this MCMC library is-    -- concerned, 'Proposal's are /elementary updates/ in that they cannot be-    -- decomposed into smaller updates.+    -- | A 'Proposal' is an instruction about how to advance a given Markov+    -- chain so that it possibly reaches a new state. That is, 'Proposal's+    -- specify how the chain traverses the state space. As far as this MCMC+    -- library is concerned, 'Proposal's are considered to be /elementary+    -- updates/ in that they cannot be decomposed into smaller updates.     --     -- 'Proposal's can be combined to form composite updates, a technique often     -- referred to as /composition/. On the other hand, /mixing/ (used in the     -- sense of mixture models) is the random choice of a 'Proposal' (or a     -- composition of 'Proposal's) from a given set.     ---    -- The __composition__ and __mixture__ of 'Proposal's allows specification of-    -- nearly all MCMC algorithms involving a single chain (i.e., population+    -- The __composition__ and __mixture__ of 'Proposal's allows specification+    -- of nearly all MCMC algorithms involving a single chain (i.e., population     -- methods such as particle filters are excluded). In particular, Gibbs-    -- samplers of all sorts can be specified using this procedure.+    -- samplers of all sorts can be specified using this procedure. For+    -- reference, please see the short [encyclopedia of MCMC+    -- methods](https://dschrempf.github.io/coding/2020-11-12-encyclopedia-of-markov-chain-monte-carlo-methods/).     ---    -- This library enables composition and mixture of 'Proposal's via the 'Cycle'-    -- data type. Essentially, a 'Cycle' is a set of 'Proposal's. The chain advances-    -- after the completion of each 'Cycle', which is called an __iteration__,-    -- and the iteration counter is increased by one.+    -- This library enables composition and mixture of 'Proposal's via the+    -- 'Cycle' data type. Essentially, a 'Cycle' is a set of 'Proposal's. The+    -- chain advances after the completion of each 'Cycle', which is called an+    -- __iteration__, and the iteration counter is increased by one.     --     -- The 'Proposal's in a 'Cycle' can be executed in the given order or in a     -- random sequence which allows, for example, specification of a fixed scan@@ -56,32 +83,35 @@     -- Proposals are named according to what they do, i.e., how they change the     -- state of a Markov chain, and not according to the intrinsically used     -- probability distributions. For example, 'slideSymmetric' is a sliding-    -- proposal. Under the hood, it uses the normal distribution with mean zero and-    -- given variance. The sampled variate is added to the current value of the-    -- variable (hence, the name slide). The same nomenclature is used by-    -- RevBayes [1]. The probability distributions and intrinsic properties of a+    -- proposal. Under the hood, it uses the normal distribution with mean zero+    -- and given variance. The sampled variate is added to the current value of+    -- the variable (hence, the name slide). The same nomenclature is used by+    -- RevBayes [4]. The probability distributions and intrinsic properties of a     -- specific proposal are specified in detail in the documentation.     --     -- The other method, which is used intrinsically, is more systematic, but     -- also a little bit more complicated: we separate between the proposal-    -- distribution and how the state is affected. And here, I am not only-    -- referring to the accessor (i.e., the lens), but also to the operator-    -- (addition, multiplication, any other binary operator). For example, the-    -- sliding proposal (without tuning information) is implemented as+    -- distribution and how the state is affected. And here, I am referring to+    -- the operator (addition, multiplication, any other binary operator). For+    -- example, the sliding proposal with mean @m@, standard deviation @s@, and+    -- tuning parameter @t@ is implemented as     --     -- @-    -- slideSimple :: Lens' a Double -> Double -> Double -> Double -> ProposalSimple a-    -- slideSimple l m s t = genericContinuous l (normalDistr m (s * t)) (+) (-)+    -- slideSimple :: Double -> Double -> Double -> ProposalSimple Double+    -- slideSimple m s t =+    --   genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing     -- @     --     -- This specification is more involved. Especially since we need to know the-    -- probability of jumping back, and so we need to know the inverse operator.-    -- However, it also allows specification of new proposals with great ease.+    -- probability of jumping back, and so we need to know the inverse operator+    -- 'negate'. However, it also allows specification of new proposals with+    -- great ease.     ---    -- [1] Höhna, S., Landis, M. J., Heath, T. A., Boussau, B., Lartillot, N., Moore,-    -- B. R., Huelsenbeck, J. P., …, Revbayes: bayesian phylogenetic inference using-    -- graphical models and an interactive model-specification language, Systematic-    -- Biology, 65(4), 726–736 (2016). http://dx.doi.org/10.1093/sysbio/syw021+    -- @[4]@ Höhna, S., Landis, M. J., Heath, T. A., Boussau, B., Lartillot, N.,+    -- Moore, B. R., Huelsenbeck, J. P., …, Revbayes: bayesian phylogenetic+    -- inference using graphical models and an interactive model-specification+    -- language, Systematic Biology, 65(4), 726–736 (2016).+    -- http://dx.doi.org/10.1093/sysbio/syw021     PName (..),     PWeight (..),     Proposal,@@ -98,41 +128,28 @@     slideBactrian,     module Mcmc.Proposal.Simplex,     Cycle,-    fromList,+    cycleFromList,     Order (..),     setOrder, -    -- * Initialization--    -- | The 'Status' contains all information to run an MCMC chain. It is-    -- constructed using the function 'status'.-    ---    -- The 'Status' of a Markov chain includes information about current state-    -- ('Mcmc.Item.Item') and iteration, the history of the chain-    -- ('Mcmc.Trace.Trace'), the 'Acceptance' ratios, and the random number-    -- generator.-    ---    -- Further, the 'Status' includes auxiliary variables and functions such as-    -- the prior and likelihood functions, instructions to move around the state-    -- space (see above) and to monitor the MCMC run, as well as some auxiliary-    -- information.-    status,-    Cleaner (..),-    cleanWith,-    saveWith,-    force,-    quiet,-    debug,-    noData,+    -- * Settings+    module Mcmc.Settings, -    -- * Monitor+    -- * Monitors      -- | A 'Monitor' describes which part of the Markov chain should be logged-    -- and where. There are three different types:-    -- - 'MonitorStdOut': Log to standard output.-    -- - 'MonitorFile': Log to a file.-    -- - 'MonitorBatch': Log summary statistics such as the mean of the last-    -- - states to a file.+    -- and where. Monitor files can directly be loaded into established MCMC+    -- analysis tools working with tab separated tables (for example,+    -- [Tracer](http://tree.bio.ed.ac.uk/software/tracer/)).+    --+    -- There are three different 'Monitor' types:+    --+    -- ['MonitorStdOut'] Log to standard output.+    --+    -- ['MonitorFile'] Log to a file.+    --+    -- ['MonitorBatch'] Log summary statistics such as the mean of the last+    -- states to a file.     Monitor (Monitor),     MonitorStdOut,     monitorStdOut,@@ -144,22 +161,33 @@     module Mcmc.Monitor.ParameterBatch,      -- * Prior distributions++    -- | Convenience functions for computing priors.     module Mcmc.Prior, +    -- * MCMC samplers+    mcmc,+    mcmcContinue,++    -- | See also 'settingsLoad', 'mhgLoad', and 'mc3Load'.+     -- * Algorithms+    module Mcmc.Algorithm.Metropolis,+    module Mcmc.Algorithm.MC3, -    -- | At the moment, the library is tailored to the Metropolis-Hastings-    -- algorithm ('mh') since it covers most use cases. However, implementation-    -- of more algorithms is planned in the future.-    mh,-    mhContinue, -    -- * Misc-    loadStatus,+    -- * Useful type synonyms+    PriorFunction,+    noPrior,+    LikelihoodFunction,+    noLikelihood,   ) where -import Mcmc.Metropolis+import Mcmc.Algorithm.MC3+import Mcmc.Algorithm.Metropolis+import Mcmc.Chain.Chain+import Mcmc.Mcmc import Mcmc.Monitor import Mcmc.Monitor.Parameter import Mcmc.Monitor.ParameterBatch@@ -169,5 +197,4 @@ import Mcmc.Proposal.Scale import Mcmc.Proposal.Simplex import Mcmc.Proposal.Slide-import Mcmc.Save-import Mcmc.Status+import Mcmc.Settings

@@ -0,0 +1,62 @@+-- |+-- Module      :  Mcmc.Algorithm+-- Description :  Algortihms for Markov chain Monte Carlo samplers+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Mon Nov 16 14:37:11 2020.+module Mcmc.Algorithm+  ( Algorithm (..),+  )+where++import qualified Data.ByteString.Lazy.Char8 as BL+import Data.Time+import Mcmc.Settings++-- | Class for algorithms used by MCMC samplers.+class Algorithm a where+  -- | Name.+  aName :: a -> String++  -- | Current iteration.+  aIteration :: a -> Int++  -- | Sample the next state.+  aIterate :: ParallelizationMode -> a -> IO a++  -- | Auto tune all proposals.+  aAutoTune :: a -> a++  -- | Reset acceptance counts.+  aResetAcceptance :: a -> a++  -- | Summarize the cycle.+  aSummarizeCycle :: a -> BL.ByteString++  -- | Open required monitor files and setup corresponding file handles.+  aOpenMonitors :: AnalysisName -> ExecutionMode -> a -> IO a++  -- | Execute file monitors and possibly return a monitor string to be written+  -- to the standard output and the log file.+  aExecuteMonitors ::+    Verbosity ->+    -- | Starting time.+    UTCTime ->+    -- | Total number of iterations including burn in.+    Int ->+    a ->+    IO (Maybe BL.ByteString)++  -- | Header of monitor to standard output.+  aStdMonitorHeader :: a -> BL.ByteString++  -- | Close monitor files and remove the file handles.+  aCloseMonitors :: a -> IO a++  -- | Save analysis.+  aSave :: AnalysisName -> a -> IO ()

@@ -0,0 +1,546 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TemplateHaskell #-}++-- |+-- Module      :  Mcmc.Algorithm.MC3+-- Description :  Metropolis-coupled Markov chain Monte Carlo algorithm+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Mon Nov 23 15:20:33 2020.+--+-- The Metropolis-coupled Markov chain Monte Carlo ('MC3') algorithm.+--+-- Also known as parallel tempering.+--+-- Like any other parallel MCMC algorithm, the 'MC3' algorithm is essentially an+-- 'Mcmc.Algorithm.Metropolis.MHG' algorithm on the product space of all+-- parallel chains.+--+-- For example, see+--+-- - Geyer, C. J., Markov chain monte carlo maximum likelihood, Computing+--   Science and Statistics, Proceedings of the 23rd Symposium on the Interface,+--   (1991).+--+-- - Altekar, G., Dwarkadas, S., Huelsenbeck, J. P., & Ronquist, F., Parallel+--   metropolis coupled markov chain monte carlo for bayesian phylogenetic+--   inference, Bioinformatics, 20(3), 407–415 (2004).+module Mcmc.Algorithm.MC3+  ( -- * Definitions+    NChains (..),+    SwapPeriod (..),+    NSwaps (..),+    MC3Settings (..),+    MHGChains,+    ReciprocalTemperatures,++    -- * Metropolis-coupled Markov chain Monte Carlo algorithm+    MC3 (..),+    mc3,+    mc3Save,+    mc3Load,+  )+where++import Codec.Compression.GZip+import Control.Monad+import qualified Control.Monad.Parallel as P+import Data.Aeson+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+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U+import Data.Word+-- import Debug.Trace hiding (trace)+import Mcmc.Algorithm+import Mcmc.Algorithm.Metropolis+import Mcmc.Chain.Chain+import Mcmc.Chain.Link+import Mcmc.Chain.Save+import Mcmc.Chain.Trace+import Mcmc.Internal.Random+import Mcmc.Internal.Shuffle+import Mcmc.Monitor+import Mcmc.Proposal+import Mcmc.Settings+import Numeric.Log hiding (sum)+import System.Random.MWC++-- | Total number of parallel chains.+--+-- Must be two or larger.+newtype NChains = NChains {fromNChains :: Int}+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''NChains)++-- | The period of proposing state swaps between chains.+--+-- Must be one or larger.+newtype SwapPeriod = SwapPeriod {fromSwapPeriod :: Int}+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''SwapPeriod)++-- | The number of proposed swaps at each swapping event.+--+-- Must be in @[1, NChains - 1]@.+newtype NSwaps = NSwaps {fromNSwaps :: Int}+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''NSwaps)++-- | MC3 settings.+data MC3Settings = MC3Settings+  { -- | The number of chains has to be larger equal two.+    mc3NChains :: NChains,+    mc3SwapPeriod :: SwapPeriod,+    mc3NSwaps :: NSwaps+  }+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''MC3Settings)++-- | Vector of MHG chains.+type MHGChains a = V.Vector (MHG a)++-- | Vector of reciprocal temperatures.+type ReciprocalTemperatures = U.Vector Double++data SavedMC3 a = SavedMC3+  { savedMC3Settings :: MC3Settings,+    savedMC3Chains :: V.Vector (SavedChain a),+    savedMC3ReciprocalTemperatures :: ReciprocalTemperatures,+    savedMC3Iteration :: Int,+    savedMC3SwapAcceptance :: Acceptance Int,+    savedMC3Generator :: U.Vector Word32+  }+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''SavedMC3)++toSavedMC3 ::+  MC3 a ->+  IO (SavedMC3 a)+toSavedMC3 (MC3 s mhgs bs i ac g) = do+  scs <- V.mapM (toSavedChain . fromMHG) mhgs+  g' <- saveGen g+  return $ SavedMC3 s scs bs i ac g'++fromSavedMC3 ::+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  SavedMC3 a ->+  IO (MC3 a)+fromSavedMC3 pr lh cc mn (SavedMC3 s scs bs i ac g') = do+  mhgs <- V.mapM (fmap MHG . fromSavedChain pr lh cc mn) scs+  g <- loadGen g'+  return $ MC3 s mhgs bs i ac g++-- | The MC3 algorithm.+data MC3 a = MC3+  { mc3Settings :: MC3Settings,+    -- | The first chain is the cold chain with temperature 1.0.+    mc3MHGChains :: MHGChains a,+    -- | Vector of reciprocal temperatures.+    mc3ReciprocalTemperatures :: ReciprocalTemperatures,+    -- | Current iteration.+    mc3Iteration :: Int,+    -- | Number of accepted and rejected swaps.+    mc3SwapAcceptance :: Acceptance Int,+    mc3Generator :: GenIO+  }++instance ToJSON a => Algorithm (MC3 a) where+  aName = const "Metropolis-coupled Markov chain Monte Carlo (MC3)"+  aIteration = mc3Iteration+  aIterate = mc3Iterate+  aAutoTune = mc3AutoTune+  aResetAcceptance = mc3ResetAcceptance+  aSummarizeCycle = mc3SummarizeCycle+  aOpenMonitors = mc3OpenMonitors+  aExecuteMonitors = mc3ExecuteMonitors+  aStdMonitorHeader = mc3StdMonitorHeader+  aCloseMonitors = mc3CloseMonitors+  aSave = mc3Save++--  The prior and likelihood values of the current link are updated.+--+-- NOTE: The trace is not changed! In particular, the prior and likelihood+-- values are not updated for any link of the trace, and no new link is added to+-- the trace.+setReciprocalTemperature ::+  -- Cold prior function.+  PriorFunction a ->+  -- Cold likelihood function.+  LikelihoodFunction a ->+  -- New reciprocal temperature.+  Double ->+  MHG a ->+  MHG a+setReciprocalTemperature prf lhf beta a =+  MHG $+    c+      { priorFunction = prf',+        likelihoodFunction = lhf',+        link = Link x (prf' x) (lhf' x)+      }+  where+    c = fromMHG a+    b' = Exp $ log beta+    -- We need twice the amount of computations compared to taking the power+    -- after calculating the posterior (pr x * lh x) ** b'. However, I don't+    -- think this is a serious problem.+    --+    -- To minimize computations, it is key to avoid modification of the+    -- reciprocal temperature for the cold chain.+    prf' = (** b') . prf+    lhf' = (** b') . lhf+    x = state $ link c++initMHG ::+  -- Cold prior function.+  PriorFunction a ->+  -- Cold likelihood function.+  LikelihoodFunction a ->+  -- Index of MHG chain.+  Int ->+  -- Reciprocal temperature.+  Double ->+  MHG a ->+  IO (MHG a)+initMHG prf lhf i beta a+  | i < 0 = error "initMHG: Chain index negative."+  -- Do not temper with the cold chain.+  | i == 0 = return a+  | 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 = i, trace = t'}+  where+    a' = setReciprocalTemperature prf lhf beta a+    c = fromMHG a'+    l = link c+    t = trace c++-- TODO: Splitmix. Initialization of the MC3 algorithm is an IO action because+-- the generators have to be split. And also because of the mutable trace.++-- | Initialize an MC3 algorithm with a given number of chains.+--+-- Call 'error' if:+--+-- - The number of chains is one or lower.+--+-- - The swap period is zero or negative.+mc3 ::+  MC3Settings ->+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  a ->+  GenIO ->+  IO (MC3 a)+mc3 s pr lh cc mn i0 g+  | n < 2 = error "mc3: The number of chains must be two or larger."+  | sp < 1 = error "mc3: The swap period must be strictly positive."+  | sn < 1 || sn > n - 1 = error "mc3: The number of swaps must be in [1, NChains - 1]."+  | otherwise = do+    -- Split random number generators.+    gs <- V.fromList <$> splitGen n g+    cs <- V.mapM (mhg pr lh cc mn i0) gs+    hcs <- V.izipWithM (initMHG pr lh) (V.convert bs) cs+    return $ MC3 s hcs bs 0 (emptyA [0 .. n - 2]) g+  where+    n = fromNChains $ mc3NChains s+    sp = fromSwapPeriod $ mc3SwapPeriod s+    sn = fromNSwaps $ mc3NSwaps s+    -- NOTE: The initial choice of reciprocal temperatures is based on a few+    -- tests but otherwise pretty arbitrary.+    --+    -- NOTE: Have to 'take n' elements, because vectors are not as lazy as+    -- lists.+    bs = U.fromList $ take n $ iterate (* 0.92) 1.0++mc3Fn :: AnalysisName -> FilePath+mc3Fn (AnalysisName nm) = nm ++ ".mc3"++-- | Save an MC3 algorithm.+mc3Save ::+  ToJSON a =>+  AnalysisName ->+  MC3 a ->+  IO ()+mc3Save nm a = do+  savedMC3 <- toSavedMC3 a+  BL.writeFile (mc3Fn nm) $ compress $ encode savedMC3++-- | Load an MC3 algorithm.+--+-- See 'Mcmc.Mcmc.mcmcContinue'.+mc3Load ::+  FromJSON a =>+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  AnalysisName ->+  IO (MC3 a)+mc3Load pr lh cc mn nm = do+  savedMC3 <- eitherDecode . decompress <$> BL.readFile (mc3Fn nm)+  either error (fromSavedMC3 pr lh cc mn) savedMC3++-- 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+-- than the right chain.+swapWith ::+  -- Index i>=0 of left chain.+  Int ->+  -- Index j>=0, j/=i of right chain.+  Int ->+  MHGChains a ->+  (MHGChains a, Log Double)+swapWith i j xs+  | i < 0 = error "swapWith: Left index is negative."+  | j < 0 = error "swapWith: Right index is negative."+  | i == j = error "swapWith: Indices are equal."+  | otherwise = (xs', q)+  where+    -- Gather information from current chains.+    cl = fromMHG $ xs V.! i+    cr = fromMHG $ xs V.! j+    ll = link cl+    lr = link cr+    prl = prior ll+    prr = prior lr+    lhl = likelihood ll+    lhr = likelihood lr+    -- Swap the states.+    xl' = state lr+    xr' = state ll+    -- Compute new priors and likelihoods.+    prl' = priorFunction cl xl'+    prr' = priorFunction cr xr'+    lhl' = likelihoodFunction cl xl'+    lhr' = likelihoodFunction cr xr'+    -- Set the new links and the proposed state.+    ll' = Link xl' prl' lhl'+    lr' = Link xr' prr' lhr'+    cl' = cl {link = ll'}+    cr' = cr {link = lr'}+    xs' = xs V.// [(i, MHG cl'), (j, MHG cr')]+    -- Compute the Metropolis ratio.+    nominator = prl' * prr' * lhl' * lhr'+    denominator = prl * prr * lhl * lhr+    q = nominator / denominator++mc3ProposeSwap ::+  MC3 a ->+  -- Index of left chain.+  Int ->+  IO (MC3 a)+mc3ProposeSwap a i = do+  -- 1. Sample new state and get the Metropolis ratio.+  let (!y, !r) = swapWith i (i + 1) $ mc3MHGChains a+  -- 2. Accept or reject.+  accept <- mhgAccept r g+  if accept+    then do+      let !ac' = pushA i True (mc3SwapAcceptance a)+      return $ a {mc3MHGChains = y, mc3SwapAcceptance = ac'}+    else do+      let !ac' = pushA i False (mc3SwapAcceptance a)+      return $ a {mc3SwapAcceptance = ac'}+  where+    g = mc3Generator a++-- TODO: Splimix. 'mc3Iterate' is actually not parallel, but concurrent because+-- of the IO constraint. Use pure parallel code when we have a pure generator.+-- However, we have honor the mutable traces.+mc3Iterate ::+  ToJSON a =>+  ParallelizationMode ->+  MC3 a ->+  IO (MC3 a)+mc3Iterate pm a = do+  -- 1. Maybe propose swaps.+  --+  -- NOTE: Swaps have to be proposed first, because the traces are automatically+  -- updated at step 2.+  let s = mc3Settings a+  a' <-+    if mc3Iteration a `mod` fromSwapPeriod (mc3SwapPeriod s) == 0+      then do+        let n = V.length $ mc3MHGChains a+            is = [0 .. n - 2]+            ns = fromNSwaps $ mc3NSwaps s+        is' <- shuffle is $ mc3Generator a+        foldM mc3ProposeSwap a (take ns is')+      else return a+  -- 2. Iterate all chains and increment iteration.+  mhgs <- case pm of+    Sequential -> V.mapM (aIterate pm) (mc3MHGChains a')+    Parallel ->+      -- See 'Control.Monad.Parallel' of package 'monad-parallel'. Go via a+      -- list, and use 'forkIO'.+      V.fromList <$> P.mapM (aIterate pm) (V.toList (mc3MHGChains a'))+  let i = mc3Iteration a'+  return $ a' {mc3MHGChains = mhgs, mc3Iteration = succ i}++tuneBeta ::+  -- The old reciprocal temperatures are needed to retrieve the old ratios.+  ReciprocalTemperatures ->+  -- Index i of left chain. Change the reciprocal temperature of chain (i+1).+  Int ->+  -- Exponent xi of the reciprocal temperature ratio.+  Double ->+  -- The new reciprocal temperatures are updated incrementally using the+  -- reciprocal temperature ratios during the fold (see 'mc3AutoTune' below).+  ReciprocalTemperatures ->+  ReciprocalTemperatures+tuneBeta bsOld i xi bsNew = bsNew U.// [(j, brNew)]+  where+    j = i + 1+    blOld = bsOld U.! i+    brOld = bsOld U.! j+    blNew = bsNew U.! i+    -- The new ratio is in (0,1).+    rNew = (brOld / blOld) ** xi+    brNew = blNew * rNew++mc3AutoTune :: ToJSON a => MC3 a -> MC3 a+mc3AutoTune a = a {mc3MHGChains = mhgs'', mc3ReciprocalTemperatures = bs'}+  where+    mhgs = mc3MHGChains a+    -- 1. Auto tune all chains.+    mhgs' = V.map aAutoTune mhgs+    -- 2. Auto tune temperatures.+    optimalRate = getOptimalRate PDimensionUnknown+    currentRates = acceptanceRates $ mc3SwapAcceptance 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+    -- ratios, which is the same. Also, by working with ratios in (0,1) of+    -- neighboring chains, we ensure the monotonicity of the reciprocal+    -- temperatures.+    --+    -- The factor (1/2) was determined by a few tests and is otherwise+    -- absolutely arbitrary.+    xi i = exp $ (/ 2) $ (currentRates M.! i) - optimalRate+    bs = mc3ReciprocalTemperatures a+    n = fromNChains $ mc3NChains $ mc3Settings a+    -- Do not change the temperature, and the prior and likelihood functions of+    -- the cold chain.+    bs' = foldl' (\xs j -> tuneBeta bs j (xi j) xs) bs [0 .. n - 2]+    coldChain = fromMHG $ V.head mhgs'+    coldPrF = priorFunction coldChain+    coldLhF = likelihoodFunction coldChain+    mhgs'' =+      V.head mhgs'+        `V.cons` V.zipWith+          (setReciprocalTemperature coldPrF coldLhF)+          (V.convert $ U.tail bs')+          (V.tail mhgs')++mc3ResetAcceptance :: ToJSON a => MC3 a -> MC3 a+mc3ResetAcceptance a = a'+  where+    -- 1. Reset acceptance of all chains.+    mhgs' = V.map aResetAcceptance (mc3MHGChains a)+    -- 2. Reset acceptance of swaps.+    ac' = resetA $ mc3SwapAcceptance a+    --+    a' = a {mc3MHGChains = mhgs', mc3SwapAcceptance = ac'}++-- Information in cycle summary:+--+-- - The complete summary of the cycle of the cold chain.+--+-- - The combined acceptance rate of proposals within the hot chains.+--+-- - The temperatures of the chains and the acceptance rates of the state swaps.+mc3SummarizeCycle :: ToJSON a => MC3 a -> BL.ByteString+mc3SummarizeCycle a =+  BL.intercalate "\n" $+    [ "MC3: Cycle of cold chain.",+      coldMHGCycleSummary+    ]+      ++ [ "MC3: Average acceptance rate across all chains: " <> BL.fromStrict (BC.toFixed 2 ar)+           | not $ isNaN ar+         ]+      ++ [ "MC3: Reciprocal temperatures of the chains: " <> BL.intercalate ", " bsB <> ".",+           "MC3: Summary of state swaps.",+           "MC3: The swap period is " <> swapPeriodB <> ".",+           "MC3: The state swaps are executed in random order.",+           proposalHeader,+           proposalHLine+         ]+      ++ [ summarizeProposal+             (PName $ show i ++ " <-> " ++ show (i + 1))+             (PDescription "Swap states between chains")+             (PWeight 1)+             (Just $ bs U.! (i + 1))+             PDimensionUnknown+             (acceptanceRate i swapAcceptance)+           | i <- [0 .. n - 2]+         ]+      ++ [proposalHLine]+  where+    mhgs = mc3MHGChains a+    coldMHGCycleSummary = aSummarizeCycle $ V.head mhgs+    cs = V.map fromMHG mhgs+    as = V.map (acceptanceRates . acceptance) cs+    vAr = V.map (\m -> sum m / fromIntegral (length m)) as+    ar = V.sum vAr / fromIntegral (V.length vAr)+    bs = mc3ReciprocalTemperatures a+    bsB = map (BL.fromStrict . BC.toFixed 2) $ U.toList bs+    swapPeriod = fromSwapPeriod $ mc3SwapPeriod $ mc3Settings a+    swapPeriodB = BB.toLazyByteString $ BB.intDec swapPeriod+    swapAcceptance = mc3SwapAcceptance a+    n = fromNChains $ mc3NChains $ mc3Settings a++-- 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)+  return $ a {mc3MHGChains = mhgs'}++mc3ExecuteMonitors ::+  ToJSON a =>+  Verbosity ->+  -- Starting time.+  UTCTime ->+  -- Total number of iterations.+  Int ->+  MC3 a ->+  IO (Maybe BL.ByteString)+mc3ExecuteMonitors vb t0 iTotal a = V.head <$> V.imapM f (mc3MHGChains a)+  where+    -- The first chain honors verbosity.+    f 0 = aExecuteMonitors vb t0 iTotal+    -- All other chains are to be quiet.+    f _ = aExecuteMonitors Quiet t0 iTotal++mc3StdMonitorHeader :: ToJSON a => MC3 a -> BL.ByteString+mc3StdMonitorHeader = aStdMonitorHeader . V.head . mc3MHGChains++mc3CloseMonitors :: ToJSON a => MC3 a -> IO (MC3 a)+mc3CloseMonitors a = do+  mhgs' <- V.mapM aCloseMonitors $ mc3MHGChains a+  return $ a {mc3MHGChains = mhgs'}

@@ -0,0 +1,254 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  Mcmc.Algorithm.Metropolis+-- Description :  Metropolis-Hastings-Green algorithm+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue May  5 20:11:30 2020.+--+-- The Metropolis-Hastings-Green ('MHG') algorithm.+--+-- For example, see Geyer, C. J., Introduction to Markov chain Monte Carlo, In+-- Handbook of Markov Chain Monte Carlo (pp. 45) (2011). CRC press.+module Mcmc.Algorithm.Metropolis+  ( MHG (..),+    mhg,+    mhgSave,+    mhgLoad,+    mhgAccept,+  )+where++import Codec.Compression.GZip+import Control.Monad+import Control.Monad.IO.Class+import Data.Aeson+import qualified Data.ByteString.Lazy.Char8 as BL+import Data.Time+import Mcmc.Algorithm+import Mcmc.Chain.Chain+import Mcmc.Chain.Link+import Mcmc.Chain.Save+import Mcmc.Chain.Trace+import Mcmc.Monitor+import Mcmc.Proposal+import Mcmc.Settings+import Numeric.Log+import System.Random.MWC+import Text.Printf+import Prelude hiding (cycle)++-- | The MHG algorithm.+newtype MHG a = MHG {fromMHG :: Chain a}++instance ToJSON a => Algorithm (MHG a) where+  aName = const "Metropolis-Hastings-Green (MHG)"+  aIteration = iteration . fromMHG+  aIterate = mhgIterate+  aAutoTune = mhgAutoTune+  aResetAcceptance = mhgResetAcceptance+  aSummarizeCycle = mhgSummarizeCycle+  aOpenMonitors = mhgOpenMonitors+  aExecuteMonitors = mhgExecuteMonitors+  aStdMonitorHeader = mhgStdMonitorHeader+  aCloseMonitors = mhgCloseMonitors+  aSave = mhgSave++-- NOTE: IO is required because the trace is mutable.++-- | Initialize an MHG algorithm.+mhg ::+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  -- | The initial state in the state space @a@.+  a ->+  -- | A source of randomness. For reproducible runs, make sure to use+  -- generators with the same seed.+  GenIO ->+  IO (MHG a)+mhg pr lh cc mn i0 g = do+  -- The trace is a mutable vector and the mutable state needs to be handled by+  -- a monad.+  tr <- replicateT traceLength l0+  return $ MHG $ Chain 0 l0 0 tr ac g 0 pr lh cc mn+  where+    l0 = Link i0 (pr i0) (lh i0)+    ac = emptyA $ ccProposals cc+    batchMonitorSizes = map getMonitorBatchSize $ mBatches mn+    traceLength = maximum $ 1 : batchMonitorSizes++mhgFn :: AnalysisName -> FilePath+mhgFn (AnalysisName nm) = nm ++ ".mhg"++-- | Save an MHG algorithm.+mhgSave ::+  ToJSON a =>+  AnalysisName ->+  MHG a ->+  IO ()+mhgSave nm (MHG c) = do+  savedChain <- toSavedChain c+  BL.writeFile (mhgFn nm) $ compress $ encode savedChain++-- | Load an MHG algorithm.+--+-- See 'Mcmc.Mcmc.mcmcContinue'.+mhgLoad ::+  FromJSON a =>+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  AnalysisName ->+  IO (MHG a)+mhgLoad pr lh cc mn nm = do+  savedChain <- eitherDecode . decompress <$> BL.readFile (mhgFn nm)+  chain <- either error (fromSavedChain pr lh cc mn) savedChain+  return $ MHG chain++-- The MHG ratio.+--+-- 'Infinity' if fX is zero. In this case, the proposal is always accepted.+--+-- 'NaN' if (fY or q) and fX are zero. In this case, the proposal is always+-- rejected.++-- There is a discrepancy between authors saying that one should (a) always+-- accept the new state when the current posterior is zero (Chapter 4 of the+-- Handbook of Markov Chain Monte Carlo), or (b) almost surely reject the+-- proposal when either fY or q are zero (Chapter 1). Since I trust the author+-- of Chapter 1 (Charles Geyer) I choose to follow option (b).+mhgRatio :: Log Double -> Log Double -> Log Double -> Log Double -> Log Double+-- q = qYX / qXY * jXY; see 'ProposalSimple'.+-- j = Jacobian.+mhgRatio fX fY q j = fY / fX * q * j+{-# INLINE mhgRatio #-}++-- | Accept or reject a proposal with given MHG ratio?+mhgAccept :: Log Double -> GenIO -> IO Bool+mhgAccept r g+  | ln r >= 0.0 = return True+  | otherwise = do+    b <- uniform g+    return $ b < exp (ln r)++mhgPropose :: MHG a -> Proposal a -> IO (MHG a)+mhgPropose (MHG c) p = do+  -- 1. Sample new state.+  (!y, !q, !j) <- liftIO $ s x g+  -- 2. Calculate Metropolis-Hastings-Green ratio.+  let !pY = pF y+      !lY = lF y+      !r = mhgRatio (pX * lX) (pY * lY) q j+  -- 3. Accept or reject.+  -- if ln r >= 0.0+  --   then do+  --     let !ac' = pushA p True ac+  --     return $ MHG $ c {link = Link y pY lY, acceptance = ac'}+  --   else do+  --     b <- uniform g+  --     if b < exp (ln r)+  --       then do+  --         let !ac' = pushA p True ac+  --         return $ MHG $ c {link = Link y pY lY, acceptance = ac'}+  --       else do+  --         let !ac' = pushA p False ac+  --         return $ MHG $ c {acceptance = pushA p False ac'}+  accept <- mhgAccept r g+  if accept+    then do+      let !ac' = pushA p True ac+      return $ MHG $ c {link = Link y pY lY, acceptance = ac'}+    else do+      let !ac' = pushA p False ac+      return $ MHG $ c {acceptance = pushA p False ac'}+  where+    s = pSimple p+    (Link x pX lX) = link c+    pF = priorFunction c+    lF = likelihoodFunction c+    ac = acceptance c+    g = generator c++mhgPush :: MHG a -> IO (MHG a)+mhgPush (MHG c) = do+  t' <- pushT i t+  return $ MHG c {trace = t', iteration = succ n}+  where+    i = link c+    t = trace c+    n = iteration c++-- Ignore the number of capabilities. I have tried a lot of stuff, but the MHG+-- algorithm is just inherently sequential. Parallelization can be achieved by+-- having parallel prior and/or likelihood functions, or by using algorithms+-- running parallel chains such as 'MC3'.+mhgIterate :: ParallelizationMode -> MHG a -> IO (MHG a)+mhgIterate _ a = do+  ps <- orderProposals cc g+  a' <- foldM mhgPropose a ps+  mhgPush a'+  where+    c = fromMHG a+    cc = cycle c+    g = generator c++mhgAutoTune :: MHG a -> MHG a+mhgAutoTune (MHG c) = MHG $ c {cycle = autoTuneCycle ac cc}+  where+    ac = acceptance c+    cc = cycle c++mhgResetAcceptance :: MHG a -> MHG a+mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac}+  where+    ac = acceptance c++mhgSummarizeCycle :: MHG a -> BL.ByteString+mhgSummarizeCycle (MHG c) = summarizeCycle ac cc+  where+    cc = cycle c+    ac = acceptance c++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'}+  where+    m = monitor c+    pre = fromAnalysisName nm+    suf = printf "%02d" $ chainId c++mhgExecuteMonitors ::+  Verbosity ->+  -- Starting time.+  UTCTime ->+  -- Total number of iterations.+  Int ->+  MHG a ->+  IO (Maybe BL.ByteString)+mhgExecuteMonitors vb t0 iTotal (MHG c) = mExec vb i i0 t0 tr iTotal m+  where+    i = iteration c+    i0 = start c+    tr = trace c+    m = monitor c++mhgStdMonitorHeader :: MHG a -> BL.ByteString+mhgStdMonitorHeader (MHG c) = msHeader (mStdOut $ monitor c)++mhgCloseMonitors :: MHG a -> IO (MHG a)+mhgCloseMonitors (MHG c) = do+  m' <- mClose m+  return $ MHG $ c {monitor = m'}+  where+    m = monitor c

@@ -0,0 +1,95 @@+-- |+-- Module      :  Mcmc.Chain.Chain+-- Description :  Simple representation of a Markov chain+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue May  5 18:01:15 2020.+module Mcmc.Chain.Chain+  ( PriorFunction,+    noPrior,+    LikelihoodFunction,+    noLikelihood,+    Chain (..),+  )+where++-- Note: It is not necessary to add another type @b@ to store supplementary+-- information about the chain. The information can just be stored in @a@+-- equally well.++import Mcmc.Chain.Link+import Mcmc.Chain.Trace+import Mcmc.Monitor+import Mcmc.Proposal+import Numeric.Log+import System.Random.MWC hiding (save)+import Prelude hiding (cycle)++-- | Prior function.+type PriorFunction a = a -> Log Double++-- | Flat prior function. Useful for testing and debugging.+noPrior :: PriorFunction a+noPrior = const 1.0++-- | Likelihood function.+type LikelihoodFunction a = a -> Log Double++-- | Flat likelihood function. Useful for testing and debugging.+noLikelihood :: LikelihoodFunction a+noLikelihood = const 1.0++-- | The chain contains all information to run an MCMC sampler.+--+-- The state of a chain has type @a@. If necessary, the type @a@ can also be+-- used to store auxiliary information.+--+-- For example, the chain stores information about the current 'Link' and+-- 'iteration', the 'Trace', the 'Acceptance' rates, and the random number+-- generator.+--+-- Further, the chain includes auxiliary variables and functions such as the+-- prior and likelihood functions, or 'Proposal's to move around the state space+-- and to 'Monitor' an MCMC run.+--+-- The 'Mcmc.Environment.Environment' of the chain is not stored externally.+data Chain a = Chain+  { -- Variables; saved.+    -- | Chain index; useful if more chains are run.+    chainId :: 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.+    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.+    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),+    -- | The random number generator.+    generator :: GenIO,+    --+    -- Variables and functions; not saved.++    -- | Starting iteration of the chain; used to calculate run time and ETA.+    start :: Int,+    -- | The prior function. The un-normalized posterior is the product of the+    -- prior and the likelihood.+    priorFunction :: PriorFunction a,+    -- | The likelihood function. The un-normalized posterior is the product of+    -- the prior and the likelihood.+    likelihoodFunction :: LikelihoodFunction a,+    -- | A set of 'Proposal's form a 'Cycle'.+    cycle :: Cycle a,+    -- | A 'Monitor' observing the chain.+    monitor :: Monitor a+  }

@@ -0,0 +1,47 @@+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  Mcmc.Chain.Link+-- Description :  The state combined with auxiliary variables+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed May 20 09:10:27 2020.+module Mcmc.Chain.Link+  ( Link (..),+  )+where++import Data.Aeson+import Data.Aeson.Types+import Numeric.Log++-- | Link of a Markov chain. For reasons of computational efficiency, each state+-- is associated with the corresponding prior and likelihood.+data Link a = Link+  { -- | The current state in the state space @a@.+    state :: a,+    -- | The current prior.+    prior :: Log Double,+    -- | The current likelihood.+    likelihood :: Log Double+  }+  deriving (Eq, Ord, Show, Read)++instance ToJSON a => ToJSON (Link a) where+  toJSON (Link x (Exp p) (Exp l)) = object ["s" .= x, "p" .= p, "l" .= l]+  toEncoding (Link x (Exp p) (Exp l)) = pairs ("s" .= x <> "p" .= p <> "l" .= l)++link :: FromJSON a => Object -> Parser (Link a)+link v = do+  s <- v .: "s"+  p <- v .: "p"+  l <- v .: "l"+  return $ Link s (Exp p) (Exp l)++instance FromJSON a => FromJSON (Link a) where+  parseJSON = withObject "Link" link

@@ -0,0 +1,108 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TemplateHaskell #-}++-- |+-- Module      :  Mcmc.Chain.Save+-- Description :  Save and load a Markov chain+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--++-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue Jun 16 10:18:54 2020.+--+-- 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+  ( SavedChain (..),+    toSavedChain,+    fromSavedChain,+  )+where++import Control.Monad+import Data.Aeson+import Data.Aeson.TH+import Data.List hiding (cycle)+import qualified Data.Map as M+import Data.Maybe+import qualified Data.Vector as VB+import qualified Data.Vector.Unboxed as VU+import Data.Word+import Mcmc.Chain.Chain+import Mcmc.Chain.Link+import Mcmc.Chain.Trace+import Mcmc.Internal.Random+import Mcmc.Monitor+import Mcmc.Proposal+import Prelude hiding (cycle)+import qualified Data.Stack.Circular as C++-- | Storable values of a Markov chain.+--+-- See 'toSavedChain'.+data SavedChain a = SavedChain+  {+    savedId :: Int,+    savedLink :: Link a,+    savedIteration :: Int,+    savedTrace :: C.Stack VB.Vector (Link a),+    savedAcceptance :: Acceptance Int,+    savedSeed :: VU.Vector Word32,+    savedTuningParameters :: [Maybe Double]+  }+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''SavedChain)++-- | Save a chain.+toSavedChain ::+  Chain a ->+  IO (SavedChain a)+toSavedChain (Chain ci it i tr ac g _ _ _ cc _) = do+  g' <- saveGen g+  tr' <- freezeT tr+  return $ SavedChain ci it i tr' ac' g' ts+  where+    ps = ccProposals cc+    ac' = transformKeysA ps [0 ..] ac+    ts = [fmap tParam mt | mt <- map pTuner ps]++-- | Load a saved chain.+--+-- Recompute and check the prior and likelihood for the last state because the+-- functions may have changed. Of course, we cannot test for the same function,+-- but having the same prior and likelihood at the last state is already a good+-- indicator.+fromSavedChain ::+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  SavedChain a ->+  IO (Chain a)+fromSavedChain pr lh cc mn (SavedChain ci it i tr ac' g' ts)+  | pr (state it) /= prior it =+    error "fromSave: Provided prior function does not match the saved prior."+  | lh (state it) /= likelihood it =+    error "fromSave: Provided likelihood function does not match the saved likelihood."+  | otherwise = do+      g <- loadGen g'+      tr' <- thawT tr+      return $ Chain ci it i tr' ac g i pr lh cc' mn+  where+    ac = transformKeysA [0 ..] (ccProposals cc) ac'+    getTuningF mt = case mt of+      Nothing -> const 1.0+      Just t -> const t+    cc' =+      tuneCycle+        ( M.map getTuningF $+            M.fromList $+              zip (ccProposals cc) ts+        )+        cc

@@ -0,0 +1,75 @@+-- |+-- Module      :  Mcmc.Chain.Trace+-- Description :  History of a Markov chain+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed May 20 09:11:25 2020.+module Mcmc.Chain.Trace+  ( Trace,+    replicateT,+    lengthT,+    pushT,+    headT,+    takeT,+    freezeT,+    thawT,+  )+where++import Control.Monad.Primitive+import qualified Data.Stack.Circular as C+import qualified Data.Vector as VB+import Mcmc.Chain.Link++-- | A 'Trace' is a mutable circular stack that passes through a list of states+-- with associated priors and likelihoods called 'Link's.+newtype Trace a = Trace {fromTrace :: C.MStack VB.Vector RealWorld (Link a)}++-- | Initialize a trace of given length by replicating the same value.+--+-- Be careful not to compute summary statistics before pushing enough values.+--+-- Call 'error' if the maximum size is zero or negative.+replicateT :: Int -> Link a -> IO (Trace a)+replicateT n l = Trace <$> C.replicate n l++-- | 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+  s' <- C.push x (fromTrace t)+  return $ Trace s'+{-# INLINEABLE pushT #-}++-- | 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'.+takeT :: Int -> Trace a -> IO (VB.Vector (Link a))+takeT k = C.take k . fromTrace++-- | 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'.+thawT :: C.Stack VB.Vector (Link a) -> IO (Trace a)+thawT t = Trace <$> C.thaw t

@@ -0,0 +1,49 @@+-- |+-- Module      :  Mcmc.Environment+-- Description :  Runtime environment+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue Nov 17 11:00:09 2020.+module Mcmc.Environment+  ( Environment (..),+    initializeEnvironment,+  )+where++import Data.Time.Clock+import Mcmc.Settings+import System.IO++-- | The environment of an MCMC run.+data Environment = Environment+  { settings :: Settings,+    -- | We have to use 'Maybe' here, because we do not want to open any log+    -- file when being 'Quiet'.+    logHandle :: Maybe Handle,+    -- | Used to calculate the ETA.+    startingTime :: UTCTime+  }+  deriving (Eq, Show)++-- | Initialize the environment.+--+-- Open log file, get current time.+initializeEnvironment ::+  Settings ->+  IO Environment+initializeEnvironment s = do+  t <- getCurrentTime+  mh <- case sVerbosity s of+    Quiet -> return Nothing+    _ -> do+      h <- openWithExecutionMode em fn+      return $ Just h+  return $ Environment s mh t+  where+    fn = fromAnalysisName (sAnalysisName s) ++ ".log"+    em = sExecutionMode s

@@ -10,7 +10,8 @@ -- -- Creation date: Mon Aug  3 10:46:27 2020. module Mcmc.Internal.ByteString-  ( alignRightWith,+  ( alignRightWithNoTrim,+    alignRightWith,     alignRight,     alignLeftWith,     alignLeft,@@ -19,11 +20,16 @@  import qualified Data.ByteString.Lazy.Char8 as BL +-- | For a given width, align string to the right; use given fill character.+alignRightWithNoTrim :: Char -> Int -> BL.ByteString -> BL.ByteString+alignRightWithNoTrim c n s = BL.replicate (fromIntegral n - l) c <> s+  where+    l = BL.length s+ -- | For a given width, align string to the right; use given fill character; -- trim on the left if string is longer. alignRightWith :: Char -> Int -> BL.ByteString -> BL.ByteString-alignRightWith c n s =-  BL.replicate (fromIntegral n - l) c <> BL.take (fromIntegral n) s+alignRightWith c n s = BL.replicate (fromIntegral n - l) c <> BL.take (fromIntegral n) s   where     l = BL.length s 

@@ -0,0 +1,48 @@+{-# LANGUAGE ScopedTypeVariables #-}++-- |+-- Module      :  Mcmc.Internal.Random+-- Description :  Tools for random calculations+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed Nov 25 07:14:52 2020.+module Mcmc.Internal.Random+  ( splitGen,+    saveGen,+    loadGen,+  )+where++import Control.Monad+import Control.Monad.Primitive+import qualified Data.Vector.Unboxed as V+import Data.Word+import System.Random.MWC++-- | Split a generator.+--+-- Splitting an MWC pseudo number generator is not good practice. However, I+-- have to go with this solution for now, and wait for proper support of+-- spittable pseudo random number generators such as @splitmix@.+splitGen :: PrimMonad m => Int -> Gen (PrimState m) -> m [Gen (PrimState m)]+splitGen n gen+  | n <= 0 = return []+  | otherwise = do+    seeds :: [V.Vector Word32] <- replicateM n $ uniformVector gen 256+    mapM initialize seeds++-- TODO: Splitmix. Remove or amend these functions as soon as split mix is used+-- and is available with the statistics package.++-- | Save a generator to a seed.+saveGen :: GenIO -> IO (V.Vector Word32)+saveGen = fmap fromSeed . save++-- | Load a generator from a seed.+loadGen :: V.Vector Word32 -> IO GenIO+loadGen = restore . toSeed

@@ -13,51 +13,37 @@ -- From https://wiki.haskell.org/Random_shuffle. module Mcmc.Internal.Shuffle   ( shuffle,-    shuffleN,     grabble,   ) where  import Control.Monad import Control.Monad.ST-import Data.Vector (Vector) import qualified Data.Vector as V import qualified Data.Vector.Mutable as M import System.Random.MWC-  ( GenIO,-    uniformR,-  ) --- | Shuffle a list.-shuffle :: [a] -> GenIO -> IO [a]-shuffle xs g = head <$> grabble xs 1 (length xs) g---- | Shuffle a list @n@ times.-shuffleN :: [a] -> Int -> GenIO -> IO [[a]]-shuffleN xs n = grabble xs n (length xs)+-- 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. --- -- Using System.Random.Shuffle. Speed is the same, so stay without additional dependency.--- -- | Shuffle a list @n@ times.--- shuffleN :: [a] -> Int -> GenIO -> IO [[a]]--- shuffleN xs n g = replicateM n $ fmap (shuffle xs) (rseqM (length xs - 1) g)---   where---     rseqM :: Int -> GenIO -> IO [Int]---     rseqM 0 _ = return []---     rseqM i gen = liftM2 (:) (uniformR (0, i) gen) (rseqM (i - 1) gen)+-- | Shuffle a vector.+shuffle :: [a] -> GenIO -> IO [a]+shuffle xs = grabble xs (length xs)  -- | @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@ -- elements from @xs@, without replacement, and that @m@ times.-grabble :: [a] -> Int -> Int -> GenIO -> IO [[a]]-grabble xs m n gen = do-  swapss <- replicateM m $-    forM [0 .. min (l - 1) n] $ \i -> do-      j <- uniformR (i, l) gen-      return (i, j)-  return $ map (V.toList . V.take n . swapElems (V.fromList xs)) swapss+grabble :: [a] -> Int -> GenIO -> IO [a]+grabble xs m gen = do+  swaps <- forM [0 .. min (l - 1) m] $ \i -> do+    j <- uniformR (i, l) gen+    return (i, j)+  return $ (V.toList . V.take m . swapElems (V.fromList xs)) swaps   where     l = length xs - 1 -swapElems :: Vector a -> [(Int, Int)] -> Vector a+swapElems :: V.Vector a -> [(Int, Int)] -> V.Vector a swapElems xs swaps = runST $ do   mxs <- V.unsafeThaw xs   mapM_ (uncurry $ M.unsafeSwap mxs) swaps

@@ -1,48 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}---- |--- Module      :  Mcmc.Item--- Description :  Links of Markov chains--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Wed May 20 09:10:27 2020.-module Mcmc.Item-  ( Item (..),-  )-where--import Data.Aeson-import Data.Aeson.Types-import Numeric.Log---- | An 'Item', or link of the Markov chain. For reasons of computational--- efficiency, each state is associated with the corresponding prior and--- likelihood.-data Item a = Item-  { -- | The current state in the state space @a@.-    state :: a,-    -- | The current prior.-    prior :: Log Double,-    -- | The current likelihood.-    likelihood :: Log Double-  }-  deriving (Eq, Ord, Show, Read)--instance ToJSON a => ToJSON (Item a) where-  toJSON (Item x (Exp p) (Exp l)) = object ["s" .= x, "p" .= p, "l" .= l]-  toEncoding (Item x (Exp p) (Exp l)) = pairs ("s" .= x <> "p" .= p <> "l" .= l)--item :: FromJSON a => Object -> Parser (Item a)-item v = do-  s <- v .: "s"-  p <- v .: "p"-  l <- v .: "l"-  return $ Item s (Exp p) (Exp l)--instance FromJSON a => FromJSON (Item a) where-  parseJSON = withObject "Item" item

@@ -2,7 +2,7 @@  -- | -- Module      :  Mcmc.Mcmc--- Description :  Mcmc helpers+-- Description :  Framework for running Markov chain Monte Carlo samplers -- Copyright   :  (c) Dominik Schrempf, 2020 -- License     :  GPL-3.0-or-later --@@ -12,265 +12,277 @@ -- -- Creation date: Fri May 29 10:19:45 2020. ----- Functions to work with the 'Mcmc' state transformer.+-- This module provides the general framework for running MCMC samplers. By+-- design choice this module is agnostic about the details of the used+-- 'Algorithm'. module Mcmc.Mcmc-  ( Mcmc,-    mcmcOutB,-    mcmcOutS,-    mcmcWarnB,-    mcmcWarnS,-    mcmcInfoB,-    mcmcInfoS,-    mcmcDebugB,-    mcmcDebugS,-    mcmcAutotune,-    mcmcClean,-    mcmcResetA,-    mcmcSummarizeCycle,-    mcmcReport,-    mcmcMonitorExec,-    mcmcRun,+  ( mcmc,+    mcmcContinue,   ) where  import Control.Monad import Control.Monad.IO.Class-import Control.Monad.Trans.State-import Data.Aeson+-- import Control.Monad.Trans.RWS.CPS+import Control.Monad.Trans.Reader import qualified Data.ByteString.Lazy.Char8 as BL import Data.Maybe import Data.Time.Clock-import Data.Time.Format-import Mcmc.Item-import Mcmc.Monitor+import Mcmc.Algorithm+import Mcmc.Environment import Mcmc.Monitor.Time-import Mcmc.Proposal-import Mcmc.Save-import Mcmc.Status hiding (debug)-import Mcmc.Verbosity-import Numeric.Log-import System.Directory+import Mcmc.Settings import System.IO+import Text.Show.Pretty import Prelude hiding (cycle) --- | An Mcmc state transformer; usually fiddling around with this type is not--- required, but it is used by the different inference algorithms.-type Mcmc a = StateT (Status a) IO--msgPrepare :: Char -> BL.ByteString -> BL.ByteString-msgPrepare c t = BL.cons c $ ": " <> t+-- The MCMC algorithm has read access to an environment and uses an algorithm+-- transforming the state @a@.+type MCMC a = ReaderT Environment IO a --- | Write to standard output and log file.-mcmcOutB :: BL.ByteString -> Mcmc a ()-mcmcOutB msg = do-  h <- fromMaybe (error "mcmcOut: Log handle is missing.") <$> gets logHandle-  liftIO $ BL.putStrLn msg >> BL.hPutStrLn h msg+msgPrepare :: BL.ByteString -> BL.ByteString -> BL.ByteString+msgPrepare pref msg = BL.intercalate "\n" $ map (BL.append pref) $ BL.lines msg --- | Write to standard output and log file.-mcmcOutS :: String -> Mcmc a ()-mcmcOutS = mcmcOutB . BL.pack+-- Write to standard output and log file.+mcmcOutB :: BL.ByteString -> BL.ByteString -> MCMC ()+mcmcOutB pref msg = do+  h <- fromMaybe (error "mcmcOut: Log handle is missing.") <$> reader logHandle+  liftIO $ BL.putStrLn msg' >> BL.hPutStrLn h msg'+  where+    msg' = msgPrepare pref msg --- Perform warning action.-mcmcWarnA :: Mcmc a () -> Mcmc a ()-mcmcWarnA a = gets verbosity >>= \v -> info v a+-- -- Perform warning action.+-- mcmcWarnA :: MCMC a () -> MCMC a ()+-- mcmcWarnA a = reader (verbosity . settings) >>= \v -> when (v >= Warn) a --- | Print warning message.-mcmcWarnB :: BL.ByteString -> Mcmc a ()-mcmcWarnB = mcmcWarnA . mcmcOutB . msgPrepare 'W'+-- -- Print warning message.+-- mcmcWarnB :: BL.ByteString -> MCMC a ()+-- mcmcWarnB = mcmcWarnA . mcmcOutB . msgPrepare 'W' --- | Print warning message.-mcmcWarnS :: String -> Mcmc a ()-mcmcWarnS = mcmcWarnB . BL.pack+-- -- Print warning message.+-- mcmcWarnS :: String -> MCMC a ()+-- mcmcWarnS = mcmcWarnB . BL.pack  -- Perform info action.-mcmcInfoA :: Mcmc a () -> Mcmc a ()-mcmcInfoA a = gets verbosity >>= \v -> info v a+mcmcInfoA :: MCMC () -> MCMC ()+mcmcInfoA a = reader (sVerbosity . settings) >>= \v -> when (v >= Info) a --- | Print info message.-mcmcInfoB :: BL.ByteString -> Mcmc a ()-mcmcInfoB = mcmcInfoA . mcmcOutB . msgPrepare 'I'+-- Print info message.+mcmcInfoB :: BL.ByteString -> MCMC ()+mcmcInfoB = mcmcInfoA . mcmcOutB "I: " --- | Print info message.-mcmcInfoS :: String -> Mcmc a ()+-- Print info message.+mcmcInfoS :: String -> MCMC () mcmcInfoS = mcmcInfoB . BL.pack  -- Perform debug action.-mcmcDebugA :: Mcmc a () -> Mcmc a ()-mcmcDebugA a = gets verbosity >>= \v -> debug v a+mcmcDebugA :: MCMC () -> MCMC ()+mcmcDebugA a = reader (sVerbosity . settings) >>= \v -> when (v == Debug) a --- | Print debug message.-mcmcDebugB :: BL.ByteString -> Mcmc a ()-mcmcDebugB = mcmcDebugA . mcmcOutB . msgPrepare 'D'+-- Print debug message.+mcmcDebugB :: BL.ByteString -> MCMC ()+mcmcDebugB = mcmcDebugA . mcmcOutB "D: " --- | Print debug message.-mcmcDebugS :: String -> Mcmc a ()+-- Print debug message.+mcmcDebugS :: String -> MCMC () mcmcDebugS = mcmcDebugB . BL.pack --- | Auto tune the 'Proposal's in the 'Cycle' of the chain. Reset acceptance counts.--- See 'autotuneCycle'.-mcmcAutotune :: Mcmc a ()-mcmcAutotune = do-  mcmcDebugB "Auto tune."-  s <- get-  let a = acceptance s-      c = cycle s-      c' = autotuneCycle a c-  put $ s {cycle = c'}+mcmcReportTime :: MCMC ()+mcmcReportTime = do+  mcmcDebugB "Report time."+  ti <- reader startingTime+  mcmcInfoS $ "Starting time of MCMC sampler: " <> renderTime ti --- | Clean the state.-mcmcClean :: Mcmc a ()-mcmcClean = do-  s <- get-  let cl = cleaner s-      i = iteration s-  case cl of-    Just (Cleaner n f) | i `mod` n == 0 -> do-      mcmcDebugB "Clean state."-      let (Item st pr lh) = item s-      mcmcDebugS $-        "Old log prior and log likelihood: " ++ show (ln pr) ++ ", " ++ show (ln lh) ++ "."-      let prF = priorF s-          lhF = likelihoodF s-          st' = f st-          pr' = prF st'-          lh' = lhF st'-      mcmcDebugS $-        "New log prior and log likelihood: " ++ show (ln pr') ++ ", " ++ show (ln lh') ++ "."-      let dLogPr = abs $ ln pr - ln pr'-          dLogLh = abs $ ln lh - ln lh'-      when-        (dLogPr > 0.01)-        (mcmcWarnS $ "Log of old and new prior differ by " ++ show dLogPr ++ ".")-      when-        (dLogPr > 0.01)-        (mcmcWarnS $ "Log of old and new likelihood differ by " ++ show dLogLh ++ ".")-      put $ s {item = Item st' pr' lh'}-    _ -> return ()+mcmcExecute :: Algorithm a => a -> MCMC a+mcmcExecute a = do+  mcmcDebugB "Executing MCMC run."+  s <- reader settings+  a' <- case sExecutionMode s of+    Fail -> mcmcNewRun a+    Overwrite -> mcmcNewRun a+    Continue -> mcmcContinueRun a+  mcmcDebugB "Executed MCMC run."+  return a' --- | Reset acceptance counts.-mcmcResetA :: Mcmc a ()-mcmcResetA = do-  mcmcDebugB "Reset acceptance ratios."-  s <- get-  let a = acceptance s-  put $ s {acceptance = resetA a}+-- Reset acceptance counts.+mcmcResetAcceptance :: Algorithm a => a -> MCMC a+mcmcResetAcceptance a = do+  mcmcDebugB "Reset acceptance rates."+  return $ aResetAcceptance a --- | Print short summary of 'Proposal's in 'Cycle'. See 'summarizeCycle'.-mcmcSummarizeCycle :: Mcmc a BL.ByteString-mcmcSummarizeCycle = do-  a <- gets acceptance-  c <- gets cycle-  return $ summarizeCycle a c+-- Execute the monitors of the chain.+mcmcExecuteMonitors :: Algorithm a => a -> MCMC ()+mcmcExecuteMonitors a = do+  e <- ask+  let s = settings e+      vb = sVerbosity s+      t0 = startingTime e+      iTotal = burnInIterations (sBurnIn s) + fromIterations (sIterations s)+  mStdLog <- liftIO (aExecuteMonitors vb t0 iTotal a)+  forM_ mStdLog (mcmcOutB "   ") -fTime :: FormatTime t => t -> String-fTime = formatTime defaultTimeLocale "%B %-e, %Y, at %H:%M %P, %Z."+mcmcIterate :: Algorithm a => Int -> a -> MCMC a+mcmcIterate n a+  | n < 0 = error "mcmcIterate: Number of iterations is negative."+  | n == 0 = return a+  | otherwise = do+    p <- sParallelizationMode . settings <$> ask+    a' <- liftIO $ aIterate p a+    mcmcExecuteMonitors a'+    mcmcIterate (n -1) a' --- Open log file.-mcmcOpenLog :: Mcmc a ()-mcmcOpenLog = do-  s <- get-  let lfn = name s ++ ".log"-      n = iteration s-      frc = forceOverwrite s-  fe <- liftIO $ doesFileExist lfn-  mh <- liftIO $ case verbosity s of-    Quiet -> return Nothing-    _ ->-      Just <$> case (fe, n, frc) of-        (False, _, _) -> openFile lfn WriteMode-        (True, 0, True) -> openFile lfn WriteMode-        (True, 0, False) -> error "mcmcInit: Log file exists; use 'force' to overwrite output files."-        (True, _, _) -> openFile lfn AppendMode-  put s {logHandle = mh}-  mcmcDebugS $ "Log file name: " ++ lfn ++ "."-  mcmcDebugB "Log file opened."+mcmcNewRun :: Algorithm a => a -> MCMC a+mcmcNewRun a = do+  s <- reader settings+  mcmcInfoB "Start new MCMC sampler."+  mcmcInfoB "Initial state."+  mcmcInfoB $ aStdMonitorHeader a+  mcmcExecuteMonitors a+  mcmcInfoB $ aSummarizeCycle a+  a' <- mcmcBurnIn a+  a'' <- mcmcResetAcceptance a'+  let i = fromIterations $ sIterations s+  mcmcInfoS $ "Run chain for " ++ show i ++ " iterations."+  mcmcInfoB $ aStdMonitorHeader a''+  mcmcIterate i a'' --- Set the total number of iterations, the current time and open the 'Monitor's--- of the chain. See 'mOpen'.-mcmcInit :: Mcmc a ()-mcmcInit = do-  mcmcOpenLog-  s <- get-  -- Start time.-  t <- liftIO getCurrentTime-  mcmcInfoS $ "Start time: " <> fTime t-  -- Monitor.-  let m = monitor s-      n = iteration s-      nm = name s-      frc = forceOverwrite s-  m' <- if n == 0 then liftIO $ mOpen nm frc m else liftIO $ mAppend nm m-  put $ s {monitor = m', start = Just (n, t)}+mcmcContinueRun :: Algorithm a => a -> MCMC a+mcmcContinueRun a = do+  s <- reader settings+  let iTotal = fromIterations (sIterations s) + burnInIterations (sBurnIn s)+  mcmcInfoB "Continuation of MCMC sampler."+  let iCurrent = aIteration a+  mcmcInfoS $ "Current iteration: " ++ show iCurrent ++ "."+  mcmcInfoS $ "Total iterations: " ++ show iTotal ++ "."+  let di = iTotal - iCurrent+  mcmcInfoB $ aSummarizeCycle a+  mcmcInfoS $ "Run chain for " ++ show di ++ " iterations."+  mcmcInfoB $ aStdMonitorHeader a+  mcmcIterate di a --- | Report what is going to be done.-mcmcReport :: ToJSON a => Mcmc a ()-mcmcReport = do-  s <- get-  let b = burnInIterations s-      t = autoTuningPeriod s-      n = iterations s-      c = cleaner s-  case b of-    Just b' -> mcmcInfoS $ "Burn in for " <> show b' <> " iterations."-    Nothing -> return ()-  case t of-    Just t' -> mcmcInfoS $ "Auto tune every " <> show t' <> " iterations (during burn in only)."-    Nothing -> return ()-  case c of-    Just (Cleaner c' _) -> mcmcInfoS $ "Clean state every " <> show c' <> " iterations."-    Nothing -> return ()-  mcmcInfoS $ "Run chain for " <> show n <> " iterations."-  mcmcInfoB "Initial state."-  mcmcMonitorExec+mcmcBurnIn :: Algorithm a => a -> MCMC a+mcmcBurnIn a = do+  s <- reader settings+  case sBurnIn s of+    NoBurnIn -> do+      mcmcInfoS "No burn in."+      return a+    BurnInWithoutAutoTuning n -> do+      mcmcInfoS $ "Burn in for " <> show n <> " iterations."+      mcmcInfoS "Auto tuning is disabled."+      mcmcInfoB $ aStdMonitorHeader a+      a' <- mcmcIterate n a+      mcmcInfoB $ aSummarizeCycle a'+      mcmcInfoB "Burn in finished."+      return a'+    BurnInWithAutoTuning n t -> do+      mcmcInfoS $ "Burn in for " ++ show n ++ " iterations."+      mcmcInfoS $ "Auto tuning is enabled with a period of " ++ show t ++ "."+      mcmcInfoB $ aStdMonitorHeader a+      a' <- mcmcBurnInWithAutoTuning n t a+      mcmcInfoB "Burn in finished."+      return a' --- Save the status of an MCMC run. See 'saveStatus'.-mcmcSave :: ToJSON a => Mcmc a ()-mcmcSave = do-  s <- get-  case save s of-    Just n -> do-      mcmcInfoB $ "Save Markov chain with trace of length " <> BL.pack (show n) <> "."+-- Auto tune the proposals.+mcmcAutotune :: Algorithm a => a -> MCMC a+mcmcAutotune a = do+  mcmcDebugB "Auto tune."+  return $ aAutoTune a++mcmcBurnInWithAutoTuning :: Algorithm a => Int -> Int -> a -> MCMC a+mcmcBurnInWithAutoTuning b t a+  | b > t = do+    a' <- mcmcResetAcceptance a+    a'' <- mcmcIterate t a'+    mcmcDebugB $ aSummarizeCycle a''+    a''' <- mcmcAutotune a''+    mcmcDebugB $ aStdMonitorHeader a''+    mcmcBurnInWithAutoTuning (b - t) t a'''+  | otherwise = do+    a' <- mcmcResetAcceptance a+    a'' <- mcmcIterate b a'+    mcmcInfoB $ aSummarizeCycle a''+    mcmcInfoS $ "Acceptance rates calculated over the last " <> show b <> " iterations."+    return a''++mcmcInitialize :: Algorithm a => a -> MCMC a+mcmcInitialize a = do+  mcmcInfoS $ aName a ++ " algorithm."+  s <- settings <$> ask+  mcmcDebugB "Opening monitors."+  a' <- liftIO $ aOpenMonitors (sAnalysisName s) (sExecutionMode s) a+  mcmcDebugB "Monitors opened."+  return a'++-- Save the MCMC run.+mcmcSave :: Algorithm a => a -> MCMC ()+mcmcSave a = do+  s <- reader settings+  case sSaveMode s of+    NoSave -> mcmcInfoB "Do not save the MCMC analysis."+    Save -> do+      mcmcInfoB "Save settings."+      liftIO $ settingsSave s+      let nm = sAnalysisName s+      mcmcInfoB "Save compressed MCMC analysis."       mcmcInfoB "For long traces, or complex objects, this may take a while."-      liftIO $ saveStatus (name s <> ".mcmc") s-      mcmcInfoB "Done saving Markov chain."-    Nothing -> mcmcInfoB "Do not save the Markov chain."+      liftIO $ aSave nm a+      mcmcInfoB "Markov chain saved." --- | Execute the 'Monitor's of the chain. See 'mExec'.-mcmcMonitorExec :: ToJSON a => Mcmc a ()-mcmcMonitorExec = do-  s <- get-  let i = iteration s-      j = iterations s + fromMaybe 0 (burnInIterations s)-      m = monitor s-      (ss, st) = fromMaybe (error "mcmcMonitorExec: Starting state and time not set.") (start s)-      tr = trace s-      vb = verbosity s-  mt <- liftIO $ mExec vb i ss st tr j m-  forM_ mt mcmcOutB+-- Report and finish up.+mcmcClose :: Algorithm a => a -> MCMC a+mcmcClose a = do+  mcmcDebugB "Closing MCMC run."+  mcmcInfoB $ aSummarizeCycle a+  mcmcInfoS $ aName a ++ " algorithm finished."+  mcmcSave a+  ti <- reader startingTime+  te <- liftIO getCurrentTime+  let dt = te `diffUTCTime` ti+  mcmcInfoB $ "Wall clock run time: " <> renderDuration dt <> "."+  mcmcInfoS $ "End time: " <> renderTime te+  a' <- liftIO $ aCloseMonitors a+  h <- reader logHandle+  liftIO $ forM_ h hClose+  return a' --- Close the 'Monitor's of the chain. See 'mClose'.-mcmcClose :: ToJSON a => Mcmc a ()-mcmcClose = do-  s <- get-  mcmcSummarizeCycle >>= mcmcInfoB-  mcmcInfoB "Metropolis-Hastings sampler finished."-  let m = monitor s-  m' <- liftIO $ mClose m-  put $ s {monitor = m'}-  mcmcSave-  t <- liftIO getCurrentTime-  let rt = case start s of-        Nothing -> error "mcmcClose: Start time not set."-        Just (_, st) -> t `diffUTCTime` st-  mcmcInfoB $ "Wall clock run time: " <> renderDuration rt <> "."-  mcmcInfoS $ "End time: " <> fTime t-  case logHandle s of-    Just h -> liftIO $ hClose h-    Nothing -> return ()+-- Initialize the run, execute the run, and close the run.+mcmcRun :: Algorithm a => a -> MCMC a+mcmcRun a = do+  mcmcDebugB "The settings are:"+  reader settings >>= mcmcDebugS . ppShow --- | Run an MCMC algorithm.-mcmcRun :: ToJSON a => Mcmc a () -> Status a -> IO (Status a)-mcmcRun algorithm = execStateT $ do-  mcmcInit-  algorithm-  mcmcClose+  -- Initialize.+  a' <- mcmcInitialize a+  mcmcReportTime++  -- Execute.+  a'' <- mcmcExecute a'++  -- Close.+  mcmcClose a''++-- | Run an MCMC algorithm with given settings.+mcmc :: Algorithm a => Settings -> a -> IO a+mcmc s a = do+  settingsCheck s $ aIteration a+  e <- initializeEnvironment s+  runReaderT (mcmcRun a) e++-- | Continue an MCMC algorithm for the given number of iterations.+--+-- Currently, it is only possible to continue MCMC algorithms that have+-- completed successfully. This restriction is necessary, because for parallel+-- chains, it is hardly possible to ensure all chains are synchronized when the+-- process is killed.+--+-- See:+--+-- - 'Mcmc.Algorithm.Metropolis.mhgLoad'+--+-- - 'Mcmc.Algorithm.MC3.mc3Load'+mcmcContinue :: Algorithm a => Int -> Settings -> a -> IO a+mcmcContinue dn s = mcmc s'+  where+    n' = Iterations $ fromIterations (sIterations s) + dn+    s' = s {sIterations = n', sExecutionMode = Continue}

@@ -1,196 +0,0 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE OverloadedStrings #-}---- |--- Module      :  Mcmc.Metropolis--- Description :  Metropolis-Hastings at its best--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Tue May  5 20:11:30 2020.------ Metropolis-Hastings algorithm.-module Mcmc.Metropolis-  ( mh,-    mhContinue,-  )-where--import Control.Monad-import Control.Monad.IO.Class-import Control.Monad.Trans.State-import Data.Aeson-import Data.Maybe-import Mcmc.Item-import Mcmc.Mcmc-import Mcmc.Proposal-import Mcmc.Status-import Mcmc.Trace-import Numeric.Log-import System.Random.MWC-import Prelude hiding (cycle)---- The Metropolis-Hastings ratio.------ 'Infinity' if fX is zero. In this case, the proposal is always accepted.------ 'NaN' if (fY or q) and fX are zero. In this case, the proposal is always--- rejected.---- There is a discrepancy between authors saying that one should (a) always--- accept the new state when the current posterior is zero (Chapter 4 of the--- Handbook of Markov Chain Monte Carlo), or (b) almost surely reject the--- proposal when either fY or q are zero (Chapter 1). Since I trust the author--- of Chapter 1 (Charles Geyer) I choose to follow option (b).-mhRatio :: Log Double -> Log Double -> Log Double -> Log Double -> Log Double--- q = qYX / qXY * jXY; see 'ProposalSimple'.--- j = Jacobian.-mhRatio fX fY q j = fY / fX * q * j-{-# INLINE mhRatio #-}--mhPropose :: Proposal a -> Mcmc a ()-mhPropose m = do-  let p = pSimple m-  s <- get-  let (Item x pX lX) = item s-      pF = priorF s-      lF = likelihoodF s-      a = acceptance s-      g = generator s-  -- 1. Sample new state.-  (!y, !q, !j) <- liftIO $ p x g-  -- 2. Calculate Metropolis-Hastings ratio.-  let !pY = pF y-      !lY = lF y-      !r = mhRatio (pX * lX) (pY * lY) q j-  -- 3. Accept or reject.-  if ln r >= 0.0-    then put $ s {item = Item y pY lY, acceptance = pushA m True a}-    else do-      b <- uniform g-      if b < exp (ln r)-        then put $ s {item = Item y pY lY, acceptance = pushA m True a}-        else put $ s {acceptance = pushA m False a}---- TODO: Splitmix. Split the generator here. See SaveSpec -> mhContinue.---- Run one iterations; perform all proposals in a Cycle.-mhIter :: ToJSON a => [Proposal a] -> Mcmc a ()-mhIter ps = do-  mapM_ mhPropose ps-  s <- get-  let i = item s-      t = trace s-      n = iteration s-  put $ s {trace = pushT i t, iteration = succ n}-  mcmcClean-  mcmcMonitorExec---- Run N iterations.-mhNIter :: ToJSON a => Int -> Mcmc a ()-mhNIter n = do-  mcmcDebugS $ "Run " <> show n <> " iterations."-  c <- gets cycle-  g <- gets generator-  cycles <- liftIO $ getNIterations c n g-  forM_ cycles mhIter---- Burn in and auto tune.-mhBurnInN :: ToJSON a => Int -> Maybe Int -> Mcmc a ()-mhBurnInN b (Just t)-  | t <= 0 = error "mhBurnInN: Auto tuning period smaller equal 0."-  | b > t = do-    mcmcResetA-    mhNIter t-    mcmcSummarizeCycle >>= mcmcDebugB-    mcmcAutotune-    mhBurnInN (b - t) (Just t)-  | otherwise = do-    mcmcResetA-    mhNIter b-    mcmcSummarizeCycle >>= mcmcInfoB-    mcmcInfoS $ "Acceptance ratios calculated over the last " <> show b <> " iterations."-mhBurnInN b Nothing = mhNIter b---- Initialize burn in for given number of iterations.-mhBurnIn :: ToJSON a => Int -> Maybe Int -> Mcmc a ()-mhBurnIn b t-  | b < 0 = error "mhBurnIn: Negative number of burn in iterations."-  | b == 0 = return ()-  | otherwise = do-    mcmcInfoS $ "Burn in for " <> show b <> " cycles."-    mcmcDebugS $ "Auto tuning period is " <> show t <> "."-    mhBurnInN b t-    mcmcInfoB "Burn in finished."---- Run for given number of iterations.-mhRun :: ToJSON a => Int -> Mcmc a ()-mhRun n = do-  mcmcResetA-  mcmcInfoS $ "Run chain for " <> show n <> " iterations."-  -- let (m, r) = n `quotRem` 100-  -- -- Print header to standard output every 100 iterations.-  -- replicateM_ m $ do-  --   mcmcMonitorStdOutHeader-  --   mhNIter 100-  -- when (r > 0) $ do-  --   mcmcMonitorStdOutHeader-  --   mhNIter r-  mhNIter n--mhT :: ToJSON a => Mcmc a ()-mhT = do-  mcmcInfoB "Metropolis-Hastings sampler."-  mcmcSummarizeCycle >>= mcmcInfoB-  mcmcReport-  s <- get-  let b = fromMaybe 0 (burnInIterations s)-  mhBurnIn b (autoTuningPeriod s)-  mhRun $ iterations s--mhContinueT :: ToJSON a => Int -> Mcmc a ()-mhContinueT dn = do-  mcmcInfoB "Continuation of Metropolis-Hastings sampler."-  mcmcInfoS $ "Run chain for " <> show dn <> " additional iterations."-  mcmcSummarizeCycle >>= mcmcInfoB-  mhRun dn---- | Continue a Markov chain for a given number of Metropolis-Hastings steps.------ At the moment, when an MCMC run is continued, the old @.mcmc@ file is--- deleted. This behavior may change in the future.------ This means that an interrupted continuation also breaks previous runs. This--- step is necessary because, otherwise, incomplete monitor files are left on--- disk, if a continuation is canceled. Subsequent continuations would append to--- the incomplete monitor files and produce garbage.-mhContinue ::-  ToJSON a =>-  -- | Additional number of Metropolis-Hastings steps.-  Int ->-  -- | Loaded status of the Markov chain.-  Status a ->-  IO (Status a)-mhContinue dn s-  | dn <= 0 = error "mhContinue: The number of iterations is zero or negative."-  | otherwise = mcmcRun (mhContinueT dn) s'-  where-    n' = iterations s + dn-    s' = s {iterations = n'}---- | Run a Markov chain for a given number of Metropolis-Hastings steps.-mh ::-  ToJSON a =>-  -- | Initial (or last) status of the Markov chain.-  Status a ->-  IO (Status a)-mh s =-  if iteration s == 0-    then mcmcRun mhT s-    else do-      putStrLn "To continue a Markov chain run, please use 'mhContinue'."-      error $ "mh: Current iteration " ++ show (iteration s) ++ " is non-zero."

@@ -16,14 +16,15 @@     Monitor (..),     MonitorStdOut,     monitorStdOut,+    msHeader,     MonitorFile,     monitorFile,     MonitorBatch,     monitorBatch,+    getMonitorBatchSize,      -- * Use monitors     mOpen,-    mAppend,     mExec,     mClose,   )@@ -33,21 +34,24 @@ import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL import Data.Int+import Data.List hiding (sum) import Data.Time.Clock+import qualified Data.Vector as VB+import Mcmc.Chain.Link+import Mcmc.Chain.Trace import Mcmc.Internal.ByteString-import Mcmc.Item import Mcmc.Monitor.Log import Mcmc.Monitor.Parameter import Mcmc.Monitor.ParameterBatch import Mcmc.Monitor.Time-import Mcmc.Trace-import Mcmc.Verbosity+import Mcmc.Settings import Numeric.Log-import System.Directory import System.IO import Prelude hiding (sum) --- | A 'Monitor' describes which part of the Markov chain should be logged and+-- | A 'Monitor' observing the chain.+--+-- A 'Monitor' describes which part of the Markov chain should be logged and -- where. Further, they allow output of summary statistics per iteration in a -- flexible way. data Monitor a = Monitor@@ -78,7 +82,7 @@   | otherwise = MonitorStdOut ps p  msIWidth :: Int-msIWidth = 12+msIWidth = 9  msWidth :: Int msWidth = 22@@ -88,6 +92,7 @@   where     vals = map (alignRight msWidth) (tail xs) +-- | Header of monitor to standard output. msHeader :: MonitorStdOut a -> BL.ByteString msHeader m = BL.intercalate "\n" [row, sep]   where@@ -96,21 +101,21 @@         ["Iteration", "Log-Prior", "Log-Likelihood", "Log-Posterior"]           ++ nms           ++ ["Runtime", "ETA"]-    sep = "   " <> BL.replicate (BL.length row - 3) '-'+    sep = BL.replicate (BL.length row) '-'     nms = [BL.pack $ mpName p | p <- msParams m]  msDataLine ::   Int ->-  Item a ->+  Link a ->   Int ->   UTCTime ->   Int ->   MonitorStdOut a ->   IO BL.ByteString-msDataLine i (Item x p l) ss st j m = do+msDataLine i (Link x p l) ss st j m = do   ct <- getCurrentTime   let dt = ct `diffUTCTime` st-      -- Careful, don't evaluate this when i == ss.+      -- 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@@ -128,7 +133,7 @@  msExec ::   Int ->-  Item a ->+  Link a ->   Int ->   UTCTime ->   Int ->@@ -136,9 +141,9 @@   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+  -- -- | 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'.@@ -165,32 +170,12 @@ mfRenderRow :: [BL.ByteString] -> BL.ByteString mfRenderRow = BL.intercalate "\t" -open' :: String -> Bool -> IO Handle-open' n frc = do-  fe <- doesFileExist n-  case (fe, frc) of-    (False, _) -> openFile n WriteMode-    (True, True) -> openFile n WriteMode-    (True, False) -> error $ "open': File \"" <> n <> "\" exists; probably use 'force'?"--mfOpen :: String -> Bool -> MonitorFile a -> IO (MonitorFile a)-mfOpen n frc m = do-  let mfn = n <> mfName m <> ".monitor"-  h <- open' mfn frc-  hSetBuffering h LineBuffering+mfOpen :: String -> String -> ExecutionMode -> MonitorFile a -> IO (MonitorFile a)+mfOpen pre suf em m = do+  let fn = intercalate "." $ filter (not . null) [pre, mfName m, suf, "monitor"]+  h <- openWithExecutionMode em fn   return $ m {mfHandle = Just h} -mfAppend :: String -> MonitorFile a -> IO (MonitorFile a)-mfAppend n m = do-  let fn = n <> mfName m <> ".monitor"-  fe <- doesFileExist fn-  if fe-    then do-      h <- openFile fn AppendMode-      hSetBuffering h LineBuffering-      return $ m {mfHandle = Just h}-    else error $ "mfAppend: Monitor file does not exist: " ++ fn ++ "."- mfHeader :: MonitorFile a -> IO () mfHeader m = case mfHandle m of   Nothing ->@@ -206,10 +191,10 @@  mfExec ::   Int ->-  Item a ->+  Link a ->   MonitorFile a ->   IO ()-mfExec i (Item x p l) m+mfExec i (Link x p l) m   | i `mod` mfPeriod m /= 0 = return ()   | otherwise = case mfHandle m of     Nothing ->@@ -231,11 +216,10 @@   Just h -> hClose h   Nothing -> error $ "mfClose: File was not opened for monitor " <> mfName m <> "." --- | Monitor to a file, but calculate batch means for the given batch size;--- constructed with 'monitorBatch'.+-- | Batch monitor to a file. ----- Batch monitors are slow at the moment because the monitored parameter has to--- be extracted from the state for each iteration.+-- Calculate summary statistics over the last given number of iterations (batch+-- size). Construct with 'monitorBatch'. data MonitorBatch a = MonitorBatch   { mbName :: String,     mbHandle :: Maybe Handle,@@ -243,12 +227,11 @@     mbSize :: Int   } --- | Monitor parameters to a file, see 'MonitorBatch'.+-- | Batch monitor parameters to a file, see 'MonitorBatch'. monitorBatch ::   -- | Name; used as part of the file name.   String ->-  -- | Instructions about which parameters to log-  -- and how to calculate the batch means.+  -- | Instructions about how to calculate the summary statistics.   [MonitorParameterBatch a] ->   -- | Batch size.   Int ->@@ -257,23 +240,17 @@   | p < 2 = error "monitorBatch: Batch size has to be 2 or larger."   | otherwise = MonitorBatch n Nothing ps p -mbOpen :: String -> Bool -> MonitorBatch a -> IO (MonitorBatch a)-mbOpen n frc m = do-  let mfn = n <> mbName m <> ".batch"-  h <- open' mfn frc-  hSetBuffering h LineBuffering-  return $ m {mbHandle = Just h}+-- | Batch monitor size.+--+-- Useful to determine the trace length.+getMonitorBatchSize :: MonitorBatch a -> Int+getMonitorBatchSize = mbSize -mbAppend :: String -> MonitorBatch a -> IO (MonitorBatch a)-mbAppend n m = do-  let fn = n <> mbName m <> ".batch"-  fe <- doesFileExist fn-  if fe-    then do-      h <- openFile fn AppendMode-      hSetBuffering h LineBuffering-      return $ m {mbHandle = Just h}-    else error $ "mbAppend: Monitor file does not exist: " ++ fn ++ "."+mbOpen :: String -> String -> ExecutionMode -> MonitorBatch a -> IO (MonitorBatch a)+mbOpen pre suf em m = do+  let fn = intercalate "." $ filter (not . null) [pre, mbName m, suf, "batch"]+  h <- openWithExecutionMode em fn+  return $ m {mbHandle = Just h}  mbHeader :: MonitorBatch a -> IO () mbHeader m = case mbHandle m of@@ -288,15 +265,15 @@         ["Iteration", "Mean log-Prior", "Mean log-Likelihood", "Mean log-Posterior"]           ++ [BL.pack $ mbpName mbp | mbp <- mbParams m] -mean :: [Log Double] -> Log Double-mean xs = sum xs / fromIntegral (length xs)+mean :: VB.Vector (Log Double) -> Log Double+mean xs = VB.sum xs / fromIntegral (VB.length xs)  mbExec ::   Int ->   Trace a ->   MonitorBatch a ->   IO ()-mbExec i t' m+mbExec i t m   | (i `mod` mbSize m /= 0) || (i == 0) = return ()   | otherwise = case mbHandle m of     Nothing ->@@ -304,22 +281,21 @@         "mbExec: No handle available for batch monitor with name "           <> mbName m           <> "."-    Just h ->+    Just h -> do+      xs <- takeT (mbSize m) t+      let lps = VB.map prior xs+          lls = VB.map likelihood xs+          los = VB.zipWith (*) lps lls+          mlps = mean lps+          mlls = mean lls+          mlos = mean los       BL.hPutStrLn h $         mfRenderRow $           BL.pack (show i) :           renderLog mlps :           renderLog mlls :           renderLog mlos :-            [BB.toLazyByteString $ mbpFunc mbp (map state t) | mbp <- mbParams m]-  where-    t = takeItems (mbSize m) t'-    lps = map prior t-    lls = map likelihood t-    los = zipWith (*) lps lls-    mlps = mean lps-    mlls = mean lls-    mlos = mean los+            [BB.toLazyByteString $ mbpFunc mbp (VB.map state xs) | mbp <- mbParams m]  mbClose :: MonitorBatch a -> IO () mbClose m = case mbHandle m of@@ -327,22 +303,22 @@   Nothing -> error $ "mfClose: File was not opened for batch monitor: " <> mbName m <> "."  -- | Open the files associated with the 'Monitor'.-mOpen :: String -> Bool -> Monitor a -> IO (Monitor a)-mOpen n frc (Monitor s fs bs) = do-  fs' <- mapM (mfOpen n frc) fs-  mapM_ mfHeader fs'-  bs' <- mapM (mbOpen n frc) bs-  mapM_ mbHeader bs'+mOpen ::+  -- Base name prefix.+  String ->+  -- Base name suffix.+  String ->+  ExecutionMode ->+  Monitor a ->+  IO (Monitor a)+mOpen pre suf em (Monitor s fs bs) = do+  fs' <- mapM (mfOpen pre suf em) fs+  unless (em == Continue) $ mapM_ mfHeader fs'+  bs' <- mapM (mbOpen pre suf em) bs+  unless (em == Continue) $ mapM_ mbHeader bs'   hSetBuffering stdout LineBuffering   return $ Monitor s fs' bs' --- | Open the files associated with the 'Monitor' in append mode.-mAppend :: String -> Monitor a -> IO (Monitor a)-mAppend n (Monitor s fs bs) = do-  fs' <- mapM (mfAppend n) fs-  bs' <- mapM (mbAppend n) bs-  return $ Monitor s fs' bs'- -- | Execute monitors; print status information to files and return text to be -- printed to standard output and log file. mExec ::@@ -362,11 +338,16 @@   Monitor a ->   IO (Maybe BL.ByteString) mExec v i ss st xs j (Monitor s fs bs) = do-  mapM_ (mfExec i $ headT xs) fs+  x <- headT xs+  mapM_ (mfExec i x) fs+  -- NOTE: Batch monitors are slow because separate batch monitors will extract+  -- separate immutable stacks from the trace. However, using folds on the+  -- mutable stack only could be an option! But then, we require two polymorphic+  -- types (for the fold).   mapM_ (mbExec i xs) bs   if v == Quiet     then return Nothing-    else msExec i (headT xs) ss st j s+    else msExec i x ss st j s  -- | Close the files associated with the 'Monitor'. mClose :: Monitor a -> IO (Monitor a)

@@ -15,7 +15,6 @@   ( -- * Parameter monitors     MonitorParameter (..),     (>$<),-    (@.),     monitorInt,     monitorDouble,     monitorDoubleF,@@ -45,19 +44,6 @@  instance Contravariant MonitorParameter where   contramap f (MonitorParameter n m) = MonitorParameter n (m . f)---- | Convert a parameter monitor from one data type to another.------ DEPRECATED.------ For example, to monitor a 'Double' value being the first entry of a tuple:------ @--- mon = fst @. monitorDouble--- @-(@.) :: (b -> a) -> MonitorParameter a -> MonitorParameter b-(@.) = contramap-{-# DEPRECATED (@.) "Superseded by the contravariant instance, use '(>$<)'." #-}  -- | Monitor 'Int'. monitorInt ::

@@ -2,7 +2,7 @@  -- | -- Module      :  Mcmc.Monitor.ParameterBatch--- Description :  Monitor parameters+-- Description :  Batch monitor parameters -- Copyright   :  (c) Dominik Schrempf, 2020 -- License     :  GPL-3.0-or-later --@@ -12,25 +12,27 @@ -- -- Creation date: Fri May 29 11:11:49 2020. ----- Batch mean monitors.+-- A batch monitor prints summary statistics of a parameter collected over a+-- specific number of last iterations. The functions provided in this module+-- calculate the mean of the monitored parameter. However, custom batch monitors+-- can use more complex functions. module Mcmc.Monitor.ParameterBatch   ( -- * Batch parameter monitors     MonitorParameterBatch (..),     (>$<),-    (@#),     monitorBatchMean,     monitorBatchMeanF,     monitorBatchMeanE,-    monitorBatchCustom,   ) 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  -- | Instruction about a parameter to monitor via batch means. Usually, the--- monitored parameter is average over the batch size. However, arbitrary+-- monitored parameter is averaged over the batch size. However, arbitrary -- functions performing more complicated analyses on the states in the batch can -- be provided. --@@ -42,35 +44,26 @@ -- mon = fst >$< monitorBatchMean -- @ ----- XXX: Batch monitors are slow at the moment because the monitored parameter--- has to be extracted from the state for each iteration.+-- Batch monitors may be slow because the monitored parameter has to be+-- extracted from the state for each iteration. data MonitorParameterBatch a = MonitorParameterBatch   { -- | Name of batch monitored parameter.     mbpName :: String,-    -- | Instruction about how to extract the batch mean from the trace.-    mbpFunc :: [a] -> BB.Builder+    -- | For a given batch, extract the summary statistics.+    mbpFunc :: VB.Vector a -> BB.Builder   } -instance Contravariant (MonitorParameterBatch) where-  contramap f (MonitorParameterBatch n m) = MonitorParameterBatch n (m . map f)---- | Convert a batch parameter monitor from one data type to another.------ For example, to batch monitor the mean of the first entry of a tuple:------ @--- mon = fst @# monitorBatchMean--- @-(@#) :: (b -> a) -> MonitorParameterBatch a -> MonitorParameterBatch b-(@#) f (MonitorParameterBatch n m) = MonitorParameterBatch n (m . map f)-{-# DEPRECATED (@#) "Superseded by the contravariant instance, use '(>$<)'." #-}+instance Contravariant MonitorParameterBatch where+  contramap f (MonitorParameterBatch n m) = MonitorParameterBatch n (m . VB.map f) -mean :: Real a => [a] -> Double-mean xs = realToFrac (sum xs) / fromIntegral (length xs)-{-# SPECIALIZE mean :: [Double] -> Double #-}-{-# SPECIALIZE mean :: [Int] -> Double #-}+mean :: Real a => VB.Vector a -> Double+mean xs = realToFrac (VB.sum xs) / fromIntegral (VB.length xs)+{-# SPECIALIZE mean :: VB.Vector Double -> Double #-}+{-# SPECIALIZE mean :: VB.Vector Int -> Double #-} --- | Batch monitor. Print the mean with eight decimal places (half precision).+-- | Batch mean monitor.+--+-- Print the mean with eight decimal places (half precision). monitorBatchMean ::   Real a =>   -- | Name.@@ -80,8 +73,10 @@ {-# SPECIALIZE monitorBatchMean :: String -> MonitorParameterBatch Int #-} {-# SPECIALIZE monitorBatchMean :: String -> MonitorParameterBatch Double #-} --- | Batch monitor. Print the mean with full precision computing the shortest--- string of digits that correctly represent the number.+-- | Batch mean monitor.+--+-- Print the mean with full precision computing the shortest string of digits+-- that correctly represent the number. monitorBatchMeanF ::   Real a =>   -- | Name.@@ -91,8 +86,10 @@ {-# SPECIALIZE monitorBatchMeanF :: String -> MonitorParameterBatch Int #-} {-# SPECIALIZE monitorBatchMeanF :: String -> MonitorParameterBatch Double #-} --- | Batch monitor real float parameters such as 'Double' with scientific--- notation and eight decimal places.+-- | Batch mean monitor.+--+-- Print the real float parameters such as 'Double' with scientific notation and+-- eight decimal places. monitorBatchMeanE ::   Real a =>   -- | Name.@@ -101,14 +98,3 @@ monitorBatchMeanE n = MonitorParameterBatch n (BB.byteString . BC.toExponential 8 . mean) {-# SPECIALIZE monitorBatchMeanE :: String -> MonitorParameterBatch Int #-} {-# SPECIALIZE monitorBatchMeanE :: String -> MonitorParameterBatch Double #-}---- | Batch monitor parameters with custom lens and builder.-monitorBatchCustom ::-  -- | Name.-  String ->-  -- | Function to calculate the batch mean.-  ([a] -> a) ->-  -- | Custom builder.-  (a -> BB.Builder) ->-  MonitorParameterBatch a-monitorBatchCustom n f b = MonitorParameterBatch n (b . f)

@@ -14,12 +14,14 @@ module Mcmc.Monitor.Time   ( renderDuration,     renderDurationS,+    renderTime,   ) where  import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL import Data.Time.Clock+import Data.Time.Format import Mcmc.Internal.ByteString  -- | Adapted from System.ProgressBar.renderDuration of package@@ -35,7 +37,7 @@     -- Total amount of seconds     ts :: Int     ts = round dt-    renderDecimal n = alignRightWith '0' 2 $ BB.toLazyByteString $ BB.intDec n+    renderDecimal n = alignRightWithNoTrim '0' 2 $ BB.toLazyByteString $ BB.intDec n  -- | Render duration in seconds. renderDurationS :: NominalDiffTime -> BL.ByteString@@ -43,3 +45,7 @@   where     ts :: Int     ts = round dt++-- | Render a time stamp.+renderTime :: FormatTime t => t -> String+renderTime = formatTime defaultTimeLocale "%B %-e, %Y, at %H:%M %P, %Z."

@@ -2,7 +2,7 @@  -- | -- Module      :  Prior--- Description :  Convenience functions to compute priors+-- Description :  Convenience functions for computing priors -- Copyright   :  (c) Dominik Schrempf, 2020 -- License     :  GPL-3.0-or-later --@@ -12,18 +12,21 @@ -- -- Creation date: Thu Jul 23 13:26:14 2020. module Mcmc.Prior-  ( -- * Continuous priors+  ( -- * Improper priors     largerThan,     positive,     lowerThan,     negative,-    uniform,-    normal,++    -- * Continuous priors     exponential,     gamma,-    -- -- * Discrete priors-    -- No discrete priors are available yet.+    normal,+    uniform, +    -- * Discrete priors+    poisson,+     -- * Auxiliary functions     product',   )@@ -36,6 +39,7 @@ import qualified Statistics.Distribution.Exponential as S import qualified Statistics.Distribution.Gamma as S import qualified Statistics.Distribution.Normal as S+import qualified Statistics.Distribution.Poisson as S  -- | Improper uniform prior; strictly larger than a given value. largerThan :: Double -> Double -> Log Double@@ -57,18 +61,27 @@ negative :: Double -> Log Double negative = lowerThan 0 --- | Uniform prior on [a, b].-uniform ::-  -- | Lower bound a.+-- | Exponential distributed prior.+exponential ::+  -- | Rate.   Double ->-  -- | Upper bound b.   Double ->+  Log Double+exponential l x = Exp $ S.logDensity d x+  where+    d = S.exponential l++-- | Gamma distributed prior.+gamma ::+  -- | Shape.   Double ->+  -- | Scale.+  Double ->+  Double ->   Log Double-uniform a b x-  | x <= a = 0-  | x >= b = 0-  | otherwise = Exp 0+gamma k t x = Exp $ S.logDensity d x+  where+    d = S.gammaDistr k t  -- | Normal distributed prior. normal ::@@ -82,27 +95,28 @@   where     d = S.normalDistr m s --- | Exponential distributed prior.-exponential ::-  -- | Rate.+-- | Uniform prior on [a, b].+uniform ::+  -- | Lower bound a.   Double ->+  -- | Upper bound b.   Double ->+  Double ->   Log Double-exponential l x = Exp $ S.logDensity d x-  where-    d = S.exponential l+uniform a b x+  | x <= a = 0+  | x >= b = 0+  | otherwise = Exp 0 --- | Gamma distributed prior.-gamma ::-  -- | Shape.-  Double ->-  -- | Scale.-  Double ->+-- | Poisson distributed prior.+poisson ::+  -- | Rate.   Double ->+  Int ->   Log Double-gamma k t x = Exp $ S.logDensity d x+poisson l x = Exp $ S.logProbability d x   where-    d = S.gammaDistr k t+    d = S.poisson l  -- | Intelligent product that stops when encountering a zero. --

@@ -1,10 +1,11 @@ {-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DerivingVia #-} {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE RankNTypes #-}  -- | -- Module      :  Mcmc.Proposal--- Description :  Proposals and cycles+-- Description :  Proposals are instruction to move around the state space -- Copyright   :  (c) Dominik Schrempf 2020 -- License     :  GPL-3.0-or-later --@@ -18,6 +19,7 @@     PName (..),     PDescription (..),     PWeight (..),+    PDimension (..),     Proposal (..),     (@~),     ProposalSimple,@@ -25,15 +27,19 @@     Tune (..),     createProposal,     tune,+    getOptimalRate,+    proposalHeader,+    proposalHLine,+    summarizeProposal,      -- * Cycle     Order (..),     Cycle (ccProposals),-    fromList,+    cycleFromList,     setOrder,-    getNIterations,+    orderProposals,     tuneCycle,-    autotuneCycle,+    autoTuneCycle,     summarizeCycle,      -- * Acceptance@@ -42,10 +48,12 @@     pushA,     resetA,     transformKeysA,-    acceptanceRatios,+    acceptanceRate,+    acceptanceRates,   ) where +import Control.DeepSeq import Data.Aeson import Data.Bifunctor import qualified Data.ByteString.Builder as BB@@ -54,7 +62,6 @@ import qualified Data.Double.Conversion.ByteString as BC import Data.Function import Data.List-import Data.Map.Strict (Map) import qualified Data.Map.Strict as M import Data.Maybe import Lens.Micro@@ -66,6 +73,7 @@ -- | Proposal name. newtype PName = PName {fromPName :: String}   deriving (Show, Eq, Ord)+  deriving (Monoid, Semigroup) via String  -- | Proposal description. newtype PDescription = PDescription {fromPDescription :: String}@@ -76,19 +84,57 @@ newtype PWeight = PWeight {fromPWeight :: Int}   deriving (Show, Eq, Ord) +-- | Proposal dimension.+--+-- The number of affected, independent parameters.+--+-- The optimal acceptance rate of low dimensional proposals is higher than for+-- high dimensional ones.+--+-- Optimal acceptance rates are still subject to controversies. As far as I+-- know, research has focused on random walk proposal with a multivariate normal+-- distribution of dimension @d@. In this case, the following acceptance rates+-- are desired:+--+-- - one dimension: 0.44 (numerical results);+--+-- - five and more dimensions: 0.234 (numerical results);+--+-- - infinite dimensions: 0.234 (theorem for specific target distributions).+--+-- See Handbook of Markov chain Monte Carlo, chapter 4.+--+-- Of course, many proposals may not be classical random walk proposals. For+-- example, the beta proposal on a simplex ('Mcmc.Proposal.Simplex.beta')+-- samples one new variable of the simplex from a beta distribution while+-- rescaling all other variables. What is the dimension of this proposal? I+-- don't know, but I set the dimension to 2. The reason is that if the dimension+-- of the simplex is 2, two variables are changed. If the dimension of the+-- simplex is high, one variable is changed substantially, while all others are+-- changed marginally.+--+-- Further, if a proposal changes a number of variables in the same way (and not+-- independently like in a random walk proposal), I still set the dimension of+-- the proposal to the number of variables changed.+--+-- Finally, I assume that proposals of unknown dimension have high dimension,+-- and use the optimal acceptance rate 0.234.+data PDimension = PDimension Int | PDimensionUnknown+ -- | A 'Proposal' is an instruction about how the Markov chain will traverse the -- state space @a@. Essentially, it is a probability mass or probability density--- conditioned on the current state (i.e., a kernel).+-- conditioned on the current state (i.e., a Markov kernel). -- -- A 'Proposal' may be tuneable in that it contains information about how to enlarge--- or shrink the step size to tune the acceptance ratio.------ No proposals with the same name and description are allowed in a 'Cycle'.+-- or shrink the step size to tune the acceptance rate. data Proposal a = Proposal   { -- | Name of the affected variable.     pName :: PName,     -- | Description of the proposal type and parameters.     pDescription :: PDescription,+    -- | Dimension of the proposal. The dimension is used to calculate the+    -- optimal acceptance rate, and does not have to be exact.+    pDimension :: PDimension,     -- | The weight determines how often a 'Proposal' is executed per iteration of     -- the Markov chain.     pWeight :: PWeight,@@ -98,10 +144,6 @@     pTuner :: Maybe (Tuner a)   } --- XXX: This should be removed.-instance Show (Proposal a) where-  show m = fromPName (pName m) <> " " <> fromPDescription (pDescription m) <> ", weight " <> show (fromPWeight $ pWeight m)- instance Eq (Proposal a) where   m == n = pName m == pName n && pDescription m == pDescription n @@ -116,15 +158,15 @@ -- scaleFirstEntryOfTuple = _1 @~ scale -- @ (@~) :: Lens' b a -> Proposal a -> Proposal b-(@~) l (Proposal n d w s t) = Proposal n d w (convertS l s) (convertT l <$> t)+(@~) l (Proposal n r d w s t) = Proposal n r d w (convertProposalSimple l s) (convertTuner l <$> t)  -- | Simple proposal without tuning information. -- -- Instruction about randomly moving from the current state to a new state, -- given some source of randomness. ----- In order to calculate the Metropolis-Hastings ratio, we need to know the--- ratio of the backward to forward kernels (i.e., the probability masses or+-- In order to calculate the Metropolis-Hastings-Green ratio, we need to know+-- the ratio of the backward to forward kernels (i.e., the probability masses or -- probability densities) and the absolute value of the determinant of the -- Jacobian matrix. --@@ -139,81 +181,82 @@ -- determinant of the Jacobian matrix differs from 1.0. type ProposalSimple a = a -> GenIO -> IO (a, Log Double, Log Double) -convertS :: Lens' b a -> ProposalSimple a -> ProposalSimple b-convertS l s = s'+convertProposalSimple :: Lens' b a -> ProposalSimple a -> ProposalSimple b+convertProposalSimple l s = s'   where     s' v g = do       (x', r, j) <- s (v ^. l) g       return (set l x' v, r, j) --- | Tune the acceptance ratio of a 'Proposal'; see 'tune', or 'autotuneCycle'.+-- | Tune the acceptance rate of a 'Proposal'; see 'tune', or 'autoTuneCycle'. data Tuner a = Tuner   { tParam :: Double,     tFunc :: Double -> ProposalSimple a   } -convertT :: Lens' b a -> Tuner a -> Tuner b-convertT l (Tuner p f) = Tuner p f'+convertTuner :: Lens' b a -> Tuner a -> Tuner b+convertTuner l (Tuner p f) = Tuner p f'   where-    f' x = convertS l $ f x+    f' x = convertProposalSimple l $ f x  -- | Tune the proposal? data Tune = Tune | NoTune   deriving (Show, Eq) --- | Create a possibly tuneable proposal.+-- | Create a tuneable proposal. createProposal ::   -- | Description of the proposal type and parameters.   PDescription ->   -- | Function creating a simple proposal for a given tuning parameter. The   -- larger the tuning parameter, the larger the proposal (and the lower the-  -- expected acceptance ratio), and vice versa.+  -- expected acceptance rate), and vice versa.   (Double -> ProposalSimple a) ->+  -- | Dimension.+  PDimension ->   -- | Name.   PName ->-  -- | PWeight.+  -- | Weight.   PWeight ->   -- | Activate tuning?   Tune ->   Proposal a-createProposal d f n w Tune = Proposal n d w (f 1.0) (Just $ Tuner 1.0 f)-createProposal d f n w NoTune = Proposal n d w (f 1.0) Nothing+createProposal r f d n w Tune = Proposal n r d w (f 1.0) (Just $ Tuner 1.0 f)+createProposal r f d n w NoTune = Proposal n r d w (f 1.0) Nothing  -- Minimal tuning parameter; subject to change. tuningParamMin :: Double tuningParamMin = 1e-12 --- | Tune a 'Proposal'. Return 'Nothing' if 'Proposal' is not tuneable. If the parameter---   @dt@ is larger than 1.0, the 'Proposal' is enlarged, if @0<dt<1.0@, it is---   shrunk. Negative tuning parameters are not allowed.-tune :: Double -> Proposal a -> Maybe (Proposal a)-tune dt m-  | dt <= 0 = error $ "tune: Tuning parameter not positive: " <> show dt <> "."-  | otherwise = do-    (Tuner t f) <- pTuner m-    -- Ensure that the tuning parameter is not too small.-    let t' = max tuningParamMin (t * dt)-    return $ m {pSimple = f t', pTuner = Just $ Tuner t' f}+-- | Tune a 'Proposal'. Return 'Nothing' if 'Proposal' is not tuneable. The size+--   of the proposal is proportional to the tuning parameter. Negative tuning+--   parameters are not allowed.+tune :: (Double -> Double) -> Proposal a -> Maybe (Proposal a)+tune f m = do+  (Tuner t g) <- pTuner m+  -- Ensure that the tuning parameter is strictly positive.+  let t' = max tuningParamMin (f t)+  return $ m {pSimple = g t', pTuner = Just $ Tuner t' g} --- The desired acceptance ratio 0.44 is optimal for one-dimensional proposals;--- one could also store the affected number of dimensions with the proposal and--- tune towards an acceptance ratio accounting for the number of dimensions.------ The optimal ratios seem to be:--- - One dimension: 0.44 (numerical result).--- - Five and more dimensions: 0.234 seems to be a good value (numerical result).--- - Infinite dimensions: 0.234 (theorem for specific target distributions).--- See Handbook of Markov chain Monte Carlo, chapter 4.-ratioOpt :: Double-ratioOpt = 0.44+-- | See 'PDimension'.+getOptimalRate :: PDimension -> Double+getOptimalRate (PDimension n)+  | n <= 0 = error "getOptimalRate: Proposal dimension is zero or negative."+  | n == 1 = 0.44+  -- Use a linear interpolation with delta 0.0515.+  | n == 2 = 0.3885+  | n == 3 = 0.337+  | n == 4 = 0.2855+  | n >= 5 = 0.234+  | otherwise = error "getOptimalRate: Proposal dimension is not an integer?"+getOptimalRate PDimensionUnknown = 0.234 --- Warn if acceptance ratio is lower.-ratioMin :: Double-ratioMin = 0.1+-- Warn if acceptance rate is lower.+rateMin :: Double+rateMin = 0.1 --- Warn if acceptance ratio is larger.-ratioMax :: Double-ratioMax = 0.9+-- Warn if acceptance rate is larger.+rateMax :: Double+rateMax = 0.9  -- | Define the order in which 'Proposal's are executed in a 'Cycle'. The total -- number of 'Proposal's per 'Cycle' may differ between 'Order's (e.g., compare@@ -241,6 +284,23 @@  instance Default Order where def = RandomO +-- Describe the order.+describeOrder :: Order -> BL.ByteString+describeOrder RandomO = "The proposals are executed in random order."+describeOrder SequentialO = "The proposals are executed sequentially."+describeOrder RandomReversibleO =+  BL.intercalate+    "\n"+    [ describeOrder RandomO,+      "A reversed copy of the shuffled proposals is appended to ensure reversibility."+    ]+describeOrder SequentialReversibleO =+  BL.intercalate+    "\n"+    [ describeOrder SequentialO,+      "A reversed copy of the sequential proposals is appended to ensure reversibility."+    ]+ -- | In brief, a 'Cycle' is a list of proposals. -- -- The state of the Markov chain will be logged only after all 'Proposal's in@@ -248,54 +308,69 @@ -- by one. The order in which the 'Proposal's are executed is specified by -- 'Order'. The default is 'RandomO'. ----- __Proposals must have unique names__, so that they can be identified.+-- No proposals with the same name and description are allowed in a 'Cycle', so+-- that they can be uniquely identified. data Cycle a = Cycle   { ccProposals :: [Proposal a],     ccOrder :: Order   }  -- | Create a 'Cycle' from a list of 'Proposal's.-fromList :: [Proposal a] -> Cycle a-fromList [] =-  error "fromList: Received an empty list but cannot create an empty Cycle."-fromList xs =+cycleFromList :: [Proposal a] -> Cycle a+cycleFromList [] =+  error "cycleFromList: Received an empty list but cannot create an empty Cycle."+cycleFromList xs =   if length (nub xs) == length xs     then Cycle xs def-    else error "fromList: Proposals are not unique."+    else error "cycleFromList: Proposals are not unique."  -- | Set the order of 'Proposal's in a 'Cycle'. setOrder :: Order -> Cycle a -> Cycle a setOrder o c = c {ccOrder = o}  -- | Replicate 'Proposal's according to their weights and possibly shuffle them.-getNIterations :: Cycle a -> Int -> GenIO -> IO [[Proposal a]]-getNIterations (Cycle xs o) n g = case o of-  RandomO -> shuffleN ps n g-  SequentialO -> return $ replicate n ps+orderProposals :: Cycle a -> GenIO -> IO [Proposal a]+orderProposals (Cycle xs o) g = case o of+  RandomO -> shuffle ps g+  SequentialO -> return ps   RandomReversibleO -> do-    psRs <- shuffleN ps n g-    return [psR ++ reverse psR | psR <- psRs]-  SequentialReversibleO -> return $ replicate n $ ps ++ reverse ps+    psR <- shuffle ps g+    return $ psR ++ reverse psR+  SequentialReversibleO -> return $ ps ++ reverse ps   where-    !ps = concat [replicate (fromPWeight $ pWeight m) m | m <- xs]+    !ps = concat [replicate (fromPWeight $ pWeight p) p | p <- xs] +-- The number of proposals depends on the order.+getNProposalsPerCycle :: Cycle a -> Int+getNProposalsPerCycle (Cycle xs o) = case o of+  RandomO -> once+  SequentialO -> once+  RandomReversibleO -> 2 * once+  SequentialReversibleO -> 2 * once+  where+    once = sum $ map (fromPWeight . pWeight) xs+ -- | Tune 'Proposal's in the 'Cycle'. See 'tune'.-tuneCycle :: Map (Proposal a) Double -> Cycle a -> Cycle a+tuneCycle :: M.Map (Proposal a) (Double -> Double) -> Cycle a -> Cycle a tuneCycle m c =   if sort (M.keys m) == sort ps     then c {ccProposals = map tuneF ps}-    else error "tuneCycle: Map contains proposals that are not in the cycle."+    else error "tuneCycle: Propoals in map and cycle do not match."   where     ps = ccProposals c     tuneF p = case m M.!? p of       Nothing -> p-      Just x -> fromMaybe p (tune x p)+      Just f -> fromMaybe p (tune f p) --- | Calculate acceptance ratios and auto tune the 'Proposal's in the 'Cycle'. For--- now, a 'Proposal' is enlarged when the acceptance ratio is above 0.44, and+-- | Calculate acceptance rates and auto tune the 'Proposal's in the 'Cycle'. For+-- now, a 'Proposal' is enlarged when the acceptance rate is above 0.44, and -- shrunk otherwise. Do not change 'Proposal's that are not tuneable.-autotuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a-autotuneCycle a = tuneCycle (M.map (\x -> exp $ x - ratioOpt) $ acceptanceRatios a)+autoTuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a+autoTuneCycle a = tuneCycle (M.mapWithKey tuningF $ acceptanceRates a)+  where+    tuningF proposal currentRate currentTuningParam =+      let optimalRate = getOptimalRate (pDimension proposal)+       in exp (currentRate - optimalRate) * currentTuningParam  renderRow ::   BL.ByteString ->@@ -306,65 +381,104 @@   BL.ByteString ->   BL.ByteString ->   BL.ByteString ->+  BL.ByteString ->   BL.ByteString-renderRow name ptype weight nAccept nReject acceptRatio tuneParam manualAdjustment = "   " <> nm <> pt <> wt <> na <> nr <> ra <> tp <> mt+renderRow name ptype weight nAccept nReject acceptRate 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 15 nAccept-    nr = alignRight 15 nReject-    ra = alignRight 15 acceptRatio+    na = alignRight 14 nAccept+    nr = alignRight 14 nReject+    ra = alignRight 14 acceptRate+    ro = alignRight 14 optimalRate     tp = alignRight 20 tuneParam     mt = alignRight 30 manualAdjustment +-- | Header of proposal summaries. proposalHeader :: BL.ByteString proposalHeader =-  renderRow "Name" "Description" "Weight" "Accepted" "Rejected" "Ratio" "Tuning parameter" "Consider manual adjustment"+  renderRow+    "Name"+    "Description"+    "Weight"+    "Accepted"+    "Rejected"+    "Rate"+    "Optimal rate"+    "Tuning parameter"+    "Consider manual adjustment" -summarizeProposal :: Proposal a -> Maybe (Int, Int, Double) -> BL.ByteString-summarizeProposal m r =+-- | Horizontal line of proposal summaries.+proposalHLine :: BL.ByteString+proposalHLine = BL.replicate (BL.length proposalHeader) '-'++-- | Proposal summary.+summarizeProposal ::+  PName ->+  PDescription ->+  PWeight ->+  -- Tuning parameter.+  Maybe Double ->+  PDimension ->+  Maybe (Int, Int, Double) ->+  BL.ByteString+summarizeProposal name description weight tuningParam dimension r =   renderRow-    (BL.pack $ fromPName $ pName m)-    (BL.pack $ fromPDescription $ pDescription m)-    weight+    (BL.pack $ fromPName name)+    (BL.pack $ fromPDescription description)+    weightStr     nAccept     nReject-    acceptRatio+    acceptRate+    optimalRate     tuneParamStr     manualAdjustmentStr   where-    weight = BB.toLazyByteString $ BB.intDec $ fromPWeight $ pWeight m+    weightStr = BB.toLazyByteString $ BB.intDec $ fromPWeight weight     nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) r     nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) r-    acceptRatio = BL.fromStrict $ maybe "" (BC.toFixed 3 . (^. _3)) r-    tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 3) (tParam <$> pTuner m)+    acceptRate = BL.fromStrict $ maybe "" (BC.toFixed 2 . (^. _3)) r+    optimalRate = BL.fromStrict $ BC.toFixed 2 $ getOptimalRate dimension+    tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 3) tuningParam     check v-      | v < ratioMin = "ratio too low"-      | v > ratioMax = "ratio too high"+      | v < rateMin = "rate too low"+      | v > rateMax = "rate too high"       | otherwise = ""     manualAdjustmentStr = BL.fromStrict $ maybe "" (check . (^. _3)) r -hLine :: BL.ByteString -> BL.ByteString-hLine s = "   " <> BL.replicate (BL.length s - 3) '-'---- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance ratios.+-- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance rates. summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> BL.ByteString summarizeCycle a c =   BL.intercalate "\n" $-    [ "Summary of proposal(s) in cycle. " <> mpi <> " proposal(s) per iteration.",+    [ "Summary of proposal(s) in cycle.",+      nProposalsFullStr,+      describeOrder (ccOrder c),       proposalHeader,-      hLine proposalHeader+      proposalHLine     ]-      ++ [summarizeProposal m (ar m) | m <- ps]-      ++ [hLine proposalHeader]+      ++ [ summarizeProposal+             (pName p)+             (pDescription p)+             (pWeight p)+             (tParam <$> pTuner p)+             (pDimension p)+             (ar p)+           | p <- ps+         ]+      ++ [proposalHLine]   where     ps = ccProposals c-    mpi = BB.toLazyByteString $ BB.intDec $ sum $ map (fromPWeight . pWeight) ps-    ar m = acceptanceRatio m a+    nProposals = getNProposalsPerCycle c+    nProposalsStr = BB.toLazyByteString $ BB.intDec nProposals+    nProposalsFullStr = case nProposals of+      1 -> nProposalsStr <> " proposal is performed per iteration."+      _ -> nProposalsStr <> " proposals are performed per iterations."+    ar m = acceptanceRate m a  -- | For each key @k@, store the number of accepted and rejected proposals.-newtype Acceptance k = Acceptance {fromAcceptance :: Map k (Int, Int)}+newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k (Int, Int)}+  deriving (Eq, Read, Show)  instance ToJSONKey k => ToJSON (Acceptance k) where   toJSON (Acceptance m) = toJSON m@@ -380,16 +494,16 @@ emptyA ks = Acceptance $ M.fromList [(k, (0, 0)) | k <- ks]  -- | For key @k@, prepend an accepted (True) or rejected (False) proposal.-pushA :: (Ord k, Show k) => k -> Bool -> Acceptance k -> Acceptance k-pushA k True = Acceptance . M.adjust (first succ) k . fromAcceptance-pushA k False = Acceptance . M.adjust (second succ) k . fromAcceptance+pushA :: Ord k => k -> Bool -> Acceptance k -> Acceptance k+pushA k True = Acceptance . M.adjust (force . first succ) k . fromAcceptance+pushA k False = Acceptance . M.adjust (force . second succ) k . fromAcceptance {-# INLINEABLE pushA #-}  -- | Reset acceptance storage. resetA :: Ord k => Acceptance k -> Acceptance k resetA = emptyA . M.keys . fromAcceptance -transformKeys :: (Ord k1, Ord k2) => [k1] -> [k2] -> Map k1 v -> Map k2 v+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   where     insrt m' (k1, k2) = M.insert k2 (m M.! k1) m'@@ -399,13 +513,13 @@ transformKeysA :: (Ord k1, Ord k2) => [k1] -> [k2] -> Acceptance k1 -> Acceptance k2 transformKeysA ks1 ks2 = Acceptance . transformKeys ks1 ks2 . fromAcceptance --- | Acceptance counts and ratio for a specific proposal.-acceptanceRatio :: (Show k, Ord k) => k -> Acceptance k -> Maybe (Int, Int, Double)-acceptanceRatio k a = case fromAcceptance a M.!? k of+-- | Acceptance counts and rate for a specific proposal.+acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, Double)+acceptanceRate k a = case fromAcceptance a M.!? k of   Just (0, 0) -> Nothing   Just (as, rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))-  Nothing -> error $ "acceptanceRatio: Key not found in map: " ++ show k ++ "."+  Nothing -> error "acceptanceRate: Key not found in map." --- | Acceptance ratios for all proposals.-acceptanceRatios :: Acceptance k -> Map k Double-acceptanceRatios = M.map (\(as, rs) -> fromIntegral as / fromIntegral (as + rs)) . fromAcceptance+-- | Acceptance rates for all proposals.+acceptanceRates :: Acceptance k -> M.Map k Double+acceptanceRates = M.map (\(as, rs) -> fromIntegral as / fromIntegral (as + rs)) . fromAcceptance

@@ -75,24 +75,24 @@ -- The [Bactrian -- kernel](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/figure/fig01) is -- a mixture of two symmetrically arranged normal distributions. The spike--- parameter (0, 1) loosely determines the standard deviations (>0.0) of the--- individual humps while the second parameter refers to the standard deviation--- of the complete Bactrian kernel.+-- parameter \(m \in (0, 1)\) loosely determines the standard deviations of the+-- individual humps while the second parameter \(s > 0\) refers to the+-- standard deviation of the complete Bactrian kernel. -- -- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/. slideBactrian ::-  -- | Spike parameter.+  -- | Spike parameter \(m\).   Double ->-  -- | Standard deviation.+  -- | Standard deviation \(s\).   Double ->   -- | Name.   PName ->-  -- | PWeight.+  -- | Weight.   PWeight ->   -- | Enable tuning.   Tune ->   Proposal Double-slideBactrian m s = createProposal description (bactrianAdditiveSimple m s)+slideBactrian m s = createProposal description (bactrianAdditiveSimple m s) (PDimension 1)   where     description = PDescription $ "Slide Bactrian; spike: " ++ show m ++ ", sd: " ++ show s @@ -135,11 +135,11 @@   Double ->   -- | Name.   PName ->-  -- | PWeight.+  -- | Weight.   PWeight ->   -- | Enable tuning.   Tune ->   Proposal Double-scaleBactrian m s = createProposal description (bactrianMultSimple m s)+scaleBactrian m s = createProposal description (bactrianMultSimple m s) (PDimension 1)   where     description = PDescription $ "Scale Bactrian; spike: " ++ show m <> ", sd: " <> show s

@@ -1,6 +1,6 @@ -- | -- Module      :  Mcmc.Proposal.Generic--- Description :  Generic interface to create proposals+-- Description :  Generic interface for creating proposals -- Copyright   :  (c) Dominik Schrempf 2020 -- License     :  GPL-3.0-or-later --

@@ -2,7 +2,7 @@  -- | -- Module      :  Mcmc.Proposal.Scale--- Description :  Scaling proposal with Gamma distribution+-- Description :  Multiplicative proposals -- Copyright   :  (c) Dominik Schrempf 2020 -- License     :  GPL-3.0-or-later --@@ -35,38 +35,32 @@   where     jac _ = Exp . log . recip --- | Multiplicative proposal with Gamma distributed kernel.+-- | Multiplicative proposal with gamma distributed kernel. scale ::   -- | Shape.   Double ->   -- | Scale.   Double ->-  -- | Name.   PName ->-  -- | PWeight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal Double-scale k th = createProposal description (scaleSimple k th)+scale k th = createProposal description (scaleSimple k th) (PDimension 1)   where     description = PDescription $ "Scale; shape: " ++ show k ++ ", scale: " ++ show th --- | Multiplicative proposal with Gamma distributed kernel.+-- | Multiplicative proposal with gamma distributed kernel. ----- The scale of the Gamma distributions is set to (shape)^{-1}, so that the mean--- of the Gamma distribution is 1.0.+-- The scale of the gamma distribution is set to (shape)^{-1}, so that the mean+-- of the gamma distribution is 1.0. scaleUnbiased ::   -- | Shape.   Double ->-  -- | Name.   PName ->-  -- | PWeight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal Double-scaleUnbiased k = createProposal description (scaleSimple k (1 / k))+scaleUnbiased k = createProposal description (scaleSimple k (1 / k)) (PDimension 1)   where     description = PDescription $ "Scale unbiased; shape: " ++ show k @@ -81,11 +75,7 @@     contra (x, y) u = (x * u, y / u)     jac _ u = Exp $ log $ recip $ u * u --- -- Determinant of Jacobian matrix.--- contraJac :: (Double, Double) -> Double--- contraJac (x, y) = x * y---- | Multiplicative proposal with Gamma distributed kernel.+-- | Multiplicative proposal with gamma distributed kernel. -- -- The two values are scaled contrarily so that their product stays constant. -- Contrary proposals are useful when parameters are confounded.@@ -94,13 +84,10 @@   Double ->   -- | Scale.   Double ->-  -- | Name.   PName ->-  -- | PWeight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal (Double, Double)-scaleContrarily k th = createProposal description (scaleContrarilySimple k th)+scaleContrarily k th = createProposal description (scaleContrarilySimple k th) (PDimension 2)   where     description = PDescription $ "Scale contrariliy; shape: " ++ show k ++ ", scale: " ++ show th

@@ -88,7 +88,7 @@     t' = t / 10000     -- Extremely small tuning parameters lead to numeric overflow. The square     -- root pulls the tuning parameter closer to 1.0. However, overflow may-    -- still occur (the involved Gamma functions grow faster than the+    -- still occur (the involved gamma functions grow faster than the     -- exponential). I did not observe numeric underflow in my tests.     t'' = sqrt t' @@ -126,15 +126,18 @@ -- -- For a given element of a K-dimensional simplex, propose a new element of the -- K-dimensional simplex. The new element is sampled from the multivariate--- Dirichlet distribution with parameter vector being the old element of the--- simplex.+-- Dirichlet distribution with parameter vector being proportional to the+-- current element of the simplex. -- -- The tuning parameter is used to determine the concentration of the Dirichlet -- distribution: the lower the tuning parameter, the higher the concentration. ----- This proposal may have low acceptance ratios. In this case, please see the--- coordinate wise 'beta' proposal.-dirichlet :: PName -> PWeight -> Tune -> Proposal Simplex+-- The proposal dimension, which is the dimension of the simplex, is used to+-- determine the optimal acceptance rate.+--+-- For high dimensional simplices, this proposal may have low acceptance ratios.+-- In this case, please see the coordinate wise 'beta' proposal.+dirichlet :: PDimension -> PName -> PWeight -> Tune -> Proposal Simplex dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple  -- The tuning parameter is the inverted mean of the shape values.@@ -197,7 +200,10 @@ -- No "out of bounds" checks are performed during compile time. Run time errors -- can occur if @i@ is negative, or if @i-1@ is larger than the length of the -- element vector of the simplex.+--+-- This proposal has been assigned a dimension of 2. See the discussion at+-- 'PDimension'. beta :: Int -> PName -> PWeight -> Tune -> Proposal Simplex-beta i = createProposal description (betaSimple i)+beta i = createProposal description (betaSimple i) (PDimension 2)   where     description = PDescription $ "Beta; coordinate: " ++ show i

@@ -2,7 +2,7 @@  -- | -- Module      :  Mcmc.Proposal.Slide--- Description :  Normally distributed proposal+-- Description :  Additive proposals -- Copyright   :  (c) Dominik Schrempf 2020 -- License     :  GPL-3.0-or-later --@@ -37,12 +37,12 @@   Double ->   -- | Name.   PName ->-  -- | PWeight.+  -- | Weight.   PWeight ->   -- | Enable tuning.   Tune ->   Proposal Double-slide m s = createProposal description (slideSimple m s)+slide m s = createProposal description (slideSimple m s) (PDimension 1)   where     description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s @@ -52,19 +52,19 @@   genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing Nothing  -- | Additive proposal with normally distributed kernel with mean zero. This--- proposal is very fast, because the Metropolis-Hastings ratio does not include--- calculation of the forwards and backwards kernels.+-- proposal is very fast, because the Metropolis-Hastings-Green ratio does not+-- include calculation of the forwards and backwards kernels. slideSymmetric ::   -- | Standard deviation.   Double ->   -- | Name.   PName ->-  -- | PWeight.+  -- | Weight.   PWeight ->   -- | Enable tuning.   Tune ->   Proposal Double-slideSymmetric s = createProposal description (slideSymmetricSimple s)+slideSymmetric s = createProposal description (slideSymmetricSimple s) (PDimension 1)   where     description = PDescription $ "Slide symmetric; sd: " ++ show s @@ -74,19 +74,19 @@   genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing Nothing  -- | Additive proposal with uniformly distributed kernel with mean zero. This--- proposal is very fast, because the Metropolis-Hastings ratio does not include--- calculation of the forwards and backwards kernels.+-- proposal is very fast, because the Metropolis-Hastings-Green ratio does not+-- include calculation of the forwards and backwards kernels. slideUniformSymmetric ::   -- | Delta.   Double ->   -- | Name.   PName ->-  -- | PWeight.+  -- | Weight.   PWeight ->   -- | Enable tuning.   Tune ->   Proposal Double-slideUniformSymmetric d = createProposal description (slideUniformSimple d)+slideUniformSymmetric d = createProposal description (slideUniformSimple d) (PDimension 1)   where     description = PDescription $ "Slide uniform symmetric; delta: " ++ show d @@ -108,11 +108,11 @@   Double ->   -- | Name.   PName ->-  -- | PWeight.+  -- | Weight.   PWeight ->   -- | Enable tuning.   Tune ->   Proposal (Double, Double)-slideContrarily m s = createProposal description (slideContrarilySimple m s)+slideContrarily m s = createProposal description (slideContrarilySimple m s) (PDimension 2)   where     description = PDescription $ "Slide contrarily; mean: " ++ show m ++ ", sd: " ++ show s

@@ -1,180 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE TemplateHaskell #-}---- |--- Module      :  Mcmc.Save--- Description :  Save the state of a Markov chain--- Copyright   :  (c) Dominik Schrempf, 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Tue Jun 16 10:18:54 2020.------ Save and load an MCMC run. 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.Save-  ( saveStatus,-    loadStatus,-  )-where--import Codec.Compression.GZip-import Control.Monad-import Data.Aeson-import Data.Aeson.TH-import qualified Data.ByteString.Lazy.Char8 as BL-import Data.List hiding (cycle)-import qualified Data.Map as M-import Data.Maybe-import Data.Vector.Unboxed (Vector)-import Data.Word--- TODO: Splitmix. Reproposal as soon as split mix is used and is available with the--- statistics package.-import Mcmc.Item-import Mcmc.Monitor-import Mcmc.Proposal-import Mcmc.Status hiding (save)-import Mcmc.Trace-import Mcmc.Verbosity-import Numeric.Log-import System.Directory-import System.IO.Unsafe (unsafePerformIO)-import System.Random.MWC-import Prelude hiding (cycle)--data Save a-  = Save-      -- Variables related to the chain.-      String -- Name.-      (Item a)-      Int -- Iteration.-      (Trace a)-      (Acceptance Int)-      (Maybe Int) -- Burn in.-      (Maybe Int) -- Auto tune.-      Int -- Iterations.-      Bool -- Force.-      (Maybe Int) -- Save.-      Verbosity-      (Vector Word32) -- Current seed.--      -- Variables related to the algorithm.-      [Maybe Double] -- Tuning parameters.--$(deriveJSON defaultOptions ''Save)--toSave :: Status a -> Save a-toSave (Status nm it i tr ac br at is f sv vb g _ _ _ _ _ c _) =-  Save-    nm-    it-    i-    tr'-    ac'-    br-    at-    is-    f-    sv-    vb-    g'-    ts-  where-    tr' = takeT (fromMaybe 0 sv) tr-    ac' = transformKeysA (ccProposals c) [0 ..] ac-    -- TODO: Splitmix. Remove as soon as split mix is used and is available with-    -- the statistics package.-    g' = fromSeed $ unsafePerformIO $ save g-    ts = [fmap tParam mt | mt <- map pTuner $ ccProposals c]---- | Save a 'Status' to file.------ Some important values have to be provided upon restoring the status. See--- 'loadStatus'.-saveStatus :: ToJSON a => FilePath -> Status a -> IO ()-saveStatus fn s = BL.writeFile fn $ compress $ encode (toSave s)---- fromSav prior lh cycle monitor save-fromSave ::-  (a -> Log Double) ->-  (a -> Log Double) ->-  Cycle a ->-  Monitor a ->-  Maybe (Cleaner a) ->-  Save a ->-  Status a-fromSave pr lh cc m cl (Save nm it i tr ac' br at is f sv vb g' ts) =-  Status-    nm-    it-    i-    tr-    ac-    br-    at-    is-    f-    sv-    vb-    g-    Nothing-    Nothing-    pr-    lh-    cl-    cc'-    m-  where-    ac = transformKeysA [0 ..] (ccProposals cc) ac'-    -- TODO: Splitmix. Remove as soon as split mix is used and is available with-    -- the statistics package.-    g = unsafePerformIO $ restore $ toSeed g'-    cc' = tuneCycle (M.mapMaybe id $ M.fromList $ zip (ccProposals cc) ts) cc---- | Load a 'Status' from file.------ Important information that cannot be saved and has to be provided again when--- a chain is restored:--- - prior function--- - likelihood function--- - cleaning function--- - cycle--- - monitor------ To avoid incomplete continued runs, the @.mcmc@ file is removed after load.-loadStatus ::-  FromJSON a =>-  -- | Prior function.-  (a -> Log Double) ->-  -- | Likelihood function.-  (a -> Log Double) ->-  Cycle a ->-  Monitor a ->-  -- | Cleaner, if needed.-  Maybe (Cleaner a) ->-  -- | Path of status to load.-  FilePath ->-  IO (Status a)-loadStatus pr lh cc mn cl fn = do-  res <- eitherDecode . decompress <$> BL.readFile fn-  let s = case res of-        Left err -> error err-        Right sv -> fromSave pr lh cc mn cl sv-  -- Check if prior and likelihood matches.-  let Item x svp svl = item s-  -- Recompute and check the prior and likelihood for the last state because the-  -- functions may have changed. Of course, we cannot test for the same-  -- function, but having the same prior and likelihood at the last state is-  -- already a good indicator.-  when-    (pr x /= svp)-    (error "loadStatus: Provided prior function does not match the saved prior.")-  when-    (lh x /= svl)-    (error "loadStatus: Provided likelihood function does not match the saved likelihood.")-  removeFile fn-  return s

@@ -0,0 +1,236 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TemplateHaskell #-}++-- |+-- Module      :  Mcmc.Settings+-- Description :  Settings of Markov chain Monte Carlo samplers+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Mon Nov 16 11:13:01 2020.+module Mcmc.Settings+  ( -- * Data types+    AnalysisName (..),+    BurnInSpecification (..),+    burnInIterations,+    Iterations (..),+    ExecutionMode (..),+    openWithExecutionMode,+    ParallelizationMode (..),+    SaveMode (..),+    Verbosity (..),++    -- * Settings+    Settings (..),+    settingsSave,+    settingsLoad,+    settingsCheck,+  )+where++import Data.Aeson+import Data.Aeson.TH+import qualified Data.ByteString.Lazy.Char8 as BL+import System.Directory+import System.IO++-- | Analysis name of the MCMC sampler.+newtype AnalysisName = AnalysisName {fromAnalysisName :: String}+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''AnalysisName)++-- | Burn in specification.+data BurnInSpecification+  = -- | No burn in.+    NoBurnIn+  | -- | Burn in for a given number of iterations.+    BurnInWithoutAutoTuning Int+  | -- | Burn in for a given number of iterations. Enable auto tuning with a+    -- given period.+    BurnInWithAutoTuning Int Int+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''BurnInSpecification)++-- | Get the number of burn in iterations.+burnInIterations :: BurnInSpecification -> Int+burnInIterations NoBurnIn = 0+burnInIterations (BurnInWithoutAutoTuning n) = n+burnInIterations (BurnInWithAutoTuning n _) = n++-- | Number of normal iterations after burn in.+--+-- Note that auto tuning only happens during burn in.+newtype Iterations = Iterations {fromIterations :: Int}+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''Iterations)++-- | Execution mode.+data ExecutionMode+  = -- | Perform new run.+    --+    -- Call 'error' if an output files exists.+    Fail+  | -- | Perform new run.+    --+    -- Overwrite existing output files.+    Overwrite+  | -- | Continue a previous run and append to output files.+    --+    -- Call 'error' if an output file does not exist.+    Continue+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''ExecutionMode)++-- | Open a file honoring the execution mode.+--+-- Call 'error' if execution mode is+--+-- - 'Continue' and file does not exist.+--+-- - 'Fail' and file exists.+openWithExecutionMode :: ExecutionMode -> FilePath -> IO Handle+openWithExecutionMode em fn = do+  fe <- doesFileExist fn+  case (em, fe) of+    (Continue, False) ->+      error $ "openWithExecutionMode: Cannot continue; file does not exist: " ++ fn ++ "."+    (Continue, True) ->+      openFile fn AppendMode+    (Fail, True) ->+      error $ "openWithExecutionMode: File exists: " ++ fn ++ "; use 'Overwrite'?"+    _ -> do+      h <- openFile fn WriteMode+      hSetBuffering h LineBuffering+      return h++-- One could automatically select 'Parallel' or 'Sequential' according to the+-- number of capabilities when initializing the environment or according to the+-- iteration time in dependence of the number of used capabilities. However, I+-- decided to opt for a manual configuration, because more capabilities may be+-- available and other parts of the program may be executed in parallel even if+-- sequential execution of the MCMC sampler is beneficial.++-- | Parallelization mode.+--+-- Parallel execution of the chains is only beneficial when the algorithm allows+-- for parallelization, and if computation of the next iteration takes a long+-- time. If the calculation of the next state is fast, sequential execution is+-- usually beneficial, even for algorithms involving parallel chains. If the+-- calculation of the next state is slow, parallel execution may be beneficial.+--+-- - The "Mcmc.Algorithm.Metropolis" algorithm is inherently sequential.+--+-- - The "Mcmc.Algorithm.MC3" algorithm works well with parallelization.+--+-- Of course, also the prior or likelihood functions can be computed in+-- parallel. However, this library is not aware of how these functions are+-- computed.+data ParallelizationMode+  = Sequential+  | Parallel+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''ParallelizationMode)++-- | Should the MCMC run be saved at the end of the run?+data SaveMode = NoSave | Save+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''SaveMode)++-- $(deriveJSON defaultOptions ''CleaningMode)++-- | Not much to say here.+data Verbosity = Quiet | Warn | Info | Debug+  deriving (Eq, Ord, Read, Show)++$(deriveJSON defaultOptions ''Verbosity)++-- | Settings of an MCMC sampler.+data Settings = Settings+  { sAnalysisName :: AnalysisName,+    sBurnIn :: BurnInSpecification,+    sIterations :: Iterations,+    sExecutionMode :: ExecutionMode,+    sParallelizationMode :: ParallelizationMode,+    sSaveMode :: SaveMode,+    sVerbosity :: Verbosity+  }+  deriving (Eq, Show)++$(deriveJSON defaultOptions ''Settings)++settingsFn :: String -> FilePath+settingsFn n = n ++ ".settings"++-- | Save settings to a file determined by the analysis name.+settingsSave :: Settings -> IO ()+settingsSave s = BL.writeFile fn $ encode s+  where+    fn = settingsFn $ fromAnalysisName $ sAnalysisName s++-- | Load settings.+settingsLoad :: AnalysisName -> IO Settings+settingsLoad (AnalysisName n) = either error id . eitherDecode <$> BL.readFile fn+  where+    fn = settingsFn n++-- Show settings and call 'error'.+settingsError :: Settings -> Int -> String -> a+settingsError s i err =+  error $+    show s+      ++ "\n"+      ++ "Current iteration: "+      ++ show i+      ++ "\n"+      ++ "settingsError: "+      ++ err++-- | Check settings.+--+-- Call 'error' if:+--+-- - The analysis name is the empty string.+--+-- - The number of burn in iterations is negative.+--+-- - Auto tuning period is zero or negative.+--+-- - The number of iterations is negative.+--+-- - The current iteration is larger than the total number of iterations.+--+-- - The current iteration is non-zero but the execution mode is not 'Continue'.+--+-- - The current iteration is zero but the execution mode is 'Continue'.+settingsCheck ::+  Settings ->+  -- | Current iteration.+  Int ->+  IO ()+settingsCheck s@(Settings nm bi i em _ _ _) iCurrent+  | null (fromAnalysisName nm) = serr "Analysis name is the empty string."+  | burnInIterations bi < 0 = serr "Number of burn in iterations is negative."+  | not $ burnInAutoTuningPeriodValid bi = serr "Auto tuning period is zero or negative."+  | fromIterations i < 0 = serr "Number of iterations is negative."+  | burnInIterations bi + fromIterations i - iCurrent < 0 =+    serr "Current iteration is larger than the total number of iterations."+  | iCurrent /= 0 && em /= Continue =+    serr "Current iteration is non-zero but execution mode is not 'Continue'."+  | iCurrent == 0 && em == Continue =+    serr "Current iteration is zero but execution mode is 'Continue'."+  | otherwise = return ()+  where+    serr = settingsError s iCurrent+    burnInAutoTuningPeriodValid :: BurnInSpecification -> Bool+    burnInAutoTuningPeriodValid (BurnInWithAutoTuning _ t) = t > 0+    burnInAutoTuningPeriodValid _ = True

@@ -1,218 +0,0 @@--- Note: It is not necessary to add another type @b@ to store supplementary--- information about the chain. The information can just be stored in @a@--- equally well.---- XXX: Status tuned exclusively to the Metropolis-Hastings algorithm. We should--- abstract the algorithm from the chain. Maybe something like:------ @--- data Status a = Status { Chain a; Algorithm a}--- @---- |--- Module      :  Mcmc.Status--- Description :  What is an MCMC?--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Tue May  5 18:01:15 2020.-module Mcmc.Status-  ( Cleaner (..),-    Status (..),-    status,-    cleanWith,-    saveWith,-    force,-    quiet,-    debug,-    noData,-  )-where--import Data.Maybe-import Data.Time.Clock-import Mcmc.Item-import Mcmc.Monitor-import Mcmc.Proposal-import Mcmc.Trace-import Mcmc.Verbosity (Verbosity (..))-import Numeric.Log-import System.IO-import System.Random.MWC hiding (save)-import Prelude hiding (cycle)---- | Clean the state periodically.------ The prior and the likelihood will be updated after the cleaning process.------ For long chains, successive numerical errors can accumulate such that the--- state diverges from honoring specific required constraints. In these cases, a--- 'Cleaner' can be used to ensure that the required constraints of the state--- are honored. For example, the branches of an ultrametric phylogeny may--- diverge slightly after successful many proposals such that the phylogeny is--- not anymore ultrametric.------ Please be aware that the Markov chain will not converge to the true posterior--- distribution if the state is changed substantially! Only apply subtle changes--- that are absolutely necessary to preserve the required properties of the--- state such as specific numerical constraints.-data Cleaner a = Cleaner-  { -- | Clean every given number of iterations.-    clEvery :: Int,-    -- | Cleaning function. Executed before monitoring the state.-    clFunction :: a -> a-  }---- | The 'Status' contains all information to run an MCMC chain. It is--- constructed using the function 'status'.------ The polymorphic type @a@ stores the state of the chain. It can also be used--- to store auxiliary information.-data Status a = Status-  { -- MCMC related variables; saved.--    -- | The name of the MCMC chain; used as file prefix.-    name :: String,-    -- | The current 'Item' of the chain combines the current state and the-    -- current likelihood.-    item :: Item a,-    -- | The iteration is the number of completed cycles.-    iteration :: Int,-    -- | The 'Trace' of the Markov chain in reverse order, the most recent-    -- 'Item' is at the head of the list.-    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),-    -- | Number of burn in iterations; deactivate burn in with 'Nothing'.-    burnInIterations :: Maybe Int,-    -- | Auto tuning period (only during burn in); deactivate auto tuning with-    -- 'Nothing'.-    autoTuningPeriod :: Maybe Int,-    -- | Number of normal iterations excluding burn in. Note that auto tuning-    -- only happens during burn in.-    iterations :: Int,-    ---    -- Auxiliary variables; saved.--    -- | Overwrite output files? Default is 'False', change with 'force'.-    forceOverwrite :: Bool,-    -- | Save the chain with trace of given length at the end of the run?-    -- Default is no save ('Nothing'). Change with 'saveWith'.-    save :: Maybe Int,-    -- | Verbosity.-    verbosity :: Verbosity,-    -- | The random number generator.-    generator :: GenIO,-    ---    -- Auxiliary variables; not saved.--    -- | Starting time and starting iteration of chain; used to calculate-    -- run time and ETA.-    start :: Maybe (Int, UTCTime),-    -- | Handle to log file.-    logHandle :: Maybe Handle,-    ---    -- Auxiliary functions; not saved.--    -- | The prior function. The un-normalized posterior is the product of the-    -- prior and the likelihood.-    priorF :: a -> Log Double,-    -- | The likelihood function. The un-normalized posterior is the product of-    -- the prior and the likelihood.-    likelihoodF :: a -> Log Double,-    -- | Clean the state periodically.-    cleaner :: Maybe (Cleaner a),-    ---    -- Variables related to the algorithm; not saved.--    -- | A set of 'Proposal's form a 'Cycle'.-    cycle :: Cycle a,-    -- | A 'Monitor' observing the chain.-    monitor :: Monitor a-  }---- | Initialize the 'Status' of a Markov chain Monte Carlo run.-status ::-  -- | Name of the Markov chain; used as file prefix.-  String ->-  -- | The prior function.-  (a -> Log Double) ->-  -- | The likelihood function.-  (a -> Log Double) ->-  -- | A list of 'Proposal's executed in forward order. The chain will be logged-  -- after each cycle.-  Cycle a ->-  -- | A 'Monitor' observing the chain.-  Monitor a ->-  -- | The initial state in the state space @a@.-  a ->-  -- | Number of burn in iterations; deactivate burn in with 'Nothing'.-  Maybe Int ->-  -- | Auto tuning period (only during burn in); deactivate auto tuning with-  -- 'Nothing'.-  Maybe Int ->-  -- | Number of normal iterations excluding burn in. Note that auto tuning only-  -- happens during burn in.-  Int ->-  -- | A source of randomness. For reproducible runs, make sure to use-  -- generators with the same, fixed seed.-  GenIO ->-  Status a-status n p l c m x mB mT nI g-  | isJust mT && isNothing mB = error "status: Auto tuning period given, but no burn in."-  | otherwise =-    Status-      n-      i-      0-      (singletonT i)-      (emptyA $ ccProposals c)-      mB-      mT-      nI-      False-      Nothing-      Info-      g-      Nothing-      Nothing-      p-      l-      Nothing-      c-      m-  where-    i = Item x (p x) (l x)---- | Clean the state every given number of generations using the given function.--- See 'Cleaner'.-cleanWith :: Cleaner a -> Status a -> Status a-cleanWith c s = s {cleaner = Just c}---- | Save the Markov chain with trace of given length.-saveWith :: Int -> Status a -> Status a-saveWith n s = s {save = Just n}---- | Overwrite existing files; it is not necessary to use 'force', when a chain--- is continued.-force :: Status a -> Status a-force s = s {forceOverwrite = True}---- | Do not print anything to standard output. Do not create log file. File--- monitors and batch monitors are executed normally.-quiet :: Status a -> Status a-quiet s = s {verbosity = Quiet}---- | Be verbose.-debug :: Status a -> Status a-debug s = s {verbosity = Debug}---- | Set the likelihood function to 1.0. Useful for debugging and testing.-noData :: Status a -> Status a-noData s = s {likelihoodF = const 1.0}

@@ -1,63 +0,0 @@--- |--- Module      :  Mcmc.Trace--- Description :  Trace of a Markov chain--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Wed May 20 09:11:25 2020.-module Mcmc.Trace-  ( Trace,-    singletonT,-    pushT,-    headT,-    takeItems,-    takeT,-  )-where--import Data.Aeson-import Mcmc.Item---- | A 'Trace' passes through a list of states with associated likelihoods which--- are called 'Item's. New 'Item's are prepended, and the path of the Markov--- chain is stored in reversed order.-newtype Trace a = Trace {fromTrace :: [Item a]}-  deriving (Show, Read, Eq)--instance Semigroup (Trace a) where-  (Trace l) <> (Trace r) = Trace (l <> r)--instance Monoid (Trace a) where-  mempty = Trace []--instance ToJSON a => ToJSON (Trace a) where-  toJSON (Trace xs) = toJSON xs-  toEncoding (Trace xs) = toEncoding xs--instance FromJSON a => FromJSON (Trace a) where-  parseJSON v = Trace <$> parseJSONList v---- | The empty trace.-singletonT :: Item a -> Trace a-singletonT i = Trace [i]---- | Prepend an 'Item' to a 'Trace'.-pushT :: Item a -> Trace a -> Trace a-pushT x = Trace . (:) x . fromTrace-{-# INLINEABLE pushT #-}---- | Get the most recent item of the trace.-headT :: Trace a -> Item a-headT = head . fromTrace---- | Get the N most recent items of the trace.-takeItems :: Int -> Trace a -> [Item a]-takeItems n = take n . fromTrace---- | Shorten the trace to given length.-takeT :: Int -> Trace a -> Trace a-takeT n = Trace . take n . fromTrace

@@ -1,40 +0,0 @@-{-# LANGUAGE TemplateHaskell #-}---- |--- Module      :  Mcmc.Verbosity--- Description :  Be quiet! Or better not?--- Copyright   :  (c) Dominik Schrempf, 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Sat Jun 27 10:49:28 2020.-module Mcmc.Verbosity-  ( Verbosity (..),-    warn,-    info,-    debug,-  )-where--import Control.Monad-import Data.Aeson.TH---- | Not much to say here.-data Verbosity = Quiet | Warn | Info | Debug deriving (Show, Eq, Ord)--$(deriveJSON defaultOptions ''Verbosity)---- | Perform action if 'Verbosity' is 'Warn' or higher.-warn :: Applicative m => Verbosity -> m () -> m ()-warn v = when (v >= Warn)---- | Perform action if 'Verbosity' is 'Info' or higher.-info :: Applicative m => Verbosity -> m () -> m ()-info v = when (v >= Info)---- | Perform action if 'Verbosity' is 'Debug'.-debug :: Applicative m => Verbosity -> m () -> m ()-debug v = when (v == Debug)

@@ -26,18 +26,18 @@ p2 = slideSymmetric 1.0 (PName "Test 2") (PWeight 3) Tune  c :: Cycle Double-c = fromList [p1, p2]+c = cycleFromList [p1, p2]  spec :: Spec spec =-  describe "getNIterations" $+  describe "orderProposals" $     it "returns the correct number of proposals in a cycle" $       do         g <- create-        l1 <- length . head <$> getNIterations c 1 g+        l1 <- length <$> orderProposals c g         l1 `shouldBe` 4-        l2 <- length . head <$> getNIterations (setOrder RandomReversibleO c) 1 g+        l2 <- length <$> orderProposals (setOrder RandomReversibleO c) g         l2 `shouldBe` 8-        o3 <- head <$> getNIterations (setOrder SequentialReversibleO c) 1 g+        o3 <- orderProposals (setOrder SequentialReversibleO c) g         length o3 `shouldBe` 8-        o3 `shouldBe` [p1, p2, p2, p2, p2, p2, p2, p1]+        o3 == [p1, p2, p2, p2, p2, p2, p2, p1] `shouldBe` True

@@ -15,15 +15,13 @@ where  import Mcmc-import Mcmc.Save-import Mcmc.Status hiding (save)+import Mcmc.Chain.Chain+import Mcmc.Chain.Save+import Mcmc.Chain.Trace import Numeric.Log-import Statistics.Distribution hiding-  ( mean,-    stdDev,-  )+import Statistics.Distribution import Statistics.Distribution.Normal-import System.Random.MWC+import qualified System.Random.MWC as R import Test.Hspec  trueMean :: Double@@ -32,12 +30,12 @@ trueStdDev :: Double trueStdDev = 4 -lh :: Double -> Log Double+lh :: LikelihoodFunction Double lh = Exp . logDensity (normalDistr trueMean trueStdDev)  proposals :: Cycle Double proposals =-  fromList+  cycleFromList     [ slideSymmetric 0.1 (PName "Small") (PWeight 5) Tune,       slideSymmetric 1.0 (PName "Medium") (PWeight 2) Tune,       slideSymmetric 5.0 (PName "Large") (PWeight 2) Tune,@@ -50,52 +48,57 @@ mon :: Monitor Double mon = Monitor monStd [] [] -nBurn :: Maybe Int-nBurn = Just 20--nAutoTune :: Maybe Int-nAutoTune = Just 10--nIter :: Int-nIter = 200- spec :: Spec-spec =-  describe "saveStatus and loadStatus" $+spec = do+  describe "save and load" $     it "doesn't change the MCMC chain" $       do-        gen <- create+        gen <- R.create         let s =-              force $-                quiet $-                  saveWith 100 $-                    status "SaveSpec" (const 1) lh proposals mon 0 nBurn nAutoTune nIter gen-        saveStatus "SaveSpec.json" s-        s' <- loadStatus (const 1) lh proposals mon Nothing "SaveSpec.json"-        r <- mh s-        r' <- mh s'-        -- Done during 'loadStatus'.-        -- removeFile "SaveSpec.json"-        item r `shouldBe` item r'+              Settings+                (AnalysisName "SaveSpec")+                (BurnInWithAutoTuning 20 10)+                (Iterations 200)+                Overwrite+                Sequential+                NoSave+                Quiet+        c <- fromMHG <$> mhg noPrior lh proposals mon 0 gen+        savedChain <- toSavedChain c+        c' <- fromSavedChain noPrior lh proposals mon savedChain+        putStrLn "@load . save@ should be @id@."+        link c `shouldBe` link c'+        iteration c `shouldBe` iteration c'+        frozenT1 <- freezeT (trace c)+        frozenT1' <- freezeT (trace c')+        frozenT1 `shouldBe` frozenT1'+        -- g1 <- R.save $ generator c+        -- g1' <- R.save $ generator c'+        -- g1 `shouldBe` g1'+        putStrLn "Sampling from the chains should be the same."+        r <- fromMHG <$> mcmc s (MHG c)+        r' <- fromMHG <$> mcmc s (MHG c')+        link r `shouldBe` link r'         iteration r `shouldBe` iteration r'-        trace r `shouldBe` trace r'-        g <- save $ generator r-        g' <- save $ generator r'-        g `shouldBe` g'+        frozenT2 <- freezeT (trace c)+        frozenT2' <- freezeT (trace c')+        frozenT2 `shouldBe` frozenT2'+        g2 <- R.save $ generator r+        g2' <- R.save $ generator r'+        g2 `shouldBe` g2' --- -- TODO: Splitmix. This will only work with a splittable generator--- -- because getNIterations changes the generator.+-- -- TODO. -- describe "mhContinue" --   $ it "mh 200 + mhContinue 200 == mh 400" --   $ do --     gen1 <- create---     let s1 = status "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 400 gen1+--     let s1 = chain "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 400 gen1 --     r1 <- mh s1 --     gen2 <- create---     let s2 = status "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 200 gen2+--     let s2 = chain "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 200 gen2 --     r2' <- mh s2 --     r2 <- mhContinue 200 r2'---     item r1 `shouldBe` item r2+--     link r1 `shouldBe` link r2 --     iteration r1 `shouldBe` iteration r2 --     trace r1 `shouldBe` trace r2 --     g <- save $ generator r1