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

-- |
-- Module      :  Mcmc.MarginalLikelihood
-- Description :  Calculate the marginal likelihood
-- Copyright   :  (c) Dominik Schrempf 2021
-- License     :  GPL-3.0-or-later
--
-- Maintainer  :  dominik.schrempf@gmail.com
-- Stability   :  unstable
-- Portability :  portable
--
-- Creation date: Mon Jan 11 16:34:18 2021.
module Mcmc.MarginalLikelihood
  ( MarginalLikelihood,
    NPoints (..),
    MLAlgorithm (..),
    MLSettings (..),
    marginalLikelihood,
  )
where

import Control.Monad
import Control.Monad.IO.Class
import qualified Control.Monad.Parallel as P
import Control.Monad.Trans.Reader
import Data.Aeson
import Data.List hiding (cycle)
import qualified Data.Map.Strict as M
import qualified Data.Vector as VB
import qualified Data.Vector.Unboxed as VU
import Mcmc.Acceptance
import Mcmc.Algorithm.MHG
import Mcmc.Chain.Chain
import Mcmc.Chain.Link
import Mcmc.Chain.Trace
import Mcmc.Cycle
import Mcmc.Environment
import Mcmc.Internal.Random
import Mcmc.Likelihood
import Mcmc.Logger
import Mcmc.Mcmc
import Mcmc.Monitor
import Mcmc.Prior
import Mcmc.Settings
import Numeric.Log hiding (sum)
import System.Directory
import System.Random.MWC
import Text.Printf
import Text.Show.Pretty
import Prelude hiding (cycle)

-- | Marginal likelihood values are stored in log domain.
type MarginalLikelihood = Log Double

-- Reciprocal temperature value traversed along the path integral.
type Point = Double

-- | The number of points used to approximate the path integral.
newtype NPoints = NPoints {NPoints -> Int
fromNPoints :: Int}
  deriving (NPoints -> NPoints -> Bool
(NPoints -> NPoints -> Bool)
-> (NPoints -> NPoints -> Bool) -> Eq NPoints
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: NPoints -> NPoints -> Bool
$c/= :: NPoints -> NPoints -> Bool
== :: NPoints -> NPoints -> Bool
$c== :: NPoints -> NPoints -> Bool
Eq, ReadPrec [NPoints]
ReadPrec NPoints
Int -> ReadS NPoints
ReadS [NPoints]
(Int -> ReadS NPoints)
-> ReadS [NPoints]
-> ReadPrec NPoints
-> ReadPrec [NPoints]
-> Read NPoints
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [NPoints]
$creadListPrec :: ReadPrec [NPoints]
readPrec :: ReadPrec NPoints
$creadPrec :: ReadPrec NPoints
readList :: ReadS [NPoints]
$creadList :: ReadS [NPoints]
readsPrec :: Int -> ReadS NPoints
$creadsPrec :: Int -> ReadS NPoints
Read, Int -> NPoints -> ShowS
[NPoints] -> ShowS
NPoints -> String
(Int -> NPoints -> ShowS)
-> (NPoints -> String) -> ([NPoints] -> ShowS) -> Show NPoints
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [NPoints] -> ShowS
$cshowList :: [NPoints] -> ShowS
show :: NPoints -> String
$cshow :: NPoints -> String
showsPrec :: Int -> NPoints -> ShowS
$cshowsPrec :: Int -> NPoints -> ShowS
Show)

-- | Algorithms to calculate the marginal likelihood.
data MLAlgorithm
  = -- | Use a classical path integral. Also known as thermodynamic integration.
    -- In particular, /Annealing-Melting Integration/ is used.
    --
    -- See Lartillot, N., & Philippe, H., Computing Bayes Factors Using
    -- Thermodynamic Integration, Systematic Biology, 55(2), 195–207 (2006).
    -- http://dx.doi.org/10.1080/10635150500433722
    ThermodynamicIntegration
  | -- | Use stepping stone sampling.
    --
    -- See Xie, W., Lewis, P. O., Fan, Y., Kuo, L., & Chen, M., Improving
    -- marginal likelihood estimation for Bayesian phylogenetic model selection,
    -- Systematic Biology, 60(2), 150–160 (2010).
    -- http://dx.doi.org/10.1093/sysbio/syq085
    --
    -- Or Fan, Y., Wu, R., Chen, M., Kuo, L., & Lewis, P. O., Choosing among
    -- partition models in bayesian phylogenetics, Molecular Biology and
    -- Evolution, 28(1), 523–532 (2010). http://dx.doi.org/10.1093/molbev/msq224
    SteppingStoneSampling
  deriving (MLAlgorithm -> MLAlgorithm -> Bool
(MLAlgorithm -> MLAlgorithm -> Bool)
-> (MLAlgorithm -> MLAlgorithm -> Bool) -> Eq MLAlgorithm
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: MLAlgorithm -> MLAlgorithm -> Bool
$c/= :: MLAlgorithm -> MLAlgorithm -> Bool
== :: MLAlgorithm -> MLAlgorithm -> Bool
$c== :: MLAlgorithm -> MLAlgorithm -> Bool
Eq, ReadPrec [MLAlgorithm]
ReadPrec MLAlgorithm
Int -> ReadS MLAlgorithm
ReadS [MLAlgorithm]
(Int -> ReadS MLAlgorithm)
-> ReadS [MLAlgorithm]
-> ReadPrec MLAlgorithm
-> ReadPrec [MLAlgorithm]
-> Read MLAlgorithm
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [MLAlgorithm]
$creadListPrec :: ReadPrec [MLAlgorithm]
readPrec :: ReadPrec MLAlgorithm
$creadPrec :: ReadPrec MLAlgorithm
readList :: ReadS [MLAlgorithm]
$creadList :: ReadS [MLAlgorithm]
readsPrec :: Int -> ReadS MLAlgorithm
$creadsPrec :: Int -> ReadS MLAlgorithm
Read, Int -> MLAlgorithm -> ShowS
[MLAlgorithm] -> ShowS
MLAlgorithm -> String
(Int -> MLAlgorithm -> ShowS)
-> (MLAlgorithm -> String)
-> ([MLAlgorithm] -> ShowS)
-> Show MLAlgorithm
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [MLAlgorithm] -> ShowS
$cshowList :: [MLAlgorithm] -> ShowS
show :: MLAlgorithm -> String
$cshow :: MLAlgorithm -> String
showsPrec :: Int -> MLAlgorithm -> ShowS
$cshowsPrec :: Int -> MLAlgorithm -> ShowS
Show)

