{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE OverloadedStrings #-}

-- |
-- Module      :  Mcmc.Algorithm.MHG
-- Description :  Metropolis-Hastings-Green algorithm
-- Copyright   :  2021 Dominik Schrempf
-- 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.MHG
  ( MHG (..),
    mhg,
    mhgSave,
    mhgLoad,
    mhgLoadUnsafe,
    MHGRatio,
    mhgAccept,
  )
where

import Codec.Compression.GZip
import Control.Monad
import Control.Monad.IO.Class
import Control.Parallel.Strategies
import Data.Aeson
import qualified Data.ByteString.Lazy.Char8 as BL
import Data.Maybe
import Data.Time
import qualified Data.Vector as VB
import Mcmc.Acceptance
import Mcmc.Algorithm
import Mcmc.Chain.Chain
import Mcmc.Chain.Link
import Mcmc.Chain.Save
import Mcmc.Chain.Trace
import Mcmc.Cycle
import Mcmc.Likelihood
import Mcmc.Monitor
import Mcmc.Posterior
import Mcmc.Prior hiding (uniform)
import Mcmc.Proposal
import Mcmc.Settings
import Numeric.Log
import System.Random.Stateful
import Prelude hiding (cycle)

-- | The MHG algorithm.
newtype MHG a = MHG {forall a. MHG a -> Chain a
fromMHG :: Chain a}

instance ToJSON a => Algorithm (MHG a) where
  aName :: MHG a -> [Char]
aName = forall a b. a -> b -> a
const [Char]
"Metropolis-Hastings-Green (MHG)"
  aIteration :: MHG a -> Int
aIteration = forall a. Chain a -> Int
iteration forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. MHG a -> Chain a
fromMHG
  aIsInvalidState :: MHG a -> Bool
aIsInvalidState = forall a. MHG a -> Bool
mhgIsInvalidState
  aIterate :: IterationMode -> ParallelizationMode -> MHG a -> IO (MHG a)
aIterate = forall a.
IterationMode -> ParallelizationMode -> MHG a -> IO (MHG a)
mhgIterate
  aAutoTune :: TuningType -> Int -> MHG a -> IO (MHG a)
aAutoTune = forall a. TuningType -> Int -> MHG a -> IO (MHG a)
mhgAutoTune
  aResetAcceptance :: ResetAcceptance -> MHG a -> MHG a
aResetAcceptance = forall a. ResetAcceptance -> MHG a -> MHG a
mhgResetAcceptance
  aCleanAfterBurnIn :: TraceLength -> MHG a -> IO (MHG a)
aCleanAfterBurnIn = forall a. TraceLength -> MHG a -> IO (MHG a)
mhgCleanAfterBurnIn
  aSummarizeCycle :: IterationMode -> MHG a -> ByteString
aSummarizeCycle = forall a. IterationMode -> MHG a -> ByteString
mhgSummarizeCycle
  aOpenMonitors :: AnalysisName -> ExecutionMode -> MHG a -> IO (MHG a)
aOpenMonitors = forall a. AnalysisName -> ExecutionMode -> MHG a -> IO (MHG a)
mhgOpenMonitors
  aExecuteMonitors :: Verbosity -> UTCTime -> Int -> MHG a -> IO (Maybe ByteString)
aExecuteMonitors = forall a.
Verbosity -> UTCTime -> Int -> MHG a -> IO (Maybe ByteString)
mhgExecuteMonitors
  aStdMonitorHeader :: MHG a -> ByteString
aStdMonitorHeader = forall a. MHG a -> ByteString
mhgStdMonitorHeader
  aCloseMonitors :: MHG a -> IO (MHG a)
aCloseMonitors = forall a. MHG a -> IO (MHG a)
mhgCloseMonitors
  aSave :: AnalysisName -> MHG a -> IO ()
aSave = forall a. ToJSON a => AnalysisName -> MHG a -> IO ()
mhgSave

-- Calculate required length of trace. The length may be larger during burn in,
-- because the tuners of some proposals (e.g., HMC, NUTS) require the states of
-- the last tuning interval.
getTraceLength ::
  Maybe BurnInSettings ->
  TraceLength ->
  Monitor a ->
  Cycle a ->
  Int
getTraceLength :: forall a.
Maybe BurnInSettings -> TraceLength -> Monitor a -> Cycle a -> Int
getTraceLength Maybe BurnInSettings
burnIn TraceLength
tl Monitor a
mn Cycle a
cc = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum forall a b. (a -> b) -> a -> b
$ Int
minimumTraceLength forall a. a -> [a] -> [a]
: Int
bi forall a. a -> [a] -> [a]
: [Int]
batchMonitorSizes
  where
    batchMonitorSizes :: [Int]
batchMonitorSizes = forall a b. (a -> b) -> [a] -> [b]
map forall a. MonitorBatch a -> Int
getMonitorBatchSize forall a b. (a -> b) -> a -> b
$ forall a. Monitor a -> [MonitorBatch a]
mBatches Monitor a
mn
    minimumTraceLength :: Int
minimumTraceLength = case TraceLength
tl of
      TraceLength
TraceAuto -> Int
1
      TraceMinimum Int
n -> Int
n
    bi :: Int
bi = case (forall a. Cycle a -> Bool
ccRequireTrace Cycle a
cc, Maybe BurnInSettings
burnIn) of
      (Bool
True, Just (BurnInWithAutoTuning Int
_ Int
n)) -> Int
n
      (Bool
True, Just (BurnInWithCustomAutoTuning [Int]
ns [Int]
ms)) -> forall a. Ord a => a -> a -> a
max (forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum forall a b. (a -> b) -> a -> b
$ Int
0 forall a. a -> [a] -> [a]
: [Int]
ns) (forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum forall a b. (a -> b) -> a -> b
$ Int
0 forall a. a -> [a] -> [a]
: [Int]
ms)
      (Bool, Maybe BurnInSettings)
_ -> Int
0

-- | Initialize an MHG algorithm.
--
-- NOTE: Computation in the 'IO' Monad is necessary because the trace is
-- mutable.
mhg ::
  Settings ->
  PriorFunction a ->
  LikelihoodFunction a ->
  Cycle a ->
  Monitor a ->
  InitialState a ->
  StdGen ->
  IO (MHG a)
mhg :: forall a.
Settings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> StdGen
-> IO (MHG a)
mhg Settings
s PriorFunction a
pr PriorFunction a
lh Cycle a
cc Monitor a
mn a
i0 StdGen
g = do
  -- The trace is a mutable vector and the mutable state needs to be handled by
  -- a monad.
  Trace a
tr <- forall a. Int -> Link a -> IO (Trace a)
replicateT Int
tl Link a
l0
  IOGenM StdGen
gm <- forall (m :: * -> *) g. MonadIO m => g -> m (IOGenM g)
newIOGenM StdGen
g
  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG forall a b. (a -> b) -> a -> b
$ forall a.
Link a
-> Int
-> Trace a
-> Acceptances (Proposal a)
-> IOGenM StdGen
-> Int
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> Chain a
Chain Link a
l0 Int
0 Trace a
tr Acceptances (Proposal a)
ac IOGenM StdGen
gm Int
0 PriorFunction a
pr PriorFunction a
lh Cycle a
cc Monitor a
mn
  where
    l0 :: Link a
l0 = forall a. a -> KernelRatio -> KernelRatio -> Link a
Link a
i0 (PriorFunction a
pr a
i0) (PriorFunction a
lh a
i0)
    ac :: Acceptances (Proposal a)
ac = forall k. Ord k => [k] -> Acceptances k
emptyA forall a b. (a -> b) -> a -> b
$ forall a. Cycle a -> [Proposal a]
ccProposals Cycle a
cc
    tl :: Int
tl = forall a.
Maybe BurnInSettings -> TraceLength -> Monitor a -> Cycle a -> Int
getTraceLength (forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ Settings -> BurnInSettings
sBurnIn Settings
s) (Settings -> TraceLength
sTraceLength Settings
s) Monitor a
mn Cycle a
cc

mhgFn :: AnalysisName -> FilePath
mhgFn :: AnalysisName -> [Char]
mhgFn (AnalysisName [Char]
nm) = [Char]
nm forall a. [a] -> [a] -> [a]
++ [Char]
".mcmc.mhg"

-- | Save an MHG algorithm.
mhgSave ::
  ToJSON a =>
  AnalysisName ->
  MHG a ->
  IO ()
mhgSave :: forall a. ToJSON a => AnalysisName -> MHG a -> IO ()
mhgSave AnalysisName
nm (MHG Chain a
c) = do
  SavedChain a
savedChain <- forall a. Chain a -> IO (SavedChain a)
toSavedChain Chain a
c
  [Char] -> ByteString -> IO ()
BL.writeFile (AnalysisName -> [Char]
mhgFn AnalysisName
nm) forall a b. (a -> b) -> a -> b
$ ByteString -> ByteString
compress forall a b. (a -> b) -> a -> b
$ forall a. ToJSON a => a -> ByteString
encode SavedChain a
savedChain

-- | Load an MHG algorithm.
--
-- Also create a backup of the save.
--
-- See 'Mcmc.Mcmc.mcmcContinue'.
mhgLoad ::
  FromJSON a =>
  PriorFunction a ->
  LikelihoodFunction a ->
  Cycle a ->
  Monitor a ->
  AnalysisName ->
  IO (MHG a)
mhgLoad :: forall a.
FromJSON a =>
PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> AnalysisName
-> IO (MHG a)
mhgLoad = forall a.
FromJSON a =>
(PriorFunction a
 -> PriorFunction a
 -> Cycle a
 -> Monitor a
 -> SavedChain a
 -> IO (Chain a))
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> AnalysisName
-> IO (MHG a)
mhgLoadWith forall a.
PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> SavedChain a
-> IO (Chain a)
fromSavedChain

-- | Like 'mhgLoad' but do not perform sanity checks.
--
-- Also create a backup of the save.
--
-- Useful when restarting a run with changed prior function, likelihood function
-- or proposals. Use with care!
mhgLoadUnsafe ::
  FromJSON a =>
  PriorFunction a ->
  LikelihoodFunction a ->
  Cycle a ->
  Monitor a ->
  AnalysisName ->
  IO (MHG a)
mhgLoadUnsafe :: forall a.
FromJSON a =>
PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> AnalysisName
-> IO (MHG a)
mhgLoadUnsafe = forall a.
FromJSON a =>
(PriorFunction a
 -> PriorFunction a
 -> Cycle a
 -> Monitor a
 -> SavedChain a
 -> IO (Chain a))
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> AnalysisName
-> IO (MHG a)
mhgLoadWith forall a.
PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> SavedChain a
-> IO (Chain a)
fromSavedChainUnsafe

-- Nice type :-).
mhgLoadWith ::
  FromJSON a =>
  (PriorFunction a -> LikelihoodFunction a -> Cycle a -> Monitor a -> SavedChain a -> IO (Chain a)) ->
  PriorFunction a ->
  LikelihoodFunction a ->
  Cycle a ->
  Monitor a ->
  AnalysisName ->
  IO (MHG a)
