{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE RecordWildCards #-}

-- |
-- Module      : Control.Monad.Bayes.Inference.SMC
-- Description : Sequential Monte Carlo (SMC)
-- Copyright   : (c) Adam Scibior, 2015-2020
-- License     : MIT
-- Maintainer  : leonhard.markert@tweag.io
-- Stability   : experimental
-- Portability : GHC
--
-- Sequential Monte Carlo (SMC) sampling.
--
-- Arnaud Doucet and Adam M. Johansen. 2011. A tutorial on particle filtering and smoothing: fifteen years later. In /The Oxford Handbook of Nonlinear Filtering/, Dan Crisan and Boris Rozovskii (Eds.). Oxford University Press, Chapter 8.
module Control.Monad.Bayes.Inference.SMC
  ( smc,
    smcPush,
    SMCConfig (..),
  )
where

import Control.Monad.Bayes.Class (MonadDistribution, MonadMeasure)
import Control.Monad.Bayes.Population
  ( Population,
    pushEvidence,
    withParticles,
  )
import Control.Monad.Bayes.Sequential.Coroutine as Coroutine

data SMCConfig m = SMCConfig
  { forall (m :: * -> *).
SMCConfig m -> forall x. Population m x -> Population m x
resampler :: forall x. Population m x -> Population m x,
    forall (m :: * -> *). SMCConfig m -> Int
numSteps :: Int,
    forall (m :: * -> *). SMCConfig m -> Int
numParticles :: Int
  }

-- | Sequential importance resampling.
-- Basically an SMC template that takes a custom resampler.
smc ::
  MonadDistribution m =>
  SMCConfig m ->
  Coroutine.Sequential (Population m) a ->
  Population m a
smc :: forall (m :: * -> *) a.
MonadDistribution m =>
SMCConfig m -> Sequential (Population m) a -> Population m a
smc SMCConfig {Int
forall x. Population m x -> Population m x
numParticles :: Int
numSteps :: Int
resampler :: forall x. Population m x -> Population m x
numParticles :: forall (m :: * -> *). SMCConfig m -> Int
numSteps :: forall (m :: * -> *). SMCConfig m -> Int
resampler :: forall (m :: * -> *).
SMCConfig m -> forall x. Population m x -> Population m x
..} =
  forall (m :: * -> *) a.
Monad m =>
(forall x. m x -> m x) -> Int -> Sequential m a -> m a
Coroutine.sequentially forall x. Population m x -> Population m x
resampler Int
numSteps
    forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (m :: * -> *) a.
(forall x. m x -> m x) -> Sequential m a -> Sequential m a
Coroutine.hoistFirst (forall (m :: * -> *) a.
Monad m =>
Int -> Population m a -> Population m a
withParticles Int
numParticles)

-- | Sequential Monte Carlo with multinomial resampling at each timestep.
-- Weights are normalized at each timestep and the total weight is pushed
-- as a score into the transformed monad.
smcPush ::
  MonadMeasure m => SMCConfig m -> Coroutine.Sequential (Population m) a -> Population m a
smcPush :: forall (m :: * -> *) a.
MonadMeasure m =>
SMCConfig m -> Sequential (Population m) a -> Population m a
smcPush SMCConfig m
config = forall (m :: * -> *) a.
MonadDistribution m =>
SMCConfig m -> Sequential (Population m) a -> Population m a
smc SMCConfig m
config {resampler :: forall x. Population m x -> Population m x
resampler = (forall (m :: * -> *) a.
MonadFactor m =>
Population m a -> Population m a
pushEvidence forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (m :: * -> *).
SMCConfig m -> forall x. Population m x -> Population m x
resampler SMCConfig m
config)}