-- | Settings of the marginal likelihood estimation.
data MLSettings = MLSettings
  { MLSettings -> AnalysisName
mlAnalysisName :: AnalysisName,
    MLSettings -> MLAlgorithm
mlAlgorithm :: MLAlgorithm,
    MLSettings -> NPoints
mlNPoints :: NPoints,
    -- | Initial burn in at the starting point of the path.
    MLSettings -> BurnInSettings
mlInitialBurnIn :: BurnInSettings,
    -- | Repetitive burn in at each point on the path.
    MLSettings -> BurnInSettings
mlPointBurnIn :: BurnInSettings,
    -- | The number of iterations performed at each point.
    MLSettings -> Iterations
mlIterations :: Iterations,
    MLSettings -> ExecutionMode
mlExecutionMode :: ExecutionMode,
    MLSettings -> LogMode
mlLogMode :: LogMode,
    MLSettings -> Verbosity
mlVerbosity :: Verbosity
  }
  deriving (MLSettings -> MLSettings -> Bool
(MLSettings -> MLSettings -> Bool)
-> (MLSettings -> MLSettings -> Bool) -> Eq MLSettings
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: MLSettings -> MLSettings -> Bool
$c/= :: MLSettings -> MLSettings -> Bool
== :: MLSettings -> MLSettings -> Bool
$c== :: MLSettings -> MLSettings -> Bool
Eq, ReadPrec [MLSettings]
ReadPrec MLSettings
Int -> ReadS MLSettings
ReadS [MLSettings]
(Int -> ReadS MLSettings)
-> ReadS [MLSettings]
-> ReadPrec MLSettings
-> ReadPrec [MLSettings]
-> Read MLSettings
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [MLSettings]
$creadListPrec :: ReadPrec [MLSettings]
readPrec :: ReadPrec MLSettings
$creadPrec :: ReadPrec MLSettings
readList :: ReadS [MLSettings]
$creadList :: ReadS [MLSettings]
readsPrec :: Int -> ReadS MLSettings
$creadsPrec :: Int -> ReadS MLSettings
Read, Int -> MLSettings -> ShowS
[MLSettings] -> ShowS
MLSettings -> String
(Int -> MLSettings -> ShowS)
-> (MLSettings -> String)
-> ([MLSettings] -> ShowS)
-> Show MLSettings
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [MLSettings] -> ShowS
$cshowList :: [MLSettings] -> ShowS
show :: MLSettings -> String
$cshow :: MLSettings -> String
showsPrec :: Int -> MLSettings -> ShowS
$cshowsPrec :: Int -> MLSettings -> ShowS
Show)

instance HasAnalysisName MLSettings where
  getAnalysisName :: MLSettings -> AnalysisName
getAnalysisName = MLSettings -> AnalysisName
mlAnalysisName

instance HasExecutionMode MLSettings where
  getExecutionMode :: MLSettings -> ExecutionMode
getExecutionMode = MLSettings -> ExecutionMode
mlExecutionMode

instance HasLogMode MLSettings where
  getLogMode :: MLSettings -> LogMode
getLogMode = MLSettings -> LogMode
mlLogMode

instance HasVerbosity MLSettings where
  getVerbosity :: MLSettings -> Verbosity
getVerbosity = MLSettings -> Verbosity
mlVerbosity

type ML a = ReaderT (Environment MLSettings) IO a

-- See 'getPoints'. Alpha=0.3 is the standard choice.
alpha :: Double
alpha :: Double
alpha = Double
0.3

-- Distribute the points according to a skewed beta distribution with given
-- 'alpha' value. If alpha is below 1.0, more points at lower values, which is
-- desired. It is inconvenient that the reciprocal temperatures are denoted as
-- beta, and we also use the beta distribution :). Don't mix them up!
--
-- See discussion in Xie, W., Lewis, P. O., Fan, Y., Kuo, L., & Chen, M.,
-- Improving marginal likelihood estimation for bayesian phylogenetic model
-- selection, Systematic Biology, 60(2), 150–160 (2010).
-- http://dx.doi.org/10.1093/sysbio/syq085
--
-- Or Figure 1 in Höhna, S., Landis, M. J., & Huelsenbeck, J. P., Parallel power
-- posterior analyses for fast computation of marginal likelihoods in
-- phylogenetics (2017). http://dx.doi.org/10.1101/104422
getPoints :: NPoints -> [Point]
getPoints :: NPoints -> [Double]
getPoints NPoints
x = [Int -> Double
forall a a. (Fractional a, Integral a) => a -> a
f Int
i Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1.0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
alpha) | Int
i <- [Int
0 .. Int
k1]]
  where
    k :: Int
