1 files changed, 42 insertions, 14 deletions
diff --git a/src/PPL/Sampling.hs b/src/PPL/Sampling.hs
index 8e0ab8f..6807439 100644
--- a/src/PPL/Sampling.hs
+++ b/src/PPL/Sampling.hs
@@ -2,6 +2,7 @@
 {-# LANGUAGE BlockArguments #-}
 {-# LANGUAGE LambdaCase #-}
 {-# LANGUAGE ViewPatterns #-}
+{-# LANGUAGE TupleSections #-}
 
 module PPL.Sampling
   ( mh,
@@ -24,25 +25,52 @@ import qualified Streaming as S
 import Streaming.Prelude (Of, Stream, yield)
 import System.IO.Unsafe
 import System.Random (StdGen, random, randoms)
-
-mh :: (Monad m) => StdGen -> Double -> Meas a -> Stream (Of (a, Log Double)) m ()
-mh g p m = step t0 t x w
+import qualified System.Random as R
+import Data.IORef
+import Control.Monad
+import Debug.Trace
+
+mh :: (MonadIO m) => StdGen -> Double -> Meas a -> Stream (Of (a, Log Double)) m ()
+mh g p m = do
+    let (g0, g1) = R.split g
+    omega <- liftIO $ newIORef (mempty, g0)
+    let (x, w) = head $ samples m $ newTree omega
+    step g1 omega x w
   where
-    (t0, t) = split $ randomTree g
-    (x, w) = head $ samples m t
 
-    step !t0 !t !x !w = do
-      let (t1 : t2 : t3 : t4 : _) = splitTrees t0
-          t' = mutateTree p t1 t2 t
-          (x', w') = head $ samples m t'
+    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
-          (Exp . log -> r) = draw t3
-          (t'', x'', w'') =
+          (omega'', x'', w'') =
             if r < ratio
-              then (t', x', w')
-              else (t, x, w)
+              then (omega', x', w')
+              else (omega, x, w)
       yield (x'', w'')
-      step t4 t'' x'' w''
+      step g2 omega'' x'' w''
+    
+    mutate :: MonadIO m => StdGen -> IORef (M.Map [Int] Double, StdGen) ->  m (IORef (M.Map [Int] Double, StdGen))
+    mutate g omega = do
+      (m, g0) <- liftIO $ readIORef omega
+      let (r:q:_) = R.randoms g
+          ks = M.keys m
+          k = ks !! floor (r * join traceShow (fromIntegral (length ks)))
+          m' = M.insert k q m
+      liftIO $ newIORef $ (m',g0)
+      
+     where
+      go x = do
+        g <- get
+        let (r, g1) = R.random g
+            (y, g2) = R.random g1
+        if r < p 
+          then do
+            put g2
+            pure y
+          else do
+            put g1
+            pure x
 
 -- Single site MH