mhgLoadWith :: forall a.
FromJSON a =>
(PriorFunction a
 -> PriorFunction a
 -> Cycle a
 -> Monitor a
 -> SavedChain a
 -> IO (Chain a))
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> AnalysisName
-> IO (MHG a)
mhgLoadWith PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> SavedChain a
-> IO (Chain a)
f PriorFunction a
pr PriorFunction a
lh Cycle a
cc Monitor a
mn AnalysisName
nm = do
  Either [Char] (SavedChain a)
savedChain <- forall a. FromJSON a => ByteString -> Either [Char] a
eitherDecode forall b c a. (b -> c) -> (a -> b) -> a -> c
. ByteString -> ByteString
decompress forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Char] -> IO ByteString
BL.readFile [Char]
fn
  Chain a
chain <- forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either forall a. HasCallStack => [Char] -> a
error (PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> SavedChain a
-> IO (Chain a)
f PriorFunction a
pr PriorFunction a
lh Cycle a
cc Monitor a
mn) Either [Char] (SavedChain a)
savedChain
  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG Chain a
chain
  where
    -- fnBak = mhgFn $ AnalysisName $ (fromAnalysisName nm ++ ".bak")
    fn :: [Char]
fn = AnalysisName -> [Char]
mhgFn AnalysisName
nm