k = NPoints -> Int
fromNPoints NPoints
x
    k1 :: Int
k1 = Int -> Int
forall a. Enum a => a -> a
pred Int
k
    f :: a -> a
f a
j = a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
j a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k1

sampleAtPoint ::
  ToJSON a =>
  Point ->
  Settings ->
  LikelihoodFunction a ->
  MHG a ->
  ML (MHG a)
sampleAtPoint :: Double -> Settings -> LikelihoodFunction a -> MHG a -> ML (MHG a)
sampleAtPoint Double
x Settings
ss LikelihoodFunction a
lhf MHG a
a = do
  MHG a
a'' <- IO (MHG a) -> ML (MHG a)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (MHG a) -> ML (MHG a)) -> IO (MHG a) -> ML (MHG a)
forall a b. (a -> b) -> a -> b
$ Settings -> MHG a -> IO (MHG a)
forall a. Algorithm a => Settings -> a -> IO a
mcmc Settings
ss' MHG a
a'
  let ch'' :: Chain a
ch'' = MHG a -> Chain a
forall a. MHG a -> Chain a
fromMHG MHG a
a''
      ac :: Acceptance (Proposal a)
ac = Chain a -> Acceptance (Proposal a)
forall a. Chain a -> Acceptance (Proposal a)
acceptance Chain a
ch''
      mAr :: Maybe (Map (Proposal a) Double)
mAr = Map (Proposal a) (Maybe Double) -> Maybe (Map (Proposal a) Double)
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence (Map (Proposal a) (Maybe Double)
 -> Maybe (Map (Proposal a) Double))
-> Map (Proposal a) (Maybe Double)
-> Maybe (Map (Proposal a) Double)
forall a b. (a -> b) -> a -> b
$ Acceptance (Proposal a) -> Map (Proposal a) (Maybe Double)
forall k. Acceptance k -> Map k (Maybe Double)
acceptanceRates Acceptance (Proposal a)
ac
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logDebugB ByteString
"sampleAtPoint: Summarize cycle."
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logDebugB (ByteString -> Logger (Environment MLSettings) ())
-> ByteString -> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ Acceptance (Proposal a) -> Cycle a -> ByteString
forall a. Acceptance (Proposal a) -> Cycle a -> ByteString
summarizeCycle Acceptance (Proposal a)
ac (Cycle a -> ByteString) -> Cycle a -> ByteString
forall a b. (a -> b) -> a -> b
$ Chain a -> Cycle a
forall a. Chain a -> Cycle a
cycle Chain a
ch''
  case Maybe (Map (Proposal a) Double)
mAr of
    Maybe (Map (Proposal a) Double)
Nothing -> ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logWarnB ByteString
"Some acceptance rates are unavailable. The tuning period may be too small."
    Just Map (Proposal a) Double
ar -> do
      Bool
-> Logger (Environment MLSettings) ()
-> Logger (Environment MLSettings) ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
unless (Map (Proposal a) Double -> Bool
forall k a. Map k a -> Bool
M.null (Map (Proposal a) Double -> Bool)
-> Map (Proposal a) Double -> Bool
forall a b. (a -> b) -> a -> b
$ (Double -> Bool)
-> Map (Proposal a) Double -> Map (Proposal a) Double
forall a k. (a -> Bool) -> Map k a -> Map k a
M.filter (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0.1) Map (Proposal a) Double
ar) (Logger (Environment MLSettings) ()
 -> Logger (Environment MLSettings) ())
-> Logger (Environment MLSettings) ()
-> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logWarnB ByteString
"Some acceptance rates are below 0.1."
      Bool
-> Logger (Environment MLSettings) ()
-> Logger (Environment MLSettings) ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
unless (Map (Proposal a) Double -> Bool
forall k a. Map k a -> Bool
M.null (Map (Proposal a) Double -> Bool)
-> Map (Proposal a) Double -> Bool
forall a b. (a -> b) -> a -> b
$ (Double -> Bool)
-> Map (Proposal a) Double -> Map (Proposal a) Double
forall a k. (a -> Bool) -> Map k a -> Map k a
M.filter (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0.9) Map (Proposal a) Double
ar) (Logger (Environment MLSettings) ()
 -> Logger (Environment MLSettings) ())
-> Logger (Environment MLSettings) ()
-> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logWarnB ByteString
"Some acceptance rates are above 0.9."
  MHG a -> ML (MHG a)
forall (m :: * -> *) a. Monad m => a -> m a
return MHG a
a''
  where
    -- For debugging set a proper analysis name.
    nm :: AnalysisName
nm = Settings -> AnalysisName
sAnalysisName Settings
ss
    getName :: Point -> AnalysisName
    getName :: Double -> AnalysisName
