{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE BlockArguments #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ViewPatterns #-}

module PPL.Sampling
  ( mh,
  )
where

import Control.Monad
import Control.Monad.IO.Class
import Data.IORef
import qualified Data.Vector.Hashtables as H
import qualified Data.Vector.Unboxed as V
import Data.Word
import Numeric.Log
import PPL.Internal
import Streaming.Prelude (Of, Stream, yield)
import System.Random (StdGen)
import qualified System.Random as R

mh :: (MonadIO m) => StdGen -> Double -> Meas a -> Stream (Of (a, Log Double)) m ()
mh g p m = do
  let (g0, g1) = R.split g
  hm <- liftIO $ H.initialize 0
  omega <- liftIO $ newIORef (hm, g0)
  let (x, w) = head $ samples m $ newTree omega
  step g1 omega x w
  where
    step !g0 !omega !x !w = do
      let (Exp . log -> r, R.split -> (g1, g2)) = R.random g0
      omega' <- mutate g1 omega
      let (!x', !w') = head $ samples m $ newTree omega'
          ratio = w' / w
          (omega'', x'', w'') =
            if r < ratio
              then (omega', x', w')
              else (omega, x, w)
      yield (x'', w'')
      step g2 omega'' x'' w''

    mutate :: (MonadIO m) => StdGen -> IORef (HashMap Word64 Double, StdGen) -> m (IORef (HashMap Word64 Double, StdGen))
    mutate g omega = liftIO $ do
      (m, g0) <- readIORef omega
      m' <- H.clone m
      ks <- H.keys m
      let (rs, qs) = splitAt (1 + floor (p * (n - 1))) (R.randoms g)
          n = fromIntegral (V.length ks)
      void $ zipWithM (\r q -> H.insert m' (ks V.! floor (r * n)) q) rs qs
      newIORef (m', g0)