-- | MHG ratios are stored in log domain.
type MHGRatio = Log Double

-- The MHG ratio. This implementation has the following properties:
--
-- - The kernel ratio and the Jacobian are checked carefully and should be
-- - strictly positive, finite numbers.
--
-- - The ratio is 'Infinity' if fX is zero. In this case, the proposal is always
--   accepted.
--
-- - The ratio is 'NaN' if fY and fX are zero. In this case, the proposal is
--   always rejected.
--
-- This means that a chain in a state with posterior probability zero (fX=0) can
-- only move if a state with non-zero posterior probability is proposed.
-- Otherwise it is stuck. Therefore, I print a warning when the posterior
-- probability is zero in the beginning of the MCMC run. This is probably not
-- the best behavior, but see below.
--
-- 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 [1],
-- [2]), or (b) almost surely reject the proposal when either fY or q are zero
-- (Chapter 1 of [1]).
--
-- Since I trust the author of Chapter 1 (Charles Geyer) I choose to follow
-- option (b). However, Option (a) is more user-friendly.
--
-- [1] Handbook of Markov chain Monte Carlo (2011), CRC press.
--
-- [2] Dellaportas, P., & Roberts, G. O., An introduction to MCMC, Lecture Notes
-- in Statistics, (), 1–41 (2003).
-- http://dx.doi.org/10.1007/978-0-387-21811-3_1.
mhgRatio :: Posterior -> Posterior -> KernelRatio -> Jacobian -> MHGRatio
mhgRatio :: KernelRatio
-> KernelRatio -> KernelRatio -> KernelRatio -> KernelRatio
mhgRatio KernelRatio
fX KernelRatio
fY KernelRatio
q KernelRatio
j
  | KernelRatio