getName Double
y = AnalysisName
nm AnalysisName -> AnalysisName -> AnalysisName
forall a. Semigroup a => a -> a -> a
<> String -> AnalysisName
AnalysisName (String
"/" String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String -> Double -> String
forall r. PrintfType r => String -> r
printf String
"point%.8f" Double
y)
    ss' :: Settings
ss' = Settings
ss {sAnalysisName :: AnalysisName
sAnalysisName = Double -> AnalysisName
getName Double
x}
    -- Amend the likelihood function. Don't calculate the likelihood when the
    -- point is 0.0.
    lhf' :: LikelihoodFunction a
lhf' = if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0 then Log Double -> LikelihoodFunction a
forall a b. a -> b -> a
const Log Double
1.0 else (Log Double -> Log Double -> Log Double
forall a. Floating a => a -> a -> a
** Double -> Log Double
forall a. a -> Log a
Exp (Double -> Double
forall a. Floating a => a -> a
log Double
x)) (Log Double -> Log Double)
-> LikelihoodFunction a -> LikelihoodFunction a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. LikelihoodFunction a
lhf
    -- Amend the MHG algorithm.
    ch :: Chain a
ch = MHG a -> Chain a
forall a. MHG a -> Chain a
fromMHG MHG a
a
    l :: Link a
l = Chain a -> Link a
forall a. Chain a -> Link a
link Chain a
ch
    ch' :: Chain a
ch' =
      Chain a
ch
        { -- Important: Update the likelihood using the new likelihood function.
          link :: Link a
link = Link a
l {likelihood :: Log Double
likelihood = LikelihoodFunction a
lhf' LikelihoodFunction a -> LikelihoodFunction a
forall a b. (a -> b) -> a -> b
$ Link a -> a
forall a. Link a -> a
state Link a
l},
          iteration :: Int
iteration = Int
0,
          start :: Int
start = Int
0,
          likelihoodFunction :: LikelihoodFunction a
likelihoodFunction = LikelihoodFunction a
lhf'
        }
    a' :: MHG a
a' = Chain a -> MHG a
forall a. Chain a -> MHG a
MHG Chain a
ch'

traversePoints ::
  ToJSON a =>
  -- Current point.
  Int ->
  NPoints ->
  [Point] ->
  Settings ->
  LikelihoodFunction a ->
  MHG a ->
  -- For each point a vector of obtained likelihoods stored in the log domain.
  ML [VU.Vector Likelihood]
traversePoints :: Int
-> NPoints
-> [Double]
-> Settings
-> LikelihoodFunction a
-> MHG a
-> ML [Vector (Log Double)]
traversePoints Int
_ NPoints
_ [] Settings
_ LikelihoodFunction a
_ MHG a
_ = [Vector (Log Double)] -> ML [Vector (Log Double)]
forall (m :: * -> *) a. Monad m => a -> m a
return []
traversePoints Int
i NPoints
k (Double
b : [Double]
bs) Settings
ss LikelihoodFunction a
lhf MHG a
a = do
  String -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS (String -> Logger (Environment MLSettings) ())
-> String -> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ String
"Point " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Int -> String
forall a. Show a => a -> String
show Int
i String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
" of " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Int -> String
forall a. Show a => a -> String
show Int
k' String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
": " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Double -> String
forall a. Show a => a -> String
show Double
b String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
"."
  MHG a
a' <- Double -> Settings -> LikelihoodFunction a -> MHG a -> ML (MHG a)
forall a.
ToJSON a =>
Double -> Settings -> LikelihoodFunction a -> MHG a -> ML (MHG a)
sampleAtPoint Double
b Settings
ss LikelihoodFunction a
lhf MHG a
a
  -- Get the links samples at this point.
  Vector (Link a)
ls <- IO (Vector (Link a))
-> ReaderT (Environment MLSettings) IO (Vector (Link a))
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (Vector (Link a))
 -> ReaderT (Environment MLSettings) IO (Vector (Link a)))
-> IO (Vector (Link a))
-> ReaderT (Environment MLSettings) IO (Vector (Link a))
forall a b. (a -> b) -> a -> b
$ Int -> Trace a -> IO (Vector (Link a))
forall a. Int -> Trace a -> IO (Vector (Link a))
takeT Int
n (Trace a -> IO (Vector (Link a)))
-> Trace a -> IO (Vector (Link a))
forall a b. (a -> b) -> a -> b
$ Chain a -> Trace a
forall a. Chain a -> Trace a
trace (Chain a -> Trace a) -> Chain a -> Trace a
forall a b. (a -> b) -> a -> b
$ MHG a -> Chain a
forall a. MHG a -> Chain a
fromMHG MHG a
a'
  -- Extract the likelihoods.
  --
  -- NOTE: This could be sped up by mapping (** -b) on the power likelihoods.
  --
  -- NOTE: This bang is an important one, because if the lhs are not strictly
  -- calculated here, the complete MCMC runs are dragged along before doing so
  -- resulting in a severe memory leak.
  let !lhs :: Vector (Log Double)
lhs = Vector (Log Double) -> Vector (Log Double)
forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
VU.convert (Vector (Log Double) -> Vector (Log Double))
-> Vector (Log Double) -> Vector (Log Double)
forall a b. (a -> b) -> a -> b
$ (Link a -> Log Double) -> Vector (Link a) -> Vector (Log Double)
forall a b. (a -> b) -> Vector a -> Vector b
VB.map (LikelihoodFunction a
lhf LikelihoodFunction a -> (Link a -> a) -> Link a -> Log Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Link a -> a
forall a. Link a -> a
state) Vector (Link a)
ls
  -- Sample the other points.
  [Vector (Log Double)]
