diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,14 +1,17 @@ -# Markov chain Monte Carlo +# Markov chain Monte Carlo sampler

-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 1; +- Metropolis-coupled Markov chain Monte Carlo (also known as parallel + tempering) 2 , 3. + ## 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 + +1 Geyer, C. J., Introduction to Markov chain Monte Carlo, In Handbook of +Markov Chain Monte Carlo (pp. 45) (2011). CRC press. + +2 Geyer, C. J., Markov chain monte carlo maximum likelihood, Computing +Science and Statistics, Proceedings of the 23rd Symposium on the Interface, +(), (1991). + +3 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). diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -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) diff --git a/bench/Normal.hs b/bench/Normal.hs --- a/bench/Normal.hs +++ b/bench/Normal.hs @@ -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 diff --git a/bench/Poisson.hs b/bench/Poisson.hs --- a/bench/Poisson.hs +++ b/bench/Poisson.hs @@ -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 diff --git a/mcmc.cabal b/mcmc.cabal --- a/mcmc.cabal +++ b/mcmc.cabal @@ -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 -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 +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 diff --git a/src/Mcmc.hs b/src/Mcmc.hs --- a/src/Mcmc.hs +++ b/src/Mcmc.hs @@ -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 diff --git a/src/Mcmc/Algorithm.hs b/src/Mcmc/Algorithm.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Algorithm.hs @@ -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 () diff --git a/src/Mcmc/Algorithm/MC3.hs b/src/Mcmc/Algorithm/MC3.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Algorithm/MC3.hs @@ -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'} diff --git a/src/Mcmc/Algorithm/Metropolis.hs b/src/Mcmc/Algorithm/Metropolis.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Algorithm/Metropolis.hs @@ -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 diff --git a/src/Mcmc/Chain/Chain.hs b/src/Mcmc/Chain/Chain.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Chain/Chain.hs @@ -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 + } diff --git a/src/Mcmc/Chain/Link.hs b/src/Mcmc/Chain/Link.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Chain/Link.hs @@ -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 diff --git a/src/Mcmc/Chain/Save.hs b/src/Mcmc/Chain/Save.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Chain/Save.hs @@ -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 diff --git a/src/Mcmc/Chain/Trace.hs b/src/Mcmc/Chain/Trace.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Chain/Trace.hs @@ -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 diff --git a/src/Mcmc/Environment.hs b/src/Mcmc/Environment.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Environment.hs @@ -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 diff --git a/src/Mcmc/Internal/ByteString.hs b/src/Mcmc/Internal/ByteString.hs --- a/src/Mcmc/Internal/ByteString.hs +++ b/src/Mcmc/Internal/ByteString.hs @@ -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 diff --git a/src/Mcmc/Internal/Random.hs b/src/Mcmc/Internal/Random.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Internal/Random.hs @@ -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 diff --git a/src/Mcmc/Internal/Shuffle.hs b/src/Mcmc/Internal/Shuffle.hs --- a/src/Mcmc/Internal/Shuffle.hs +++ b/src/Mcmc/Internal/Shuffle.hs @@ -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 diff --git a/src/Mcmc/Item.hs b/src/Mcmc/Item.hs deleted file mode 100644 --- a/src/Mcmc/Item.hs +++ /dev/null @@ -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 diff --git a/src/Mcmc/Mcmc.hs b/src/Mcmc/Mcmc.hs --- a/src/Mcmc/Mcmc.hs +++ b/src/Mcmc/Mcmc.hs @@ -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} diff --git a/src/Mcmc/Metropolis.hs b/src/Mcmc/Metropolis.hs deleted file mode 100644 --- a/src/Mcmc/Metropolis.hs +++ /dev/null @@ -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." diff --git a/src/Mcmc/Monitor.hs b/src/Mcmc/Monitor.hs --- a/src/Mcmc/Monitor.hs +++ b/src/Mcmc/Monitor.hs @@ -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) diff --git a/src/Mcmc/Monitor/Parameter.hs b/src/Mcmc/Monitor/Parameter.hs --- a/src/Mcmc/Monitor/Parameter.hs +++ b/src/Mcmc/Monitor/Parameter.hs @@ -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 :: diff --git a/src/Mcmc/Monitor/ParameterBatch.hs b/src/Mcmc/Monitor/ParameterBatch.hs --- a/src/Mcmc/Monitor/ParameterBatch.hs +++ b/src/Mcmc/Monitor/ParameterBatch.hs @@ -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) diff --git a/src/Mcmc/Monitor/Time.hs b/src/Mcmc/Monitor/Time.hs --- a/src/Mcmc/Monitor/Time.hs +++ b/src/Mcmc/Monitor/Time.hs @@ -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." diff --git a/src/Mcmc/Prior.hs b/src/Mcmc/Prior.hs --- a/src/Mcmc/Prior.hs +++ b/src/Mcmc/Prior.hs @@ -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. -- diff --git a/src/Mcmc/Proposal.hs b/src/Mcmc/Proposal.hs --- a/src/Mcmc/Proposal.hs +++ b/src/Mcmc/Proposal.hs @@ -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 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 diff --git a/src/Mcmc/Proposal/Bactrian.hs b/src/Mcmc/Proposal/Bactrian.hs --- a/src/Mcmc/Proposal/Bactrian.hs +++ b/src/Mcmc/Proposal/Bactrian.hs @@ -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 diff --git a/src/Mcmc/Proposal/Generic.hs b/src/Mcmc/Proposal/Generic.hs --- a/src/Mcmc/Proposal/Generic.hs +++ b/src/Mcmc/Proposal/Generic.hs @@ -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 -- diff --git a/src/Mcmc/Proposal/Scale.hs b/src/Mcmc/Proposal/Scale.hs --- a/src/Mcmc/Proposal/Scale.hs +++ b/src/Mcmc/Proposal/Scale.hs @@ -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 diff --git a/src/Mcmc/Proposal/Simplex.hs b/src/Mcmc/Proposal/Simplex.hs --- a/src/Mcmc/Proposal/Simplex.hs +++ b/src/Mcmc/Proposal/Simplex.hs @@ -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 diff --git a/src/Mcmc/Proposal/Slide.hs b/src/Mcmc/Proposal/Slide.hs --- a/src/Mcmc/Proposal/Slide.hs +++ b/src/Mcmc/Proposal/Slide.hs @@ -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 diff --git a/src/Mcmc/Save.hs b/src/Mcmc/Save.hs deleted file mode 100644 --- a/src/Mcmc/Save.hs +++ /dev/null @@ -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 diff --git a/src/Mcmc/Settings.hs b/src/Mcmc/Settings.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Settings.hs @@ -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 diff --git a/src/Mcmc/Status.hs b/src/Mcmc/Status.hs deleted file mode 100644 --- a/src/Mcmc/Status.hs +++ /dev/null @@ -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} diff --git a/src/Mcmc/Trace.hs b/src/Mcmc/Trace.hs deleted file mode 100644 --- a/src/Mcmc/Trace.hs +++ /dev/null @@ -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 diff --git a/src/Mcmc/Verbosity.hs b/src/Mcmc/Verbosity.hs deleted file mode 100644 --- a/src/Mcmc/Verbosity.hs +++ /dev/null @@ -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) diff --git a/test/Mcmc/ProposalSpec.hs b/test/Mcmc/ProposalSpec.hs --- a/test/Mcmc/ProposalSpec.hs +++ b/test/Mcmc/ProposalSpec.hs @@ -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 diff --git a/test/Mcmc/SaveSpec.hs b/test/Mcmc/SaveSpec.hs --- a/test/Mcmc/SaveSpec.hs +++ b/test/Mcmc/SaveSpec.hs @@ -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