q forall a. Eq a => a -> a -> Bool
== KernelRatio
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"mhgRatio: Kernel ratio is negative infinity. Use 'ForceReject'."
  | KernelRatio
q forall a. Eq a => a -> a -> Bool
== KernelRatio
1.0 forall a. Fractional a => a -> a -> a
/ KernelRatio
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"mhgRatio: Kernel ratio is infinity. Use 'ForceAccept'."
  | KernelRatio
q forall a. Eq a => a -> a -> Bool
== KernelRatio
0.0 forall a. Fractional a => a -> a -> a
/ KernelRatio
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"mhgRatio: Kernel ratio is NaN."
  | KernelRatio
j forall a. Eq a => a -> a -> Bool
== KernelRatio
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"mhgRatio: Jacobian is negative infinity. Use 'ForceReject'."
  | KernelRatio
j forall a. Eq a => a -> a -> Bool
== KernelRatio
1.0 forall a. Fractional a => a -> a -> a
/ KernelRatio
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"mhgRatio: Jacobian is infinity. Use 'ForceAccept'."
  | KernelRatio
j forall a. Eq a => a -> a -> Bool
== KernelRatio
0.0 forall a. Fractional a => a -> a -> a
/ KernelRatio
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"mhgRatio: Jacobian is NaN."
  | Bool