lhss <- Int
-> NPoints
-> [Double]
-> Settings
-> LikelihoodFunction a
-> MHG a
-> ML [Vector (Log Double)]
forall a.
ToJSON a =>
Int
-> NPoints
-> [Double]
-> Settings
-> LikelihoodFunction a
-> MHG a
-> ML [Vector (Log Double)]
traversePoints (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) NPoints
k [Double]
bs Settings
ss LikelihoodFunction a
lhf MHG a
a'
  [Vector (Log Double)] -> ML [Vector (Log Double)]
forall (m :: * -> *) a. Monad m => a -> m a
return ([Vector (Log Double)] -> ML [Vector (Log Double)])
-> [Vector (Log Double)] -> ML [Vector (Log Double)]
forall a b. (a -> b) -> a -> b
$ Vector (Log Double)
lhs Vector (Log Double)
-> [Vector (Log Double)] -> [Vector (Log Double)]
forall a. a -> [a] -> [a]
: [Vector (Log Double)]
lhss
  where
    n :: Int
n = Iterations -> Int
fromIterations (Iterations -> Int) -> Iterations -> Int
forall a b. (a -> b) -> a -> b
$ Settings -> Iterations
sIterations Settings
ss
    (NPoints Int
k') = NPoints
k

mlRun ::
  ToJSON a =>
  NPoints ->
  [Point] ->
  ExecutionMode ->
  Verbosity ->
  PriorFunction a ->
  LikelihoodFunction a ->
  Cycle a ->
  Monitor a ->
  a ->
  GenIO ->
  -- For each point a vector of likelihoods stored in log domain.
  ML [VU.Vector Likelihood]
mlRun :: NPoints
-> [Double]
-> ExecutionMode
-> Verbosity
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML [Vector (Log Double)]
mlRun NPoints
k [Double]
xs ExecutionMode
em Verbosity
vb PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 GenIO
g = do
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logDebugB ByteString
"mlRun: Begin."
  MLSettings
s <- (Environment MLSettings -> MLSettings)
-> ReaderT (Environment MLSettings) IO MLSettings
forall (m :: * -> *) r a. Monad m => (r -> a) -> ReaderT r m a
reader Environment MLSettings -> MLSettings
forall s. Environment s -> s
settings
  let nm :: AnalysisName
nm = MLSettings -> AnalysisName
mlAnalysisName MLSettings
s
      is :: Iterations
is = MLSettings -> Iterations
mlIterations MLSettings
s
      biI :: BurnInSettings
biI = MLSettings -> BurnInSettings
mlInitialBurnIn MLSettings
s
      biP :: BurnInSettings
biP = MLSettings -> BurnInSettings
mlPointBurnIn MLSettings
s
      -- Only log sub MCMC samplers when debugging.
      vb' :: Verbosity
vb' = case Verbosity
vb of
        Verbosity
Debug -> Verbosity
Debug
        Verbosity
_ -> Verbosity
Quiet
      trLen :: TraceLength
trLen = Int -> TraceLength
TraceMinimum (Int -> TraceLength) -> Int -> TraceLength
forall a b. (a -> b) -> a -> b
$ Iterations -> Int
fromIterations Iterations
is
      ssI :: Settings
ssI = AnalysisName
-> BurnInSettings
-> Iterations
-> TraceLength
-> ExecutionMode
-> ParallelizationMode
-> SaveMode
-> LogMode
-> Verbosity
-> Settings
Settings AnalysisName
nm BurnInSettings
biI (Int -> Iterations
Iterations Int
0) TraceLength
trLen ExecutionMode
em ParallelizationMode
Sequential SaveMode
NoSave LogMode
LogFileOnly Verbosity
vb'
      ssP :: Settings
ssP = AnalysisName
-> BurnInSettings
-> Iterations
-> TraceLength
-> ExecutionMode
-> ParallelizationMode
-> SaveMode
-> LogMode
-> Verbosity
-> Settings
Settings AnalysisName
nm BurnInSettings
biP Iterations
is TraceLength
trLen ExecutionMode
em ParallelizationMode
Sequential SaveMode
NoSave LogMode
LogFileOnly Verbosity
vb'
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logDebugB ByteString
"mlRun: Initialize MHG algorithm."
  MHG a
a0 <- IO (MHG a) -> ReaderT (Environment MLSettings) IO (MHG a)
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (MHG a) -> ReaderT (Environment MLSettings) IO (MHG a))
-> IO (MHG a) -> ReaderT (Environment MLSettings) IO (MHG a)
forall a b. (a -> b) -> a -> b
$ Settings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> IO (MHG a)
forall a.
Settings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> IO (MHG a)
mhg Settings
ssI PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 GenIO
g
  String -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logDebugS (String -> Logger (Environment MLSettings) ())
-> String -> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ String
"mlRun: Perform initial burn in at first point " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Double -> String
forall a. Show a => a -> String
show Double
x0 String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
"."
  MHG a
a1 <- Double
-> Settings
-> PriorFunction a
-> MHG a
-> ReaderT (Environment MLSettings) IO (MHG a)
forall a.
ToJSON a =>
Double -> Settings -> LikelihoodFunction a -> MHG a -> ML (MHG a)
sampleAtPoint Double
x0 Settings
ssI PriorFunction a
lhf MHG a
a0
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logDebugB ByteString
"mlRun: Traverse points."
  Int
-> NPoints
-> [Double]
-> Settings
-> PriorFunction a
-> MHG a
-> ML [Vector (Log Double)]
forall a.
ToJSON a =>
Int
-> NPoints
-> [Double]
-> Settings
-> LikelihoodFunction a
-> MHG a
-> ML [Vector (Log Double)]
traversePoints Int
1 NPoints
k [Double]
xs Settings
ssP PriorFunction a
lhf MHG a
a1
  where
    x0 :: Double
x0 = [Double] -> Double
forall a. [a] -> a
head [Double]
xs

-- Use lists since the number of points is expected to be low.
integrateSimpsonTriangle ::
  -- X values.
  [Point] ->
  -- Y values.
  [Double] ->
  -- Integral.
  Double
integrateSimpsonTriangle :: [Double] -> [Double] -> Double
integrateSimpsonTriangle [Double]
xs [Double]
ys = Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
forall p. Num p => [p] -> [p] -> p
go [Double]
xs [Double]
ys
  where
    go :: [p] -> [p] -> p
go (p
p0 : p
p1 : [p]
ps) (p
z0 : p
z1 : [p]
zs) = (p
z0 p -> p -> p
forall a. Num a => a -> a -> a
+ p
z1) p -> p -> p
forall a. Num a => a -> a -> a
* (p
p1 p -> p -> p
forall a. Num a => a -> a -> a
- p
p0) p -> p -> p
forall a. Num a => a -> a -> a
+ [p] -> [p] -> p
go (p
p1 p -> [p] -> [p]
forall a. a -> [a] -> [a]
: [p]
ps) (p
z1 p -> [p] -> [p]
forall a. a -> [a] -> [a]
: [p]
zs)
    go [p]
_ [p]
_ = p
0

tiWrapper ::
  ToJSON a =>
  MLSettings ->
  PriorFunction a ->
  LikelihoodFunction a ->
  Cycle a ->
  Monitor a ->
  a ->
  GenIO ->
  ML MarginalLikelihood
tiWrapper :: MLSettings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML (Log Double)
tiWrapper MLSettings
s PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 GenIO
g = do
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB ByteString
"Path integral (thermodynamic integration)."
  [Gen RealWorld
g0, Gen RealWorld
g1] <- Int
-> Gen (PrimState (ReaderT (Environment MLSettings) IO))
-> ReaderT
     (Environment MLSettings)
     IO
     [Gen (PrimState (ReaderT (Environment MLSettings) IO))]
forall (m :: * -> *).
PrimMonad m =>
Int -> Gen (PrimState m) -> m [Gen (PrimState m)]
splitGen Int
2 GenIO
Gen (PrimState (ReaderT (Environment MLSettings) IO))
g

  -- Parallel execution of both path integrals.
  [[Vector (Log Double)]
lhssForward, [Vector (Log Double)]
lhssBackward] <-
    [ML [Vector (Log Double)]]
-> ReaderT (Environment MLSettings) IO [[Vector (Log Double)]]
forall (m :: * -> *) a. MonadParallel m => [m a] -> m [a]
P.sequence
      [ NPoints
-> [Double]
-> ExecutionMode
-> Verbosity
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML [Vector (Log Double)]
forall a.
ToJSON a =>
NPoints
-> [Double]
-> ExecutionMode
-> Verbosity
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML [Vector (Log Double)]
mlRun NPoints
k [Double]
bsForward ExecutionMode
em Verbosity
vb PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 Gen RealWorld
GenIO
g0,
        NPoints
-> [Double]
-> ExecutionMode
-> Verbosity
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML [Vector (Log Double)]
forall a.
ToJSON a =>
NPoints
-> [Double]
-> ExecutionMode
-> Verbosity
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML [Vector (Log Double)]
mlRun NPoints
k [Double]
bsBackward ExecutionMode
em Verbosity
vb PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 Gen RealWorld
GenIO
g1
      ]
  Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasStartingTime e, HasVerbosity e) =>
Logger e ()
logInfoEndTime

  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logDebugB ByteString
"tiWrapper: Calculate mean log likelihoods."
  -- It is important to average across the log likelihoods here (and not the
  -- likelihoods). I am not exactly sure why this is.
  let getMeanLogLhs :: [Vector (Log Double)] -> [Double]
getMeanLogLhs = (Vector (Log Double) -> Double)
-> [Vector (Log Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Vector (Log Double)
x -> Vector Double -> Double
forall a. (Unbox a, Num a) => Vector a -> a
VU.sum ((Log Double -> Double) -> Vector (Log Double) -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map Log Double -> Double
forall a. Log a -> a
ln Vector (Log Double)
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Vector (Log Double) -> Int
forall a. Unbox a => Vector a -> Int
VU.length Vector (Log Double)
x))
      mlForward :: Double
mlForward = [Double] -> [Double] -> Double
integrateSimpsonTriangle [Double]
bsForward ([Vector (Log Double)] -> [Double]
getMeanLogLhs [Vector (Log Double)]
lhssForward)
      mlBackward :: Double
mlBackward = Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ [Double] -> [Double] -> Double
integrateSimpsonTriangle [Double]
bsBackward ([Vector (Log Double)] -> [Double]
getMeanLogLhs [Vector (Log Double)]
lhssBackward)
  String -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logDebugS (String -> Logger (Environment MLSettings) ())
-> String -> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ String
"tiWrapper: Marginal log likelihood of forward integral: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show Double
mlForward
  String -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logDebugS (String -> Logger (Environment MLSettings) ())
-> String -> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ String
"tiWrapper: Marginal log likelihood of backward integral: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show Double
mlBackward
  let mean :: Double
mean = Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
mlForward Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
mlBackward)
  String -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logDebugS (String -> Logger (Environment MLSettings) ())