otherwise = KernelRatio
fY forall a. Fractional a => a -> a -> a
/ KernelRatio
fX forall a. Num a => a -> a -> a
* KernelRatio
q forall a. Num a => a -> a -> a
* KernelRatio
j
{-# INLINE mhgRatio #-}

-- | Accept or reject a proposal with given MHG ratio?
mhgAccept :: MHGRatio -> IOGenM StdGen -> IO Bool
mhgAccept :: KernelRatio -> IOGenM StdGen -> IO Bool
mhgAccept KernelRatio
r IOGenM StdGen
g
  | forall a. Log a -> a
ln KernelRatio
r forall a. Ord a => a -> a -> Bool
>= Double
0.0 = forall (m :: * -> *) a. Monad m => a -> m a
return Bool
True
  | Bool
otherwise = do
      Double
b <- forall a g (m :: * -> *).
(UniformRange a, StatefulGen g m) =>
(a, a) -> g -> m a
uniformRM (Double
0, Double
1) IOGenM StdGen
g
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Double
b forall a. Ord a => a -> a -> Bool
< forall a. Floating a => a -> a
exp (forall a. Log a -> a
ln KernelRatio
r)

mhgPropose :: MHG a -> Proposal a -> IO (MHG a)
mhgPropose :: forall a. MHG a -> Proposal a -> IO (MHG a)
mhgPropose (MHG Chain a
c) Proposal a
p = do
  -- 1. Sample new state.
  (!PResult a
pres, !Maybe AcceptanceRates
mcs) <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ PFunction a
s a
x IOGenM StdGen
g
  -- 2. Define new prior and likelihood calculation functions. Avoid actual
  -- calculation of the values.
  --
  -- Most often, parallelization is not helpful, because the prior and
  -- likelihood functions are too fast; see
  -- https://stackoverflow.com/a/46603680/3536806.
  let calcPrLh :: a -> (KernelRatio, KernelRatio)
calcPrLh a
y = (PriorFunction a
pF a
y, PriorFunction a
lF a
y) forall a. a -> Strategy a -> a
`using` forall a b. Strategy a -> Strategy b -> Strategy (a, b)
parTuple2 forall a. NFData a => Strategy a
rdeepseq forall a. NFData a => Strategy a
rdeepseq
      accept :: a -> KernelRatio -> KernelRatio -> f (MHG a)
accept a
y KernelRatio
pr KernelRatio
lh =
        let !ac' :: Acceptances (Proposal a)
ac' = forall k.
Ord k =>
Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k
pushAccept Maybe AcceptanceRates
mcs Proposal a
p Acceptances (Proposal a)
ac
         in forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG forall a b. (a -> b) -> a -> b
$ Chain a
c {link :: Link a
link = forall a. a -> KernelRatio -> KernelRatio -> Link a
Link a
y KernelRatio
pr KernelRatio
lh, acceptances :: Acceptances (Proposal a)
acceptances = Acceptances (Proposal a)
ac'}
      reject :: IO (MHG a)
reject =
        let !ac' :: Acceptances (Proposal a)
ac' = forall k.
Ord k =>
Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k
pushReject Maybe AcceptanceRates
mcs Proposal a
p Acceptances (Proposal a)
ac
         in forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG forall a b. (a -> b) -> a -> b
$ Chain a
c {acceptances :: Acceptances (Proposal a)
acceptances = Acceptances (Proposal a)
ac'}
  -- 3. Accept or reject.
  --
  -- 3a. When rejection is inevitable, avoid calculation of the prior, the
  -- likelihood and the MHG ratio.
  case PResult a
pres of
    PResult a
ForceReject -> IO (MHG a)
reject
    ForceAccept a
y -> let (KernelRatio
pY, KernelRatio
lY) = a -> (KernelRatio, KernelRatio)
calcPrLh a
y in forall {f :: * -> *}.
Applicative f =>
a -> KernelRatio -> KernelRatio -> f (MHG a)
accept a
y KernelRatio
pY KernelRatio
lY
    (Propose a
y KernelRatio
q KernelRatio
j) ->
      if KernelRatio
q forall a. Ord a => a -> a -> Bool
<= KernelRatio
0.0 Bool -> Bool -> Bool
|| KernelRatio
j forall a. Ord a => a -> a -> Bool
<= KernelRatio
0.0
        then IO (MHG a)
reject
        else do
          -- 3b. Calculate Metropolis-Hastings-Green ratio.
          let (KernelRatio
pY, KernelRatio
lY) = a -> (KernelRatio, KernelRatio)
calcPrLh a
y
              !r :: KernelRatio
r = KernelRatio
-> KernelRatio -> KernelRatio -> KernelRatio -> KernelRatio
mhgRatio (KernelRatio
pX forall a. Num a => a -> a -> a
* KernelRatio
lX) (KernelRatio
pY forall a. Num a => a -> a -> a
* KernelRatio
lY) KernelRatio
q KernelRatio
j
          Bool
isAccept <- KernelRatio -> IOGenM StdGen -> IO Bool
mhgAccept KernelRatio
r IOGenM StdGen
g
          if Bool
isAccept
            then forall {f :: * -> *}.
Applicative f =>
a -> KernelRatio -> KernelRatio -> f (MHG a)
accept a
y KernelRatio
pY KernelRatio
lY
            else IO (MHG a)
reject
  where
    s :: PFunction a
s = forall a. Proposal a -> PFunction a
prFunction Proposal a
p
    (Link a
x KernelRatio
pX KernelRatio
lX) = forall a. Chain a -> Link a
link Chain a
c
    pF :: PriorFunction a
pF = forall a. Chain a -> PriorFunction a
priorFunction Chain a
c
    lF :: PriorFunction a
lF = forall a. Chain a -> PriorFunction a
likelihoodFunction Chain a
c
    ac :: Acceptances (Proposal a)
ac = forall a. Chain a -> Acceptances (Proposal a)
acceptances Chain a
c
    g :: IOGenM StdGen
g = forall a. Chain a -> IOGenM StdGen
generator Chain a
c

mhgPush :: MHG a -> IO (MHG a)
mhgPush :: forall a. MHG a -> IO (MHG a)
mhgPush (MHG Chain a
c) = do
  Trace a
t' <- forall a. Link a -> Trace a -> IO (Trace a)
pushT Link a
i Trace a
t
  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG Chain a
c {trace :: Trace a
trace = Trace a
t', iteration :: Int
iteration = forall a. Enum a => a -> a
succ Int
n}
  where
    i :: Link a
i = forall a. Chain a -> Link a
link Chain a
c
    t :: Trace a
t = forall a. Chain a -> Trace a
trace Chain a
c
    n :: Int
n = forall a. Chain a -> Int
iteration Chain a
c

-- Check if the current state is invalid.
--
-- At the moment this just checks whether the prior, likelihood, or posterior
-- are NaN or infinite.
mhgIsInvalidState :: MHG a -> Bool
mhgIsInvalidState :: forall a. MHG a -> Bool
mhgIsInvalidState MHG a
a = forall {a}. RealFloat a => Log a -> Bool
checkSoft KernelRatio
p Bool -> Bool -> Bool
|| forall {a}. RealFloat a => Log a -> Bool
check KernelRatio
l Bool -> Bool -> Bool
|| forall {a}. RealFloat a => Log a -> Bool
check (KernelRatio
p forall a. Num a => a -> a -> a
* KernelRatio
l)
  where
    x :: Link a
x = forall a. Chain a -> Link a
link forall a b. (a -> b) -> a -> b
$ forall a. MHG a -> Chain a
fromMHG MHG a
a
    p :: KernelRatio
p = forall a. Link a -> KernelRatio
prior Link a
x
    l :: KernelRatio
l = forall a. Link a -> KernelRatio
likelihood Link a
x
    check :: Log a -> Bool
check Log a
v = let v' :: a
v' = forall a. Log a -> a
ln Log a
v in forall a. RealFloat a => a -> Bool
isNaN a
v' Bool -> Bool -> Bool
|| forall a. RealFloat a => a -> Bool
isInfinite a
v' Bool -> Bool -> Bool
|| a
v' forall a. Eq a => a -> a -> Bool
== a
0
    checkSoft :: Log a -> Bool
checkSoft Log a
v = let v' :: a
v' = forall a. Log a -> a
ln Log a
v in forall a. RealFloat a => a -> Bool
isNaN a
v' Bool -> Bool -> Bool
|| forall a. RealFloat a => a -> Bool
isInfinite a
v'

-- 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 :: IterationMode -> ParallelizationMode -> MHG a -> IO (MHG a)
mhgIterate :: forall a.
IterationMode -> ParallelizationMode -> MHG a -> IO (MHG a)
mhgIterate IterationMode
m ParallelizationMode
_ MHG a
a = do
  [Proposal a]
ps <- forall g (m :: * -> *) a.
StatefulGen g m =>
IterationMode -> Cycle a -> g -> m [Proposal a]
prepareProposals IterationMode
m Cycle a
cc IOGenM StdGen
g
  MHG a
a' <- forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM forall a. MHG a -> Proposal a -> IO (MHG a)
mhgPropose MHG a
a [Proposal a]
ps
  forall a. MHG a -> IO (MHG a)
mhgPush MHG a
a'
  where
    c :: Chain a
c = forall a. MHG a -> Chain a
fromMHG MHG a
a
    cc :: Cycle a
cc = forall a. Chain a -> Cycle a
cycle Chain a
c
    g :: IOGenM StdGen
g = forall a. Chain a -> IOGenM StdGen
generator Chain a
c

mhgAutoTune :: TuningType -> Int -> MHG a -> IO (MHG a)
mhgAutoTune :: forall a. TuningType -> Int -> MHG a -> IO (MHG a)
mhgAutoTune TuningType
tt Int
n (MHG Chain a
c)
  | Bool
isIntermediate =
      forall (f :: * -> *) a. Applicative f => a -> f a
pure forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Chain a -> MHG a
MHG forall a b. (a -> b) -> a -> b
$
        if forall a. Cycle a -> Bool
ccHasIntermediateTuners Cycle a
cc
          then -- Do not provide trace when tuning intermediately.
            Chain a
c {cycle :: Cycle a
cycle = forall a.
TuningType
-> Acceptances (Proposal a)
-> Maybe (Vector a)
-> Cycle a
-> Cycle a
autoTuneCycle TuningType
tt Acceptances (Proposal a)
ac forall a. Maybe a
Nothing Cycle a
cc}
          else -- Skip intermediate tuning completely when unnecessary.
            Chain a
c
  | Bool
otherwise = do
      Maybe (Vector a)
mxs <-
        -- Provide the trace if required.
        if forall a. Cycle a -> Bool
ccRequireTrace Cycle a
cc
          then forall a. a -> Maybe a
Just forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> Vector a -> Vector b
VB.map forall a. Link a -> a
state forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall a. Int -> Trace a -> IO (Vector (Link a))
takeT Int
n Trace a
tr
          else forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a. Maybe a
Nothing
      forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG Chain a
c {cycle :: Cycle a
cycle = forall a.
TuningType
-> Acceptances (Proposal a)
-> Maybe (Vector a)
-> Cycle a
-> Cycle a
autoTuneCycle TuningType
tt Acceptances (Proposal a)
ac Maybe (Vector a)
mxs Cycle a
cc}
  where
    isIntermediate :: Bool
isIntermediate = TuningType
tt forall a. Eq a => a -> a -> Bool
== TuningType
IntermediateTuningFastProposalsOnly Bool -> Bool -> Bool
|| TuningType
tt forall a. Eq a => a -> a -> Bool
== TuningType
IntermediateTuningAllProposals
    ac :: Acceptances (Proposal a)
ac = forall a. Chain a -> Acceptances (Proposal a)
acceptances Chain a
c
    cc :: Cycle a
cc = forall a. Chain a -> Cycle a
cycle Chain a
c
    tr :: Trace a
tr = forall a. Chain a -> Trace a
trace Chain a
c

mhgResetAcceptance :: ResetAcceptance -> MHG a -> MHG a
mhgResetAcceptance :: forall a. ResetAcceptance -> MHG a -> MHG a
mhgResetAcceptance ResetAcceptance
a (MHG Chain a
c) = forall a. Chain a -> MHG a
MHG forall a b. (a -> b) -> a -> b
$ Chain a
c {acceptances :: Acceptances (Proposal a)
acceptances = forall k.
Ord k =>
ResetAcceptance -> Acceptances k -> Acceptances k
resetA ResetAcceptance
a Acceptances (Proposal a)
ac}
  where
    ac :: Acceptances (Proposal a)
ac = forall a. Chain a -> Acceptances (Proposal a)
acceptances Chain a
c

mhgCleanAfterBurnIn :: TraceLength -> MHG a -> IO (MHG a)
mhgCleanAfterBurnIn :: forall a. TraceLength -> MHG a -> IO (MHG a)
mhgCleanAfterBurnIn TraceLength
tl (MHG Chain a
c) = do
  Vector (Link a)
xs <- forall a. Int -> Trace a -> IO (Vector (Link a))
takeT Int
l Trace a
tr
  Trace a
tr' <- forall a. Vector (Link a) -> IO (Trace a)
fromVectorT Vector (Link a)
xs
  let c' :: Chain a
c' = Chain a
c {trace :: Trace a
trace = Trace a
tr'}
  forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG Chain a
c'
  where
    mn :: Monitor a
mn = forall a. Chain a -> Monitor a
monitor Chain a
c
    cc :: Cycle a
cc = forall a. Chain a -> Cycle a
cycle Chain a
c
    tr :: Trace a
tr = forall a. Chain a -> Trace a
trace Chain a
c
    l :: Int
l = forall a.
Maybe BurnInSettings -> TraceLength -> Monitor a -> Cycle a -> Int
getTraceLength forall a. Maybe a
Nothing TraceLength
tl Monitor a
mn Cycle a
cc

mhgSummarizeCycle :: IterationMode -> MHG a -> BL.ByteString
mhgSummarizeCycle :: forall a. IterationMode -> MHG a -> ByteString
mhgSummarizeCycle IterationMode
m (MHG Chain a
c) = forall a.
IterationMode -> Acceptances (Proposal a) -> Cycle a -> ByteString
summarizeCycle IterationMode
m Acceptances (Proposal a)
ac Cycle a
cc
  where
    cc :: Cycle a
cc = forall a. Chain a -> Cycle a
cycle Chain a
c
    ac :: Acceptances (Proposal a)
ac = forall a. Chain a -> Acceptances (Proposal a)
acceptances Chain a
c

mhgOpenMonitors ::
  AnalysisName ->
  ExecutionMode ->
  MHG a ->
  IO (MHG a)
mhgOpenMonitors :: forall a. AnalysisName -> ExecutionMode -> MHG a -> IO (MHG a)
mhgOpenMonitors AnalysisName
nm ExecutionMode
em (MHG Chain a
c) = do
  Monitor a
m' <- forall a.
[Char] -> [Char] -> ExecutionMode -> Monitor a -> IO (Monitor a)
mOpen [Char]
pre [Char]
"" ExecutionMode
em Monitor a
m
  forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG Chain a
c {monitor :: Monitor a
monitor = Monitor a
m'}
  where
    m :: Monitor a
m = forall a. Chain a -> Monitor a
monitor Chain a
c
    pre :: [Char]
pre = AnalysisName -> [Char]
fromAnalysisName AnalysisName
nm

mhgExecuteMonitors ::
  Verbosity ->
  -- Starting time.
  UTCTime ->
  -- Total number of iterations.
  Int ->
  MHG a ->
  IO (Maybe BL.ByteString)
mhgExecuteMonitors :: forall a.
Verbosity -> UTCTime -> Int -> MHG a -> IO (Maybe ByteString)
mhgExecuteMonitors Verbosity
vb UTCTime
t0 Int
iTotal (MHG Chain a
c) = forall a.
Verbosity
-> Int
-> Int
-> UTCTime
-> Trace a
-> Int
-> Monitor a
-> IO (Maybe ByteString)
mExec Verbosity
vb Int
i Int
i0 UTCTime
t0 Trace a
tr Int
iTotal Monitor a
m
  where
    i :: Int
i = forall a. Chain a -> Int
iteration Chain a
c
    i0 :: Int
i0 = forall a. Chain a -> Int
start Chain a
c
    tr :: Trace a
tr = forall a. Chain a -> Trace a
trace Chain a
c
    m :: Monitor a
m = forall a. Chain a -> Monitor a
monitor Chain a
c

mhgStdMonitorHeader :: MHG a -> BL.ByteString
mhgStdMonitorHeader :: forall a. MHG a -> ByteString
mhgStdMonitorHeader (MHG Chain a
c) = forall a. MonitorStdOut a -> ByteString
msHeader (forall a. Monitor a -> MonitorStdOut a
mStdOut forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> Monitor a
monitor Chain a
c)

mhgCloseMonitors :: MHG a -> IO (MHG a)
mhgCloseMonitors :: forall a. MHG a -> IO (MHG a)
mhgCloseMonitors (MHG Chain a
c) = do
  Monitor a
m' <- forall a. Monitor a -> IO (Monitor a)
mClose Monitor a
m
  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Chain a -> MHG a
MHG forall a b. (a -> b) -> a -> b
$ Chain a
c {monitor :: Monitor a
monitor = Monitor a
m'}
  where
    m :: Monitor a
m = forall a. Chain a -> Monitor a
monitor Chain a
c