-> String -> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ String
"tiWrapper: The mean is: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show Double
mean
  Log Double -> ML (Log Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Log Double -> ML (Log Double)) -> Log Double -> ML (Log Double)
forall a b. (a -> b) -> a -> b
$ Double -> Log Double
forall a. a -> Log a
Exp Double
mean
  where
    k :: NPoints
k = MLSettings -> NPoints
mlNPoints MLSettings
s
    bsForward :: [Double]
bsForward = NPoints -> [Double]
getPoints NPoints
k
    bsBackward :: [Double]
bsBackward = [Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
bsForward
    em :: ExecutionMode
em = MLSettings -> ExecutionMode
mlExecutionMode MLSettings
s
    vb :: Verbosity
vb = MLSettings -> Verbosity
mlVerbosity MLSettings
s

-- Helper function to exponentiate log domain values with a double value.
pow' :: Log Double -> Double -> Log Double
pow' :: Log Double -> Double -> Log Double
pow' Log Double
x Double
p = Double -> Log Double
forall a. a -> Log a
Exp (Double -> Log Double) -> Double -> Log Double
forall a b. (a -> b) -> a -> b
$ Log Double -> Double
forall a. Log a -> a
ln Log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p

-- See Xie2010 p. 153, bottom left.
sssCalculateMarginalLikelihood :: [Point] -> [VU.Vector Likelihood] -> MarginalLikelihood
sssCalculateMarginalLikelihood :: [Double] -> [Vector (Log Double)] -> Log Double
sssCalculateMarginalLikelihood [Double]
xs [Vector (Log Double)]
lhss = [Log Double] -> Log Double
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ([Log Double] -> Log Double) -> [Log Double] -> Log Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Vector (Log Double) -> Log Double)
-> [Double] -> [Double] -> [Vector (Log Double)] -> [Log Double]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 Double -> Double -> Vector (Log Double) -> Log Double
f [Double]
xs ([Double] -> [Double]
forall a. [a] -> [a]
tail [Double]
xs) [Vector (Log Double)]
lhss
  where
    f :: Point -> Point -> VU.Vector Likelihood -> MarginalLikelihood
    -- f beta_{k-1} beta_k lhs_{k-1}
    f :: Double -> Double -> Vector (Log Double) -> Log Double
f Double
bkm1 Double
bk Vector (Log Double)
lhs = Log Double
n1 Log Double -> Log Double -> Log Double
forall a. Num a => a -> a -> a
* Vector (Log Double) -> Log Double
forall a. (Unbox a, Num a) => Vector a -> a
VU.sum Vector (Log Double)
lhsPowered
      where
        n1 :: Log Double
n1 = Log Double -> Log Double
forall a. Fractional a => a -> a
recip (Log Double -> Log Double) -> Log Double -> Log Double
forall a b. (a -> b) -> a -> b
$ Int -> Log Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Log Double) -> Int -> Log Double
forall a b. (a -> b) -> a -> b
$ Vector (Log Double) -> Int
forall a. Unbox a => Vector a -> Int
VU.length Vector (Log Double)
lhs
        dbeta :: Double
dbeta = Double
bk Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bkm1
        lhsPowered :: Vector (Log Double)
lhsPowered = (Log Double -> Log Double)
-> Vector (Log Double) -> Vector (Log Double)
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map (Log Double -> Double -> Log Double
`pow'` Double
dbeta) Vector (Log Double)
lhs

-- -- Numerical stability by factoring out lhMax. But no observed
-- -- improvement towards the standard version.
--
-- f bkm1 bk lhs = n1 * pow' lhMax dbeta * VU.sum lhsNormedPowered
--   where n1 = recip $ fromIntegral $ VU.length lhs
--         lhMax = VU.maximum lhs
--         dbeta = bk - bkm1
--         lhsNormed = VU.map (/lhMax) lhs
--         lhsNormedPowered = VU.map (`pow'` dbeta) lhsNormed

-- -- Computation of the log of the marginal likelihood. According to the paper,
-- -- this estimator is biased and I did not observe any improvements compared
-- -- to the direct estimator implemented above.
--
-- -- See Xie2010 p. 153, top right.
-- sssCalculateMarginalLikelihood' :: [Point] -> [VU.Vector Likelihood] -> MarginalLikelihood
-- sssCalculateMarginalLikelihood' xs lhss = Exp $ sum $ zipWith3 f xs (tail xs) lhss
--   where f :: Point -> Point -> VU.Vector Likelihood -> Double
--         -- f beta_{k-1} beta_k lhs_{k-1}
--         f bkm1 bk lhs = dbeta * llhMax + log (n1 * VU.sum lhsNormedPowered)
--           where dbeta = bk - bkm1
--                 llhMax = ln $ VU.maximum lhs
--                 n1 = recip $ fromIntegral $ VU.length lhs
--                 llhs = VU.map ln lhs
--                 llhsNormed = VU.map (\x -> x - llhMax) llhs
--                 lhsNormedPowered = VU.map (\x -> exp $ dbeta * x) llhsNormed
sssWrapper ::
  ToJSON a =>
  MLSettings ->
  PriorFunction a ->
  LikelihoodFunction a ->
  Cycle a ->
  Monitor a ->
  a ->
  GenIO ->
  ML MarginalLikelihood
sssWrapper :: MLSettings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML (Log Double)
sssWrapper MLSettings
s PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 GenIO
g = do
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB ByteString
"Stepping stone sampling."
  [Vector (Log Double)]
logLhss <- NPoints
-> [Double]
-> ExecutionMode
-> Verbosity
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML [Vector (Log Double)]
forall a.
ToJSON a =>
NPoints
-> [Double]
-> ExecutionMode
-> Verbosity
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML [Vector (Log Double)]
mlRun NPoints
k [Double]
bsForward' ExecutionMode
em Verbosity
vb PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 GenIO
g
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB ByteString
"The last point does not need to be sampled with stepping stone sampling."
  ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logDebugB ByteString
"sssWrapper: Calculate marginal likelihood."
  Log Double -> ML (Log Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Log Double -> ML (Log Double)) -> Log Double -> ML (Log Double)
forall a b. (a -> b) -> a -> b
$ [Double] -> [Vector (Log Double)] -> Log Double
sssCalculateMarginalLikelihood [Double]
bsForward [Vector (Log Double)]
logLhss
  where
    k :: NPoints
k = MLSettings -> NPoints
mlNPoints MLSettings
s
    bsForward :: [Double]
bsForward = NPoints -> [Double]
getPoints NPoints
k
    bsForward' :: [Double]
bsForward' = [Double] -> [Double]
forall a. [a] -> [a]
init [Double]
bsForward
    em :: ExecutionMode
em = MLSettings -> ExecutionMode
mlExecutionMode MLSettings
s
    vb :: Verbosity
vb = MLSettings -> Verbosity
mlVerbosity MLSettings
s

-- | Estimate the marginal likelihood.
marginalLikelihood ::
  ToJSON a =>
  MLSettings ->
  PriorFunction a ->
  LikelihoodFunction a ->
  Cycle a ->
  Monitor a ->
  InitialState a ->
  -- | A source of randomness. For reproducible runs, make sure to use
  -- generators with the same seed.
  GenIO ->
  IO MarginalLikelihood
marginalLikelihood :: MLSettings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> IO (Log Double)
marginalLikelihood MLSettings
s PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 GenIO
g = do
  -- Initialize.
  Environment MLSettings
e <- MLSettings -> IO (Environment MLSettings)
forall s.
(HasAnalysisName s, HasExecutionMode s, HasLogMode s,
 HasVerbosity s) =>
s -> IO (Environment s)
initializeEnvironment MLSettings
s

  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (MLSettings -> Verbosity
mlVerbosity MLSettings
s Verbosity -> Verbosity -> Bool
forall a. Eq a => a -> a -> Bool
== Verbosity
Debug) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
    let n :: String
n = AnalysisName -> String
fromAnalysisName (AnalysisName -> String) -> AnalysisName -> String
forall a b. (a -> b) -> a -> b
$ MLSettings -> AnalysisName
mlAnalysisName MLSettings
s
    Bool -> String -> IO ()
createDirectoryIfMissing Bool
True String
n

  -- Run.
  ML (Log Double) -> Environment MLSettings -> IO (Log Double)
forall r (m :: * -> *) a. ReaderT r m a -> r -> m a
runReaderT
    ( do
        Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasStartingTime e, HasVerbosity e) =>
Logger e ()
logInfoStartingTime
        ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB ByteString
"Estimate marginal likelihood."
        ByteString -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logDebugB ByteString
"marginalLikelihood: The marginal likelihood settings are:"
        String -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logDebugS (String -> Logger (Environment MLSettings) ())
-> String -> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ MLSettings -> String
forall a. Show a => a -> String
ppShow MLSettings
s
        Log Double
val <- case MLSettings -> MLAlgorithm
mlAlgorithm MLSettings
s of
          MLAlgorithm
ThermodynamicIntegration -> MLSettings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML (Log Double)
forall a.
ToJSON a =>
MLSettings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML (Log Double)
tiWrapper MLSettings
s PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 GenIO
g
          MLAlgorithm
SteppingStoneSampling -> MLSettings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML (Log Double)
forall a.
ToJSON a =>
MLSettings
-> PriorFunction a
-> PriorFunction a
-> Cycle a
-> Monitor a
-> a
-> GenIO
-> ML (Log Double)
sssWrapper MLSettings
s PriorFunction a
prf PriorFunction a
lhf Cycle a
cc Monitor a
mn a
i0 GenIO
g
        String -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS (String -> Logger (Environment MLSettings) ())
-> String -> Logger (Environment MLSettings) ()
forall a b. (a -> b) -> a -> b
$ String
"Marginal log likelihood: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show (Log Double -> Double
forall a. Log a -> a
ln Log Double
val)
        -- TODO: Simulation variance.
        String -> Logger (Environment MLSettings) ()
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
"The simulation variance is not yet available."
        Log Double -> ML (Log Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Log Double
val
    )
    Environment MLSettings
e