packages feed

hakaru (empty) → 0.1

raw patch · 21 files changed

+1658/−0 lines, 21 filesdep +aesondep +basedep +bytestringsetup-changed

Dependencies added: aeson, base, bytestring, cassava, containers, hmatrix, logfloat, math-functions, parsec, pretty, random, statistics, text, vector, zlib

Files

+ Examples/Examples.hs view
@@ -0,0 +1,44 @@+{-# LANGUAGE RankNTypes, DataKinds, NoMonomorphismRestriction, BangPatterns #-}++module Examples where++import Types+import Data.Dynamic+import Control.Monad++import InterpreterMH hiding (main)+import Visual++bayesian_polynomial_regression = undefined++sparse_linear_regression = undefined++logistic_regression = undefined++outlier_detection = undefined++change_point_model = undefined++friends_who_smoke = undefined++latent_dirichelt_allocation = undefined++categorical_mixture = undefined++gaussian_mixture = undefined++naive_bayes = undefined++hidden_markov_model = undefined++matrix_factorization = undefined++rvm = undefined++item_response_theory = undefined++gaussian_process = undefined++hawkes_process = undefined++bayesian_neural_network = undefined
+ Examples/Tests.hs view
@@ -0,0 +1,106 @@+{-# LANGUAGE RankNTypes, NoMonomorphismRestriction, BangPatterns #-}++module Tests where++import Types+import Data.Dynamic+import Language.Hakaru.ImportanceSampler++-- Some example/test programs in our language+test :: Measure Bool+test = do+  c <- unconditioned (bern 0.5)+  _ <- conditioned (ifThenElse c (normal (lit (1 :: Double)) (lit 1))+                                 (uniformC (lit 0) (lit 3)))+  return c++test_dup :: Measure (Bool, Bool)+test_dup = do+  let c = unconditioned (bern 0.5)+  x <- c+  y <- c+  return (x,y)++test_dbn :: Measure Bool+test_dbn = do+  s0 <- unconditioned (bern 0.75)+  s1 <- unconditioned (if s0 then bern 0.75 else bern 0.25)+  _  <- conditioned (if s1 then bern 0.9 else bern 0.1)+  s2 <- unconditioned (if s1 then bern 0.75 else bern 0.25)+  _  <- conditioned (if s2 then bern 0.9 else bern 0.1)+  return s2++test_hmm :: Integer -> Measure Bool+test_hmm n = do+  s <- unconditioned (bern 0.75) +  loop_hmm n s++loop_hmm :: Integer -> (Bool -> Measure Bool)+loop_hmm !numLoops s = do+    _ <- conditioned (if s then bern 0.9 else bern 0.1)+    u <- unconditioned (if s then bern 0.75 else bern 0.25)+    if (numLoops > 1) then loop_hmm (numLoops - 1) u +                      else return s++test_carRoadModel :: Measure (Double, Double)+test_carRoadModel = do+  speed <- unconditioned (uniformC (lit (5::Double)) (lit (15::Double)))+  let z0 = lit 0 +  _ <- conditioned (normal (z0 :: Double) (lit 1))+  z1 <- unconditioned (normal (z0 + speed) (lit 1))+  _ <- conditioned (normal z1 (lit 1))+  z2 <- unconditioned (normal (z1 + speed) (lit 1))	+  _ <- conditioned (normal z2 (lit 1))+  z3 <- unconditioned (normal (z2 + speed) (lit 1))	+  _ <- conditioned (normal z3 (lit 1))+  z4 <- unconditioned (normal (z3 + speed) (lit 1))	+  return (z4, z3)++test_categorical :: Measure Bool+test_categorical = do +  rain <- unconditioned (categorical [(lit True, 0.2), (lit False, 0.8)]) +  sprinkler <- unconditioned (if rain then bern 0.01 else bern 0.4)+  _ <- conditioned (if rain then (if sprinkler then bern 0.99 else bern 0.8)+	                else (if sprinkler then bern 0.9 else bern 0.1))+  return rain++-- printing test results+main :: IO ()+main = sample_ 3 test conds >>+       putChar '\n' >>+       sample 1000 test conds >>=+       print+  where conds = [Lebesgue (toDyn (2 :: Double))]++main_dbn :: IO ()+main_dbn = sample_ 10 test_dbn conds >>+           putChar '\n' >>+           sample 1000 test_dbn conds >>=+           print +  where conds = [Discrete (toDyn (True :: Bool)),+                 Discrete (toDyn (True :: Bool))]++main_hmm :: IO ()+main_hmm = sample_ 10 (test_hmm 2) conds >>+           putChar '\n' >>+           sample 1000 (test_hmm 2) conds >>=+           print +  where conds = [Discrete (toDyn (True :: Bool)),+                 Discrete (toDyn (True :: Bool))]++main_carRoadModel :: IO ()+main_carRoadModel = sample_ 10 test_carRoadModel conds >>+                    putChar '\n' >>+                    sample 1000 test_carRoadModel conds >>=+                    print +  where conds = [Lebesgue (toDyn (0 :: Double)),+                 Lebesgue (toDyn (11 :: Double)), +                 Lebesgue (toDyn (19 :: Double)),+                 Lebesgue (toDyn (33 :: Double))]++main_categorical :: IO ()+main_categorical = sample_ 10 test_categorical conds >>+           putChar '\n' >>+           sample 1000 test_categorical conds >>=+           print +  where conds = [Discrete (toDyn (True :: Bool))]
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, The Hakaru Team++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of The Hakaru Team nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Language/Hakaru/ImportanceSampler.hs view
@@ -0,0 +1,174 @@+{-# LANGUAGE RankNTypes, NoMonomorphismRestriction, BangPatterns #-}+{-# OPTIONS -W #-}++module Language.Hakaru.ImportanceSampler where++-- This is an interpreter that's like Interpreter except conditioning is+-- checked at run time rather than by static types.  In other words, we allow+-- models to be compiled whose conditioned parts do not match the observation+-- inputs.  In exchange, we get to make Measure an instance of Monad, and we+-- can express models whose number of observations is unknown at compile time.++import Types (Cond(..), CSampler(CSampler))+import RandomChoice (normal_rng, chooseIndex)+import Mixture (Prob, empty, point, Mixture(..))+import Sampler (Sampler, deterministic, smap, sbind)++import System.Random+import Data.Monoid+import Data.Ix+import Data.Dynamic+import Data.List+import Control.Monad+import qualified Data.Map.Strict as M++import qualified Data.Number.LogFloat as LF++dirac :: (Eq a, Typeable a) => a -> CSampler a+dirac x = CSampler c where+  c Unconditioned = deterministic (point x 1)+  c (Discrete y) = case fromDynamic y of+    Just y  -> deterministic (if x == y then point x 1 else empty)+    Nothing -> error "dirac: did not get data from dynamic source"+  c _ = error "dirac: got a non-discrete sampler"++bern :: Double -> CSampler Bool+bern theta | 0 <= theta && theta <= 1 = CSampler c where+  c Unconditioned = \g0 -> case randomR (0, 1) g0 of+    (x, g) -> (point (x <= theta) 1, g)+  c (Discrete y) = case fromDynamic y of+    Just y -> deterministic (point y (LF.logFloat (if y then theta else 1 - theta)))+    Nothing -> error "bern: did not get data from dynamic source"+  c _ = error "bern: got a non-discrete sampler"+bern theta = error ("bernoulli: invalid parameter " ++ show theta)++uniformC :: (Fractional a, Real a, Random a, Typeable a) => a -> a -> CSampler a+uniformC lo hi | lo < hi = CSampler c where+  c Unconditioned = \g0 -> case randomR (lo,hi) g0 of+    (x, g) -> (point x 1, g)+  c (Lebesgue y) = case fromDynamic y of+    Just y -> deterministic (if lo < y && y < hi then point y density else empty)+    Nothing -> error "uniformC: did not get data from dynamic source"+  c _ = error "uniformC: got a discrete sampler"+  density = fromRational (toRational (recip (hi - lo)))+uniformC _ _ = error "uniformC: invalid parameters"++uniformD :: (Ix a, Random a, Typeable a) => a -> a -> CSampler a+uniformD lo hi | lo <= hi = CSampler c where+  c Unconditioned = \g0 -> case randomR (lo,hi) g0 of+    (x, g) -> (point x 1, g)+  c (Discrete y) = case fromDynamic y of+    Just y -> deterministic (if lo <= y && y <= hi then point y density else empty)+    Nothing -> error "uniformD: did not get data from dynamic source"+  c _ = error "uniformD: got a non-discrete sampler"+  density = recip (fromInteger (toInteger (rangeSize (lo,hi))))+uniformD _ _ = error "uniformD: invalid parameters"++poisson :: (Integral a, Typeable a) => Double -> CSampler a+poisson !l | 0 <= l = CSampler c where+  c Unconditioned = \g0 ->+    let probs = exp (-l) : zipWith (\k p -> p * l / k) [1..] probs+        (k, g) = chooseIndex probs g0+    in (point (fromInteger (toInteger k)) 1, g)+  c (Discrete k) = case fromDynamic k of+    Just k ->+      deterministic+        (if 0 <= k then point k (LF.logToLogFloat (-l)+                                 * LF.logFloat l ^ k+                                 / product (map fromIntegral [1..k]))+                   else empty)+    Nothing -> error "poisson: did not get data from dynamic source"+  c _ = error "poisson: got a non-discrete sampler"+poisson _ = error "poisson: invalid parameter"++normal :: (Real a, Floating a, Random a, Typeable a) => a -> a -> CSampler a+normal !mean !std | std > 0 = CSampler c where+  c Unconditioned = \g0 -> let (x, g) = normal_rng mean std g0+                           in (point (mean + std * x) 1, g)+  c (Lebesgue y) = case fromDynamic y of+    Just y ->+      let density  = exp (square ((y - mean) / std) / (-2)) / std / sqrt (2 * pi)+          square y = y * y+      in deterministic (point y (fromRational (toRational density))) -- TODO: use log-density and LogFloat directly+    Nothing -> error "normal: did not get data from dynamic source"+  c _ = error "normal: got a discrete sampler"+normal _ _ = error "normal: invalid parameters"++categorical :: (Typeable a, Eq a) => [(a, Prob)] -> CSampler a+categorical list = CSampler c where+  peak :: LF.LogFloat+  peak = maximum (map snd list)+  total :: Double+  (total, list') = mapAccumL f 0 list+    where f acc (a,b) = (acc', (a, (b', acc')))+            where b' = b/peak+                  acc' :: Double+                  acc' = acc + LF.fromLogFloat b'+  c Unconditioned =+    \g0 -> let (p, g) = randomR (0, total) g0+               (elem, _) : _ = filter (\(_,(_,p0)) -> p <= p0) list' in+           (point elem 1, g)+  c (Discrete y) = case fromDynamic y of+    Just y -> deterministic (maybe empty (point y . (/ LF.logFloat total) . fst)+                                         (lookup y list'))+    Nothing -> error "categorical: did not get data from dynamic source"+  c _ = error "categorical: got a non-discrete sampler"++-- Conditioned sampling+newtype Measure a = Measure { unMeasure :: [Cond] -> Sampler (a, [Cond]) }++bind :: Measure a -> (a -> Measure b) -> Measure b+bind measure continuation =+  Measure (\conds ->+    sbind (unMeasure measure conds)+          (\(a,conds) -> unMeasure (continuation a) conds))++instance Monad Measure where+  return x = Measure (\conds -> deterministic (point (x,conds) 1))+  (>>=)    = bind++conditioned, unconditioned :: CSampler a -> Measure a+conditioned   (CSampler f) = Measure (\(cond:conds) -> smap (\a->(a,conds)) (f cond         ))+unconditioned (CSampler f) = Measure (\      conds  -> smap (\a->(a,conds)) (f Unconditioned))++factor :: Prob -> Measure ()+factor p = Measure (\conds -> deterministic (point ((), conds) p))++-- Our language also includes the usual goodies of a lambda calculus+var :: a -> a+var = id++lit :: a -> a+lit = id++lam :: (a -> b) -> (a -> b)+lam f = f++app :: (a -> b) -> a -> b+app f x = f x++fix :: ((a -> b) -> (a -> b)) -> (a -> b)+fix g = f where f = g f++ifThenElse :: Bool -> a -> a -> a+ifThenElse True  t _ = t+ifThenElse False _ e = e++-- Drivers for testing+finish :: Mixture (a, [Cond]) -> Mixture a+finish (Mixture m) = Mixture (M.mapKeysMonotonic (\(a,[]) -> a) m)++sample :: (Ord a) => Int -> Measure a -> [Cond] -> IO (Mixture a)+sample !n measure conds = go n empty where+  once = getStdRandom (unMeasure measure conds)+  go 0 m = return m+  go n m = once >>= \result -> go (n - 1) $! mappend m (finish result)++sample_ :: (Ord a, Show a) => Int -> Measure a -> [Cond] -> IO ()+sample_ !n measure conds = replicateM_ n (once >>= pr) where+  once = getStdRandom (unMeasure measure conds)+  pr   = print . finish++logit :: Floating a => a -> a+logit !x = 1 / (1 + exp (- x))+
+ Language/Hakaru/Metropolis.hs view
@@ -0,0 +1,294 @@+{-# LANGUAGE RankNTypes, NoMonomorphismRestriction, BangPatterns,+  DeriveDataTypeable, GADTs, ScopedTypeVariables,+  ExistentialQuantification, StandaloneDeriving #-}++module Language.Hakaru.Metropolis where++import System.Random (RandomGen, StdGen, randomR, getStdGen)+import System.IO++import Control.Monad+import Data.Dynamic+import Data.Function (on)+import Data.Maybe++import qualified Data.Map.Strict as M++import RandomChoice+import Visual++{-++Shortcomings of this implementation++* uses parent-conditional sampling for proposal distribution+* re-evaluates entire program at every sample+* lacks way to block sample groups of variables++-}++type DistVal = Dynamic++data Dist a = Dist {logDensity :: a -> Likelihood,+                    sample :: forall g. RandomGen g => g -> (a, g)}+deriving instance Typeable1 Dist+ +data XRP = forall e. Typeable e => XRP (e, Dist e)++unXRP :: Typeable a => XRP -> Maybe (a, Dist a)+unXRP (XRP (e,f)) = cast (e,f)++type Var a = Int++type Likelihood = Double+type Visited = Bool+type Observed = Bool+type Cond = Maybe DistVal++type Subloc = Int+type Name = [Subloc]+type Database = M.Map Name (XRP, Likelihood, Visited, Observed)+newtype Measure a = Measure {unMeasure :: (RandomGen g) =>+                              (Name+                              ,Database+                              ,(Likelihood, Likelihood)+                              ,[Cond]+                              ,g+                           ) -> (a+                                ,Database+                                ,(Likelihood, Likelihood)+                                ,[Cond]+                                ,g)}+  deriving (Typeable)++-- n  is structural_name+-- d  is database+-- ll is likelihood of expression+-- conds is the observed data+-- g  is the random seed+++lit :: (Eq a, Typeable a) => a -> a+lit = id++return_ :: a -> Measure a+return_ x = Measure (\ (n, d, l, conds, g) -> (x, d, l, conds, g))++makeXRP :: (Typeable a, RandomGen g) => Cond -> Dist a+        -> Name -> Database -> g+        -> (a, Database, Likelihood, Likelihood, g)+makeXRP obs dist' n db g =+    case M.lookup n db of+      Just (xd, lb, b, ob) ->+        let Just (xb, dist) = unXRP xd+            (x,l) = case obs of+                      Just xd ->+                          let Just x = fromDynamic xd+                          in (x, logDensity dist x)+                      Nothing -> (xb, lb)+            l' = logDensity dist' x+            d1 = M.insert n (XRP (x,dist),+                             l',+                             True,+                             ob) db+        in (x, d1, l', 0, g)+      Nothing ->+        let (xnew, l, g1) = case obs of+             Just xdnew ->+                 let Just xnew = fromDynamic xdnew+                 in (xnew, logDensity dist' xnew, g)+             Nothing ->+                 (xnew, logDensity dist' xnew, g1)+                where (xnew, g1) = sample dist' g+            d1 = M.insert n (XRP (xnew, dist'),+                             l,+                             True,+                             isJust obs) db+        in (xnew, d1, l, l, g1)++updateLikelihood :: (Typeable a, RandomGen g) => +                    Likelihood -> Likelihood ->+                    (a, Database, Likelihood, Likelihood, g) ->+                    [Cond] ->+                    (a, Database, (Likelihood, Likelihood), [Cond], g)+updateLikelihood llTotal llFresh (x,d,l,lf,g) conds =+    (x, d, (llTotal+l, llFresh+lf), conds, g)++dirac :: (Eq a, Typeable a) => a -> Cond -> Measure a+dirac theta obs = Measure $ \(n, d, (llTotal,llFresh), conds, g) ->+    let dist' = Dist {logDensity = (\ x -> if x == theta then 0 else log 0),+                      sample = (\ g -> (theta,g))}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++bern :: Double -> Cond -> Measure Bool+bern p obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->+    let dist' = Dist {logDensity = (\ x -> log (if x then p else 1 - p)),+                      sample = (\ g -> case randomR (0, 1) g of+                                         (t, g') -> (t <= p, g'))}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++poisson :: Double -> Cond -> Measure Int+poisson l obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->+    let poissonLogDensity l x | l > 0 && x> 0 = (fromIntegral x)*(log l) - lnFact x - l+        poissonLogDensity l x | x==0 = -l+        poissonLogDensity _ _ = log 0+        dist' = Dist {logDensity = poissonLogDensity l,+                      sample = poisson_rng l}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++gamma :: Double -> Double -> Cond -> Measure Double+gamma shape scale obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->+    let dist' = Dist {logDensity = gammaLogDensity shape scale,+                      sample = gamma_rng shape scale}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++beta :: Double -> Double -> Cond -> Measure Double+beta a b obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->+    let dist' = Dist {logDensity = betaLogDensity a b,+                      sample = beta_rng a b}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++uniform :: Double -> Double -> Cond -> Measure Double+uniform lo hi obs = Measure $ \(n, d, (llTotal,llFresh), conds, g) ->+    let uniformLogDensity lo hi x | lo <= x && x <= hi = log (recip (hi - lo))+        uniformLogDensity _ _  x = log 0+        dist' = Dist {logDensity = uniformLogDensity lo hi,+                      sample = (\ g -> randomR (lo, hi) g)}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++normal :: Double -> Double -> Cond -> Measure Double+normal mu sd obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->+    let dist' = Dist {logDensity = normalLogDensity mu sd,+                      sample = normal_rng mu sd}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++laplace :: Double -> Double -> Cond -> Measure Double+laplace mu sd obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->+    let dist' = Dist {logDensity = laplaceLogDensity mu sd,+                      sample = laplace_rng mu sd}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++categorical :: (Eq a, Typeable a) => [(a,Double)] +            -> Cond -> Measure a+categorical list obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->+    let categoricalLogDensity list x = log $ fromMaybe 0 (lookup x list)+        categoricalSample list g = (elem, g1)+           where+             (p, g1) = randomR (0, total) g+             elem = fst $ head $ filter (\(_,p0) -> p <= p0) sumList+             sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list+             total = sum $ map snd list+        dist' = Dist {logDensity = categoricalLogDensity list,+                      sample = categoricalSample list}+        xrp = makeXRP obs dist' n d g+    in updateLikelihood llTotal llFresh xrp conds++factor :: Likelihood -> Measure ()+factor l = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->+   ((), d, (llTotal + l, llFresh), conds, g)++resample :: RandomGen g => XRP -> g ->+            (XRP, Likelihood, Likelihood, Likelihood, g)+resample (XRP (x, dist)) g =+    let (x', g1) = sample dist g+        fwd = logDensity dist x'+        rvs = logDensity dist x+        l' = fwd+    in (XRP (x', dist), l', fwd, rvs, g1)++bind :: Measure a -> (a -> Measure b) -> Measure b+bind (Measure m) cont = Measure $ \ (n,d,ll,conds,g) ->+    let (v, d1, ll1, conds1, g1) = m (0:n, d, ll, conds, g)+    in unMeasure (cont v) (1:n, d1, ll1, conds1, g1)++conditioned :: (Cond -> Measure a) -> Measure a+conditioned f = Measure $ \ (n,d,ll,cond:conds,g) ->+    unMeasure (f cond) (n, d, ll, conds, g)++unconditioned :: (Cond -> Measure a) -> Measure a+unconditioned f = f Nothing++instance Monad Measure where+  return = return_+  (>>=)    = bind++lam :: (a -> b) -> (a -> b)+lam f = f++app :: (a -> b) -> a -> b+app f x = f x++fix :: ((a -> b) -> (a -> b)) -> (a -> b)+fix g = f where f = g f++ifThenElse :: Bool -> a -> a -> a+ifThenElse True  t _ = t+ifThenElse False _ f = f++run :: Measure a -> [Cond] -> IO (a, Database, Likelihood)+run (Measure prog) conds = do+  g <- getStdGen+  let (v, d, ll, conds1, g') =+          prog ([0], M.empty, (0,0), conds, g)+  return (v, d, fst ll)++traceUpdate :: RandomGen g => Measure a -> Database -> [Cond] -> g+            -> (a, Database, Likelihood, Likelihood, Likelihood, g)+traceUpdate (Measure prog) d conds g = do+  let d1 = M.map (\ (x, l, _, ob) -> (x, l, False, ob)) d+  let (v, d2, (llTotal, llFresh), conds1, g1) =+          prog ([0], d1, (0,0), conds, g)+  let (d3, stale_d) = M.partition (\ (_, _, v, _) -> v) d2+  let llStale = M.foldl' (\ llStale (_,l,_,_) -> llStale + l)+                0 stale_d+  (v, d3, llTotal, llFresh, llStale, g1)++initialStep :: Measure a -> [Cond] ->+               IO (a, Database,+                   Likelihood, Likelihood, Likelihood, StdGen)+initialStep prog conds = do+  g <- getStdGen+  return $ traceUpdate prog M.empty conds g++-- TODO: Make a way of passing user-provided proposal distributions+updateDB :: (RandomGen g) => +            Name -> Database -> Observed -> XRP -> g+         -> (Database, Likelihood, Likelihood, Likelihood, g)+updateDB name db ob xd g = (db', l', fwd, rvs, g)+    where db' = M.insert name (x', l', True, ob) db+          (x', l', fwd, rvs, g1) = resample xd g++transition :: (Typeable a, RandomGen g) => Measure a -> [Cond]+           -> a -> Database -> Likelihood -> g -> [a]+transition prog conds v db ll g =+  let dbSize = M.size db+      -- choose an unconditioned choice+      (condDb, uncondDb) = M.partition (\ (_, _, _, ob) -> ob) db+      (choice, g1) = randomR (0, (M.size uncondDb) -1) g+      (name, (xd, l, _, ob))  = M.elemAt choice uncondDb+      (db', l', fwd, rvs, g2) = updateDB name db ob xd g1+      (v', db2, llTotal, llFresh, llStale, g3) = traceUpdate prog db' conds g2+      a = llTotal - ll+          + rvs - fwd+          + log (fromIntegral dbSize) - log (fromIntegral $ M.size db2)+          + llStale - llFresh+      (u, g4) = randomR (0 :: Double, 1) g3 in++  if (log u < a) then+      v' : (transition prog conds v' db2 llTotal g4)+  else+      v : (transition prog conds v db ll g4)++mcmc :: Typeable a => Measure a -> [Cond] -> IO [a]+mcmc prog conds = do+  (v, d, llTotal, llFresh, llStale, g) <- initialStep prog conds+  return $ transition prog conds v d llTotal g+
+ Language/Hakaru/Symbolic.hs view
@@ -0,0 +1,78 @@+{-# LANGUAGE GADTs, TypeFamilies #-}+{-# OPTIONS -W #-}++module Language.Hakaru.Symbolic where++data Prob+data Measure a+data Dist a++-- Symbolic AST (from Syntax.hs)+class Symbolic repr where+  real 			    :: Double -> repr Double+  bool 			    :: Bool -> repr Bool+  add, minus, mul, exp  	:: repr Double -> repr Double -> repr Double+  sqrt, cos, sin	:: repr Double -> repr Double+  bind	 		    :: repr (Measure a) -> (repr a -> repr (Measure a)) +	   		           -> repr (Measure a)+  ret 			    :: repr a -> repr (Measure a)+  uniformD, uniformC, normal :: repr Double -> repr Double -> repr (Dist Double)+  conditioned, unconditioned :: repr (Dist a) -> repr (Measure a)++-- Printer (to Maple)+type VarCounter = Int+newtype Maple a = Maple { unMaple :: Bool -> VarCounter -> String }++instance Symbolic Maple where+  real x 	= Maple $ \_ _ -> show x+  bool x 	= Maple $ \_ _ -> show x+  add e1 e2 	= Maple $ \f h -> unMaple e1 f h ++ "+" ++ unMaple e2 f h+  minus e1 e2 	= Maple $ \f h -> unMaple e1 f h ++ "-" ++ unMaple e2 f h  +  mul e1 e2 	= Maple $ \f h -> unMaple e1 f h ++ "*" ++ unMaple e2 f h+  exp e1 e2     = Maple $ \f h -> unMaple e1 f h ++ "^" ++ unMaple e2 f h+  sqrt e	= Maple $ \f h -> "sqrt(" ++ unMaple e f h ++ ")"+  cos e		= Maple $ \f h -> "cos(" ++ unMaple e f h ++ ")"+  sin e		= Maple $ \f h -> "sin(" ++ unMaple e f h ++ ")"+  bind m c 	= Maple $ \f h -> unMaple m True h ++ +  		          unMaple (c (Maple $ \_ _ -> ("x" ++ show h))) (f) (succ h)+  		          ++ unMaple m False h +  uniformC e1 e2 = Maple $ \f h -> if f == True then  +		          show (1/((read (unMaple e2 f h) :: Double) - +		          (read (unMaple e1 f h) :: Double))) ++ " * Int (" +		          else ", x" ++ show h ++ "=" ++ unMaple e1 f h ++ ".." ++ +		          unMaple e2 f h ++ ")"+  uniformD e1 e2 = Maple $ \f h -> if f == True then  +		          show (1/((read (unMaple e2 f h) :: Double) - +		          (read (unMaple e1 f h) :: Double))) ++ " * Sum (" +		          else ", x" ++ show h ++ "=" ++ unMaple e1 f h ++ ".." ++ +		          unMaple e2 f h ++ ")"				  +  normal e1 e2 	= Maple $ \f h -> if f == True then  +		          "Int (PDF (Normal (" ++ unMaple e1 f h ++ ", " +++		          unMaple e2 f h ++ ", x" ++ show h ++ ") * "  +		          else ", x" ++ show h ++ "=" ++ unMaple e1 f h ++ ".." ++ +		          unMaple e2 f h ++ ")"			  +  unconditioned e = Maple $ \f h -> unMaple e f h+  conditioned e   = Maple $ \f h -> unMaple e f h  +  ret e 	      = Maple $ \f h -> "g(" ++ unMaple e f h ++ ")"++view e = unMaple e True 0++-- TEST CASES+exp1 = unconditioned (uniformC (real 1) (real 3)) `bind` \s ->+       ret s++-- Borel's Paradox Simplified+exp2 = unconditioned (uniformD (real 1) (real 3)) `bind` \s ->+       unconditioned (uniformC (real (-1)) (real 1)) `bind` \x ->+       let y = s `mul` x in ret y++-- Borel's Paradox+exp3 = unconditioned (uniformD (real 1) (real 2)) `bind` \s ->+       unconditioned (uniformC (real (-1)) (real 1)) `bind` \x ->+       let y = (InterpreterSymbolic.sqrt ((real 1 ) `minus` +			   (InterpreterSymbolic.exp s (real 2)))) `mul`+	           (InterpreterSymbolic.sin x) in ret y  ++test = view exp1+test2 = view exp2+test3 = view exp3
+ Mixture.hs view
@@ -0,0 +1,61 @@+{-# LANGUAGE RankNTypes, BangPatterns #-}+{-# OPTIONS -W #-}++module Mixture (Prob, point, empty, scale,+  Mixture(..), toList, mnull, mmap, cross, mode) where++import Data.Monoid+import Data.Ord (comparing)+import Data.List (maximumBy)+import qualified Data.Map.Strict as M+import Data.Number.LogFloat hiding (isInfinite)+import Text.Show (showListWith)+import Numeric (showFFloat)++type Prob = LogFloat++-- Mixtures (the results of importance sampling)++newtype Mixture k = Mixture { unMixture :: M.Map k Prob }++instance (Show k) => Show (Mixture k) where+  showsPrec d (Mixture m) = showParen (d > 0) $+    showString "Mixture $ fromList " . showListWith s (M.toList m)+    where s (k,p) = showChar '('+                  . shows k+                  . showChar ','+                  . (if isInfinite l || -42 < l && l < 42+                     then showFFloat Nothing (fromLogFloat p :: Double)+                     else showString "logToLogFloat " . showsPrec 11 l)+                  . showChar ')'+            where l = logFromLogFloat p :: Double++instance (Ord k) => Monoid (Mixture k) where+  mempty        = empty+  mappend m1 m2 = Mixture (M.unionWith (+) (unMixture m1) (unMixture m2))+  mconcat ms    = Mixture (M.unionsWith (+) (map unMixture ms))++empty :: Mixture k+empty = Mixture M.empty++toList :: Mixture k -> [(k, Prob)]+toList = M.toList . unMixture++mnull :: Mixture k -> Bool+mnull = all (0>=) . M.elems . unMixture++point :: k -> Prob -> Mixture k+point k !v = Mixture (M.singleton k v)++scale :: Prob -> Mixture k -> Mixture k+scale !v = Mixture . M.map (v *) . unMixture++mmap :: (Ord k2) => (k1 -> k2) -> Mixture k1 -> Mixture k2+mmap f = Mixture . M.mapKeysWith (+) f . unMixture++cross :: (Ord k) => (k1 -> k2 -> k) -> Mixture k1 -> Mixture k2 -> Mixture k+cross f m1 m2 = mconcat [ mmap (`f` k) (scale v m1)+                        | (k,v) <- M.toList (unMixture m2) ]++mode :: Mixture k -> (k, Prob)+mode (Mixture m) = maximumBy (comparing snd) (M.toList m)
+ RandomChoice.hs view
@@ -0,0 +1,145 @@+{-# LANGUAGE BangPatterns #-}+module RandomChoice where++import System.Random+import Mixture+import Data.Maybe (fromMaybe)+import Data.List (findIndex, foldl')+import Numeric.SpecFunctions+import qualified Data.Vector.Unboxed as U+import qualified Data.Map.Strict as M+import qualified Data.Number.LogFloat as LF++marsaglia :: (RandomGen g, Random a, Ord a, Floating a) => g -> ((a, a), g)+marsaglia g0 = -- "Marsaglia polar method"+  let (x, g1) = randomR (-1,1) g0+      (y, g ) = randomR (-1,1) g1+      s       = x * x + y * y+      q       = sqrt ((-2) * log s / s)+  in if 1 >= s && s > 0 then ((x * q, y * q), g) else marsaglia g++choose :: (RandomGen g) => Mixture k -> g -> (k, Prob, g)+choose (Mixture m) g0 =+  let peak = maximum (M.elems m)+      unMix = M.map (LF.fromLogFloat . (/peak)) m+      total = M.foldl' (+) (0::Double) unMix+      (p, g) = randomR (0, total) g0+      f !k !v b !p0 = let p1 = p0 + v in if p <= p1 then k else b p1+      err p0 = error ("choose: failure p0=" ++ show p0 +++                      " total=" ++ show total +++                      " size=" ++ show (M.size m))+  in (M.foldrWithKey f err unMix 0, LF.logFloat total * peak, g)++chooseIndex :: (RandomGen g) => [Double] -> g -> (Int, g)+chooseIndex probs g0 =+  let (p, g) = random g0+      k = fromMaybe (error ("chooseIndex: failure p=" ++ show p))+                    (findIndex (p <=) (scanl1 (+) probs))+  in (k, g)++normal_rng :: (Real a, Floating a, Random a, RandomGen g) =>+              a -> a -> g -> (a, g)+normal_rng mu sd g | sd > 0 = case marsaglia g of+                                ((x, _), g1) -> (mu + sd * x, g1)+normal_rng _ _ _ = error "normal: invalid parameters"++normalLogDensity mu sd x = (-tau * square (x - mu)+                            + log (tau / pi / 2)) / 2+  where square y = y * y+        tau = 1 / square sd++lnFact = logFactorial++-- Makes use of Atkinson's algorithm as described in:+-- Monte Carlo Statistical Methods pg. 55+--+-- Further discussion at:+-- http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/+poisson_rng :: (RandomGen g) => Double -> g -> (Int, g)+poisson_rng lambda g0 = make_poisson g0+   where smu = sqrt lambda+         b  = 0.931 + 2.53*smu+         a  = -0.059 + 0.02483*b+         vr = 0.9277 - 3.6224/(b - 2)+         arep  = 1.1239 + 1.1368/(b-3.4)+         lnlam = log lambda++         make_poisson :: (RandomGen g) => g -> (Int,g)+         make_poisson g = let (u, g1) = randomR (-0.5,0.5) g+                              (v, g2) = randomR (0,1) g1+                              us = 0.5 - abs u+                              k = floor $ (2*a / us + b)*u + lambda + 0.43 in+                          case (us, v, k) of+                            (us,v,k) | us >= 0.07 && v <= vr -> (k, g2)+                            (_,_, k) | k < 0 -> make_poisson g2+                            (us,v,k) | us <= 0.013 && v > us -> make_poisson g2+                            (us,v,k) | accept_region us v k -> (k, g2)+                            _        -> make_poisson g2++         accept_region us v k = log (v * arep / (a/(us*us)+b)) <=+                                -lambda + (fromIntegral k)*lnlam - lnFact k++-- Direct implementation of  "A Simple Method for Generating Gamma Variables"+-- by George Marsaglia and Wai Wan Tsang.+gamma_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)+gamma_rng shape scale g | shape <= 0.0  = error "gamma: got a negative shape paramater"+gamma_rng shape scale g | scale <= 0.0  = error "gamma: got a negative scale paramater"+gamma_rng shape scale g | shape <  1.0  = (gvar2, g2)+                      where (gvar1, g1) = gamma_rng (shape + 1) scale g+                            (w,  g2) = randomR (0,1) g1+                            gvar2 = scale * gvar1 * (w ** recip shape) +gamma_rng shape scale g = +    let d = shape - 1/3+        c = recip $ sqrt $ 9*d+        -- Algorithm recommends inlining normal generator+        n g = normal_rng 1 c g+        (v, g2) = until (\x -> fst x > 0.0) (\ (_, g) -> normal_rng 1 c g) (n g)+        x = (v - 1) / c+        sqr = x * x+        v3 = v * v * v+        (u, g3) = randomR (0.0, 1.0) g2+        accept  = u < 1.0 - 0.0331*(sqr*sqr) || log u < 0.5*sqr + d*(1.0 - v3 + log v3)+    in case accept of+         True -> (scale*d*v3, g3)+         False -> gamma_rng shape scale g3++gammaLogDensity shape scale x | x>= 0 && shape > 0 && scale > 0 =+     scale * log shape - scale * x + (shape - 1) * log x - logGamma shape+gammaLogDensity _ _ _ = log 0++beta_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)+-- Consider adding case for a <= 1 && b <= 1+beta_rng a b g = let (ga, g1) = gamma_rng a 1 g+                     (gb, g2) = gamma_rng b 1 g1+                 in (ga / (ga + gb), g2)++betaLogDensity a b x | x < 0 || x > 1 = error "beta: value must be between 0 and 1"+betaLogDensity a b x | a <= 0 || b <= 0 = error "beta: parameters must be positve" +betaLogDensity a b x = (logGamma (a + b)+                        - logGamma a+                        - logGamma b+                        + x * log (a - 1)+                        + (1 - x) * log (b - 1))++laplace_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)+laplace_rng mu sd g = sample (randomR (0.0, 1.0) g)+   where sample (u, g1) = case u < 0.5 of+                            True  -> (mu + sd * log (u + u), g1)+                            False -> (mu - sd * log (2.0 - u - u), g1)++laplaceLogDensity mu sd x = - log (2 * sd) - abs (x - mu) / sd++-- Consider having dirichlet return Vector+-- Note: This is acutally symmetric dirichlet+dirichlet_rng :: (RandomGen g) => Int ->  Double -> g -> ([Double], g)+dirichlet_rng n a g = normalize (gammas g n)+  where gammas g 0 = ([], 0, g)+        gammas g n = let (xs, total, g1) = gammas g (n-1)+                         ( x, g2) = gamma_rng a 1 g1 +                     in ((x : xs), x+total, g2)+        normalize (a, total, g) = (map (/ total) a, g)++dirichletLogDensity a x | all (> 0) x = sum (zipWith logTerm a x) + logGamma (sum a)+  where sum a = foldl' (+) 0 a+        logTerm a x = (a-1) * log x - logGamma a+dirichletLogDensity _ _ = error "dirichlet: all values must be between 0 and 1"
+ Sampler.hs view
@@ -0,0 +1,26 @@+{-# LANGUAGE RankNTypes #-}+{-# OPTIONS -W #-}++module Sampler (Sampler, deterministic, sbind, smap) where++import Mixture (Mixture, mnull, empty, scale, point)+import RandomChoice (choose)+import System.Random (RandomGen)++-- Sampling procedures that produce one sample++type Sampler a = forall g. (RandomGen g) => g -> (Mixture a, g)++deterministic :: Mixture a -> Sampler a+deterministic m g = (m, g)++sbind :: Sampler a -> (a -> Sampler b) -> Sampler b+sbind s k g0 =+  case s g0 of { (m1, g1) ->+    if mnull m1 then (empty, g1) else+      case choose m1 g1 of { (a, v, g2) ->+        case k a g2 of { (m2, g) ->+          (scale v m2, g) } } }++smap :: (a -> b) -> Sampler a -> Sampler b+smap f s = sbind s (\a -> deterministic (point (f a) 1))
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ Syntax.hs view
@@ -0,0 +1,185 @@+{-# LANGUAGE TypeFamilies, ConstraintKinds, GADTs, FlexibleContexts #-}++module Syntax where++-- The syntax++import GHC.Exts (Constraint)++-- TODO: The pretty-printing semantics++import qualified Text.PrettyPrint as PP++-- The importance-sampling semantics++import Types (Cond, CSampler)+import Data.Dynamic (Typeable)+import qualified Data.Number.LogFloat as LF+import qualified InterpreterDynamic as IS++-- The Metropolis-Hastings semantics++import qualified InterpreterMH as MH++-- The syntax++data Prob+data Measure a+data Dist a++class Mochastic repr where+  type Type repr a :: Constraint+  real        :: Double -> repr Double+  bool        :: Bool -> repr Bool+  add, mul    :: repr Double -> repr Double -> repr Double+  neg         :: repr Double -> repr Double+  neg         =  mul (real (-1))+  logFloat, logToLogFloat+              :: repr Double -> repr Prob+  unbool      :: repr Bool -> repr c -> repr c+              -> repr c+  pair        :: repr a -> repr b -> repr (a, b)+  unpair      :: repr (a, b) -> (repr a -> repr b -> repr c)+              -> repr c+  inl         :: repr a -> repr (Either a b)+  inr         :: repr b -> repr (Either a b)+  uneither    :: repr (Either a b) -> (repr a -> repr c) -> (repr b -> repr c)+              -> repr c+  nil         :: repr [a]+  cons        :: repr a -> repr [a] -> repr [a]+  unlist      :: repr [a] -> repr c -> (repr a -> repr [a] -> repr c)+              -> repr c+  ret         :: repr a -> repr (Measure a)+  bind        :: repr (Measure a) -> (repr a -> repr (Measure b))+              -> repr (Measure b)+  conditioned, unconditioned :: repr (Dist a) -> repr (Measure a)+  factor      :: repr Prob -> repr (Measure ())+  dirac       :: (Type repr a) => repr a -> repr (Dist a)+  categorical :: (Type repr a) => repr [(a, Prob)] -> repr (Dist a)+  bern        :: (Type repr Bool) => repr Double -> repr (Dist Bool)+  bern p      =  categorical $+                 cons (pair (bool True) (logFloat p)) $+                 cons (pair (bool False) (logFloat (add (real 1) (neg p)))) $+                 nil+  normal, uniform+              :: repr Double -> repr Double -> repr (Dist Double)+  poisson     :: repr Double -> repr (Dist Int)++-- TODO: The initial (AST) "semantics"+-- (Hey Oleg, is there any better way to deal with the Type constraint, so that+-- the AST constructor doesn't have to take a repr constructor argument?)++data AST repr a where+  Real :: Double -> AST repr Double+  Unbool :: AST repr Bool -> AST repr c -> AST repr c -> AST repr c+  Categorical :: (Type repr a) => AST repr [(a, Prob)] -> AST repr (Dist a)+  -- ...++instance (Mochastic repr) => Mochastic (AST repr) where+  type Type (AST repr) a = Type repr a+  real = Real+  unbool = Unbool+  categorical = Categorical+  -- ...++eval :: (Mochastic repr) => AST repr a -> repr a+eval (Real x) = real x+eval (Unbool b x y) = unbool (eval b) (eval x) (eval y)+eval (Categorical xps) = categorical (eval xps)+-- ...++-- TODO: The pretty-printing semantics++newtype PP a = PP (Int -> PP.Doc)++-- The importance-sampling semantics++newtype IS a = IS (IS' a)+type family IS' a+type instance IS' (Measure a)  = IS.Measure (IS' a)+type instance IS' (Dist a)     = CSampler (IS' a)+type instance IS' [a]          = [IS' a]+type instance IS' (a, b)       = (IS' a, IS' b)+type instance IS' (Either a b) = Either (IS' a) (IS' b)+type instance IS' ()           = ()+type instance IS' Bool         = Bool+type instance IS' Double       = Double+type instance IS' Prob         = LF.LogFloat+type instance IS' Int          = Int++instance Mochastic IS where+  type Type IS a = (Eq (IS' a), Typeable (IS' a))+  real                    = IS+  bool                    = IS+  add (IS x) (IS y)       = IS (x + y)+  mul (IS x) (IS y)       = IS (x * y)+  neg (IS x)              = IS (-x)+  logFloat (IS x)         = IS (LF.logFloat x)+  logToLogFloat (IS x)    = IS (LF.logToLogFloat x)+  unbool (IS b) x y       = if b then x else y+  pair (IS x) (IS y)      = IS (x, y)+  unpair (IS (x, y)) c    = c (IS x) (IS y)+  inl (IS x)              = IS (Left x)+  inr (IS x)              = IS (Right x)+  uneither (IS e) c d     = either (c . IS) (d . IS) e+  nil                     = IS []+  cons (IS x) (IS xs)     = IS (x:xs)+  unlist (IS []) n c      = n+  unlist (IS (x:xs)) n c  = c (IS x) (IS xs)+  ret (IS x)              = IS (return x)+  bind (IS m) k           = IS (m >>= \x -> case k (IS x) of IS n -> n)+  conditioned (IS dist)   = IS (IS.conditioned dist)+  unconditioned (IS dist) = IS (IS.unconditioned dist)+  factor (IS p)           = IS (IS.factor p)+  dirac (IS x)            = IS (IS.dirac x)+  categorical (IS xps)    = IS (IS.categorical xps)+  bern (IS p)             = IS (IS.bern p)+  normal (IS m) (IS s)    = IS (IS.normal m s)+  uniform (IS lo) (IS hi) = IS (IS.uniformC lo hi)+  poisson (IS l)          = IS (IS.poisson l)++-- The Metropolis-Hastings semantics++newtype MH a = MH (MH' a)+type family MH' a+type instance MH' (Measure a)  = MH.Measure (MH' a)+type instance MH' (Dist a)     = MH.Cond -> MH.Measure (MH' a)+type instance MH' [a]          = [MH' a]+type instance MH' (a, b)       = (MH' a, MH' b)+type instance MH' (Either a b) = Either (MH' a) (MH' b)+type instance MH' ()           = ()+type instance MH' Bool         = Bool+type instance MH' Double       = Double+type instance MH' Prob         = MH.Likelihood+type instance MH' Int          = Int++instance Mochastic MH where+  type Type MH a = (Eq (MH' a), Typeable (MH' a), Show (MH' a))+  real                    = MH+  bool                    = MH+  add (MH x) (MH y)       = MH (x + y)+  mul (MH x) (MH y)       = MH (x * y)+  neg (MH x)              = MH (-x)+  logFloat (MH x)         = MH (LF.logFromLogFloat (LF.logFloat x))+  logToLogFloat (MH x)    = MH (LF.logFromLogFloat (LF.logToLogFloat x))+  unbool (MH b) x y       = if b then x else y+  pair (MH x) (MH y)      = MH (x, y)+  unpair (MH (x, y)) c    = c (MH x) (MH y)+  inl (MH x)              = MH (Left x)+  inr (MH x)              = MH (Right x)+  uneither (MH e) c d     = either (c . MH) (d . MH) e+  nil                     = MH []+  cons (MH x) (MH xs)     = MH (x:xs)+  unlist (MH []) n c      = n+  unlist (MH (x:xs)) n c  = c (MH x) (MH xs)+  ret (MH x)              = MH (return x)+  bind (MH m) k           = MH (m >>= \x -> case k (MH x) of MH n -> n)+  conditioned (MH dist)   = MH (MH.conditioned dist)+  unconditioned (MH dist) = MH (MH.unconditioned dist)+  factor (MH p)           = MH (MH.factor p)+  dirac (MH x)            = MH (MH.dirac x)+  categorical (MH xps)    = MH (MH.categorical xps)+  bern (MH p)             = MH (MH.bern p)+  normal (MH m) (MH s)    = MH (MH.normal m s)+  uniform (MH lo) (MH hi) = MH (MH.uniform lo hi)+  poisson                 = error "poisson: not implemented for MH" -- TODO
+ Types.hs view
@@ -0,0 +1,13 @@+{-# LANGUAGE RankNTypes, BangPatterns #-}+{-# OPTIONS -W #-}++module Types where++import Sampler (Sampler)++import Data.Dynamic++-- Basic types for conditioning and conditioned sampler+data Cond = Unconditioned | Lebesgue !Dynamic | Discrete !Dynamic+  deriving (Show)+newtype CSampler a = CSampler (Cond -> Sampler a)
+ Util/Coda.hs view
@@ -0,0 +1,12 @@+module Util.Coda where++import Statistics.Autocorrelation+import qualified Data.Packed.Vector as V+import qualified Data.Vector.Generic as G++effectiveSampleSize :: [Double] -> Double+effectiveSampleSize samples = n / (1 + 2*(G.sum rho))+  where n = fromIntegral (V.dim vec)+        vec = V.fromList samples+        cov = autocovariance vec+        rho = G.map (/ G.head cov) cov
+ Util/Csv.hs view
@@ -0,0 +1,40 @@+{-# LANGUAGE TypeOperators #-}++module Util.Csv ((:::)((:::)), decodeFile, decodeGZipFile,+                 decodeFileStream, decodeGZipFileStream) where++import Data.Csv ( HasHeader(..), FromRecord(..), FromField(..)+                , ToRecord(..), ToField(..), decode)+import qualified Data.Csv.Streaming as CS (decode)+import Codec.Compression.GZip (decompress)+import qualified Data.Foldable as F+import qualified Data.ByteString.Lazy as B+import qualified Data.Vector as V+import Control.Applicative ((<*>), (<$>))++data a ::: b = a ::: b deriving (Eq, Ord, Read, Show)+infixr 5 :::++instance (FromField a, FromRecord b) => FromRecord (a ::: b) where+  parseRecord v | V.null v  = fail "too few fields in input record"+                | otherwise = (:::) <$> parseField (V.head v) <*> parseRecord (V.tail v)++instance (ToField a, ToRecord b) => ToRecord (a ::: b) where+  toRecord (a ::: b) = V.cons (toField a) (toRecord b)++decodeBytes :: FromRecord a => B.ByteString -> [a]+decodeBytes bs = case decode HasHeader bs of+                   Left _ -> []+                   Right v -> V.toList v++decodeFile :: FromRecord a => FilePath -> IO [a]+decodeFile = fmap decodeBytes . B.readFile++decodeGZipFile :: FromRecord a => FilePath -> IO [a]+decodeGZipFile = fmap (decodeBytes . decompress) . B.readFile++decodeFileStream :: FromRecord a => FilePath -> IO [a]+decodeFileStream = fmap (F.toList . CS.decode HasHeader) . B.readFile++decodeGZipFileStream :: FromRecord a => FilePath -> IO [a]+decodeGZipFileStream = fmap (F.toList . CS.decode HasHeader . decompress) . B.readFile
+ Util/Extras.hs view
@@ -0,0 +1,84 @@+{-|+  Functions on lists and sequences.+  Some of the functions follow the style of Data.Random.Extras +  (from the random-extras package), but are written for use with+  PRNGs from System.Random rather than from the random-fu package.+-}++module Util.Extras where++import qualified Data.Sequence as S+import System.Random+import Data.Maybe+import qualified Data.Foldable as F++import Data.Dynamic+import Types++extract :: S.Seq a -> Int -> Maybe (S.Seq a, a)+extract s i | S.null r = Nothing+            | otherwise  = Just (a S.>< c, b)+    where (a, r) = S.splitAt i s +          (b S.:< c) = S.viewl r++randomExtract :: S.Seq a -> IO (Maybe (S.Seq a, a))+randomExtract s = do+  g <- newStdGen+  let (i,_) = randomR (0, S.length s - 1) g+  return $ extract s i++{-| +  Given a sequence, return a *sorted* sequence of+  n randomly selected elements from *distinct positions* in the sequence+-}++randomElems :: Ord a => S.Seq a -> Int -> IO (S.Seq a)+randomElems = randomElemsTR S.empty++randomElemsTR :: Ord a => S.Seq a -> S.Seq a -> Int -> IO (S.Seq a)+randomElemsTR ixs s n+    | n == S.length s = return $ S.unstableSort s+    | n == 1 = do (_,i) <- fmap fromJust (randomExtract s)+                  return.S.unstableSort $ i S.<| ixs+    | otherwise = do (s',i) <- fmap fromJust (randomExtract s)+                     (randomElemsTR $! (i S.<| ixs)) s' (n-1)++{-|+  Chop a sequence at the given indices. +  Assume number of indices given < length of sequence to be chopped+-}++pieces :: S.Seq a -> S.Seq Int -> [S.Seq a]+pieces s ixs = let f (ps,r,x) y = let (p,r') = S.splitAt (y-x) r+                                  in (p:ps,r',y)+                   g (a,b,_) = b:a+               in g $ F.foldl f ([],s,0) ixs++{-|+  Given n, chop a sequence at m random points+  where m = min (length-1, n-1)+-}++randomPieces :: Int -> S.Seq a -> IO [S.Seq a]+randomPieces n s+    | n >= l = return $ F.toList $ fmap S.singleton s+    | otherwise = do ixs <- randomElems (S.fromList [1..l-1]) (n-1)+                     return $ pieces s ixs+    where l = S.length s++{-|+  > pairs [1,2,3,4]+  [(1,2),(1,3),(1,4),(2,3),(2,4),(3,4)]+  > pairs [1,2,4,4]+  [(1,2),(1,4),(1,4),(2,4),(2,4),(4,4)]+-}++pairs :: [a] -> [(a,a)]+pairs [] = []+pairs (x:xs) = (zip (repeat x) xs) ++ pairs xs++l2Norm :: Floating a => [a] -> a+l2Norm l = sqrt.sum $ zipWith (*) l l++dataLoad []     = []+dataLoad (x:xs) = Lebesgue (toDyn (x :: Double)) : dataLoad xs
+ Util/FileInterpolater.hs view
@@ -0,0 +1,115 @@+{-# LANGUAGE RankNTypes, NoMonomorphismRestriction #-}+{-# OPTIONS -W #-}++module Main where++import System.IO+import Text.ParserCombinators.Parsec+import System.Environment (getArgs, getProgName)+import System.Exit (exitFailure)++data ControlMeas = CM { time1 :: Double, vel :: Double, steer :: Double }+data Sensor = S { time2 :: Double, lat :: Double, long :: Double, orient :: Double }+data Merged = M {cm :: ControlMeas, sens :: Sensor}++type Table1 = [ ControlMeas ]+type Table2 = [ Sensor ]+type MergedT = [ Merged ]++-- prev will be Nothing on startup+data Tracking a = Tr {prev :: Maybe a, curr :: a, rest :: [a]}++-- interpolate and merge 2 files+main :: IO ()+main = do+  args <- getArgs+  case args of+    [fileName1, fileName2] -> do+      handle1 <- openFile fileName1 ReadMode+      handle2 <- openFile fileName2 ReadMode+      contents1 <- hGetContents handle1+      contents2 <- hGetContents handle2+      let list1 = tail $ convertS (parseCSV contents1)+      let list2 = tail $ convertS (parseCSV contents2)+      let tableM = interpolation (constructTab1 list1) (constructTab2 list2)+      writeFile "output.csv" (unlines (convertTable tableM)) +      putStrLn "Interpolated data written to output.csv"+    _ -> do+      progName <- getProgName+      hPutStrLn stderr ("Usage: " ++ progName ++ " <control filename> <gps filename>")+      hPutStrLn stderr ("Example: " ++ progName ++ " \"slam_control.csv\" \"slam_gps.csv\"")+      exitFailure      ++-- parsec (CSV)+csvFile = endBy line eol+line = sepBy cell (char ',')+cell = many (noneOf ",\n")+eol = char '\n'++parseCSV :: String -> Either ParseError [[String]]+parseCSV input = parse csvFile "(unknown)" input++-- convert table to output format (CSV)+convertTable :: MergedT -> [String]+convertTable = map convertCSV++convertCSV :: Merged -> String+convertCSV m = show (time1 control) ++ "," ++ +               show (vel control) ++ "," ++ +               show (steer control) ++ "," ++ +               show (lat sensors) ++ "," ++ +               show (long sensors) ++ "," ++ +               show (orient sensors)+    where control = cm m+          sensors = sens m++-- Perform Interpolation+interpolation :: Table1 -> Table2 -> MergedT+interpolation [] [] = []+interpolation [] _ = error "not enough data to interpolate"+interpolation _ [] = error "not enough data to interpolate"+interpolation (x1:xs) (y1:ys) =+  go (Tr Nothing x1 xs) (Tr Nothing y1 ys)++-- interpolation using current and previous tracking+go :: Tracking ControlMeas -> Tracking Sensor -> MergedT+go (Tr _ _ []) (Tr _ _ _) = []+go (Tr pr1 cur1 rst1) (Tr pr2 cur2 rst2) = +  let t1 = time1 cur1+      t2 = time2 cur2 in+  case compare t1 t2 of+    LT -> M cur1 res : go (Tr (Just cur1) (head rst1) (tail rst1)) (Tr pr2 cur2 rst2)+          where+            S t2p latp longp orientp = maybe (S t1 0 0 0) id pr2+            S t2c latc longc orientc = cur2+            interp x y = ((x-y) / (t2c-t2p)) * (t1-t2p) + y+            res = S t1 (interp latc latp) (interp longc longp) (interp orientc orientp)+    EQ -> M cur1 cur2 : go (Tr (Just cur1) (head rst1) (tail rst1)) (Tr (Just cur2) (head rst2) (tail rst2))+    GT -> M res cur2 : if null rst2 then [] else go (Tr pr1 cur1 rst1) (Tr (Just cur2) (head rst2) (tail rst2))+          where+            CM t1p velp steerp = maybe (CM t2 0 0) id pr1+            CM t1c velc steerc = cur1+            interp x y = ((x-y) / (t1c-t1p)) * (t2-t1p) + y+            res = CM t2 (interp velc velp) (interp steerc steerp) ++-- helper functions+readD :: String -> Double+readD x = read x :: Double++convertS :: Either ParseError a -> a+convertS (Left _) = error "something went wrong in the parsing -- FIXME"+convertS (Right s) = s++constructTab1 :: [[String]] -> Table1+constructTab1 = map read3+  where +    read3 :: [String] -> ControlMeas+    read3 (x : y : z : []) = CM (readD x) (readD y) (readD z)+    read3 _ = error "Table 1 should have exactly 3 entries per row"++constructTab2 :: [[String]] -> Table2+constructTab2 = map read4+  where+    read4 :: [String] -> Sensor+    read4 (x : y : z : w : []) = S (readD x) (readD y) (readD z) (readD w)+    read4 _ = error "Table 2 should have exactly 4 entries per row"
+ Util/Finite.hs view
@@ -0,0 +1,85 @@+module Util.Finite (Finite(..), enumEverything, enumCardinality, suchThat) where++import Data.List (tails)+import Data.Maybe (fromJust)+import Data.Bits (shiftL)+import qualified Data.Set as S+import qualified Data.Map as M++class (Ord a) => Finite a where+    everything :: [a]+    cardinality :: a -> Integer++enumEverything :: (Enum a, Bounded a) => [a]+enumEverything = [minBound..maxBound]++enumCardinality :: (Enum a, Bounded a) => a -> Integer+enumCardinality dummy = succ+                      $ fromIntegral (fromEnum (maxBound `asTypeOf` dummy))+                      - fromIntegral (fromEnum (minBound `asTypeOf` dummy))++instance Finite () where+    everything = enumEverything+    cardinality = enumCardinality++instance Finite Bool where+    everything = enumEverything+    cardinality = enumCardinality++instance Finite Ordering where+    everything = enumEverything+    cardinality = enumCardinality++instance (Finite a) => Finite (Maybe a) where+    everything = Nothing : map Just everything+    cardinality = succ . cardinality . fromJust++instance (Finite a, Finite b) => Finite (Either a b) where+    everything = map Left everything ++ map Right everything+    cardinality x = cardinality l + cardinality r where+        (Left l, Right r) = (x, x)++instance (Finite a, Finite b) => Finite (a, b) where+    everything = [ (a, b) | a <- everything, b <- everything ]+    cardinality ~(a, b) = cardinality a * cardinality b++instance (Finite a, Finite b, Finite c) => Finite (a, b, c) where+    everything = [ (a, b, c) | a <- everything, b <- everything, c <- everything ]+    cardinality ~(a, b, c) = cardinality a * cardinality b * cardinality c++instance (Finite a, Finite b, Finite c, Finite d) => Finite (a, b, c, d) where+    everything = [ (a, b, c, d) | a <- everything, b <- everything, c <- everything, d <- everything ]+    cardinality ~(a, b, c, d) = cardinality a * cardinality b * cardinality c * cardinality d++instance (Finite a, Finite b, Finite c, Finite d, Finite e) => Finite (a, b, c, d, e) where+    everything = [ (a, b, c, d, e) | a <- everything, b <- everything, c <- everything, d <- everything, e <- everything ]+    cardinality ~(a, b, c, d, e) = cardinality a * cardinality b * cardinality c * cardinality d * cardinality e++instance (Finite a) => Finite (S.Set a) where+    everything = loop everything S.empty where+        loop candidates set = set+            : concat [ loop xs (S.insert x set) | x:xs <- tails candidates ]+    cardinality set = shiftL 1 (fromIntegral (cardinality (S.findMin set)))++instance (Finite a, Eq b) => Eq (a -> b) where+    f == g = all (\x -> f x == g x) everything+    f /= g = any (\x -> f x /= g x) everything++instance (Finite a, Ord b) => Ord (a -> b) where+    f `compare` g = map f everything `compare` map g everything+    f <         g = map f everything <         map g everything+    f >         g = map f everything >         map g everything+    f <=        g = map f everything <=        map g everything+    f >=        g = map f everything >=        map g everything++instance (Finite a, Finite b) => Finite (a -> b) where+    everything = [ (M.!) (M.fromDistinctAscList m)+                 | m <- loop everything ] where+        loop []     = [[]]+        loop (a:as) = [ (a,b):rest | b <- everything, rest <- loop as ]+    cardinality f = cardinality y ^ cardinality x where+        (x, y) = (x, f x)++suchThat :: (Finite a) => (a -> Bool) -> S.Set a+suchThat p = S.fromDistinctAscList (filter p everything)+
+ Util/HList.hs view
@@ -0,0 +1,22 @@+{-# LANGUAGE TypeFamilies, DataKinds, TypeOperators #-}+{-# OPTIONS -W #-}++module Util.HList where++class TList (xs :: [*]) where+  data VList (xs :: [*]) :: *+  type Append (xs :: [*]) (ys :: [*]) :: [*]+  append :: VList xs -> VList ys -> VList (Append xs ys)+  vsplit :: VList (Append xs ys) -> (VList xs, VList ys)++instance TList '[] where+  data VList '[] = VNil+  type Append '[] ys = ys+  append VNil ys = ys+  vsplit ys = (VNil, ys)++instance TList xs => TList (x ': xs) where+  data VList (x ': xs) = VCons x (VList xs)+  type Append (x ': xs) ys = x ': Append xs ys+  append (VCons x xs) ys = VCons x (append xs ys)+  vsplit (VCons x zs) = let (xs, ys) = vsplit zs in (VCons x xs, ys)
+ Util/Pretty.hs view
@@ -0,0 +1,74 @@+module Util.Pretty (Pretty(..)) where++import Text.PrettyPrint+import Text.Show.Functions+import Data.Ratio (Ratio, numerator, denominator)+import qualified Data.Map as M+import qualified Data.Set as S+import Util.Finite++class (Show a) => Pretty a where+    pretty :: a -> Doc+    pretty = text . show+    prettyList :: [a] -> Doc+    prettyList = brackets . nest 1 . fsep . punctuate comma . map pretty++instance Pretty Bool+instance Pretty Int+instance Pretty Integer+instance Pretty Float+instance Pretty Double+instance Pretty ()+instance Pretty Ordering+instance Pretty Char where prettyList = text . show++instance (Pretty a, Integral a) => Pretty (Ratio a) where+    pretty r | denom == 1 = prnum+	     | otherwise  = cat [prnum, char '/' <> pretty denom]+	where denom = denominator r+	      prnum = pretty (numerator r)++instance (Pretty a) => Pretty [a] where+    pretty = prettyList++instance (Pretty a) => Pretty (Maybe a) where+    pretty Nothing  = text "Nothing"+    pretty (Just x) = text "Just" <+> pretty x++instance (Pretty a, Pretty b) => Pretty (Either a b) where+    pretty (Left  x) = text "Left"  <+> pretty x+    pretty (Right x) = text "Right" <+> pretty x++instance (Finite a, Pretty a, Pretty b) => Pretty (a -> b)+  where+    pretty f = braces $ nest 1 $ sep $ punctuate comma $+               [ hang (pretty x <> colon) 1 (pretty (f x)) | x <- everything ]++instance (Pretty a, Pretty b) => Pretty (M.Map a b)+  where+    pretty m = braces . nest 1 . sep . punctuate comma+        $ [ hang (pretty k <> colon) 1 (pretty v) | (k,v) <- M.assocs m ]++instance (Pretty a) => Pretty (S.Set a)+  where+    pretty = braces . nest 1 . fsep . punctuate comma . map pretty . S.elems++tuple :: [Doc] -> Doc+tuple = parens . nest 1 . fsep . punctuate comma++{- The Haskell code below is generated by the following Perl program.+@a = 'a'..'z';+$" = ", ";+print <<END foreach 2..5;+instance (@{[map "Pretty $_", @a[0..$_-1]]}) => Pretty (@a[0..$_-1]) where+    pretty (@a[0..$_-1]) = tuple [@{[map "pretty $_", @a[0..$_-1]]}]+END+-}+instance (Pretty a, Pretty b) => Pretty (a, b) where+    pretty (a, b) = tuple [pretty a, pretty b]+instance (Pretty a, Pretty b, Pretty c) => Pretty (a, b, c) where+    pretty (a, b, c) = tuple [pretty a, pretty b, pretty c]+instance (Pretty a, Pretty b, Pretty c, Pretty d) => Pretty (a, b, c, d) where+    pretty (a, b, c, d) = tuple [pretty a, pretty b, pretty c, pretty d]+instance (Pretty a, Pretty b, Pretty c, Pretty d, Pretty e) => Pretty (a, b, c, d, e) where+    pretty (a, b, c, d, e) = tuple [pretty a, pretty b, pretty c, pretty d, pretty e]
+ Visual.hs view
@@ -0,0 +1,43 @@+{-# LANGUAGE OverloadedStrings #-}++module Visual where++import System.IO+import Control.Monad++import Data.Aeson+import Data.List+import qualified Data.Text as T+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.ByteString.Char8 as BS++plot :: Show a => [a] -> String -> IO ()+plot samples filename = do+  h <- openFile filename WriteMode+  hPrint h samples+  hClose h++batchPrint :: Show a => Int -> [a] -> IO ()+batchPrint n l = do+  let batch = take n l+  print batch+  when (length batch == n) $ batchPrint n (drop n l)++viz :: ToJSON a => Int -> [String] -> [[a]] -> IO ()+viz n name samples = viz' n 50 0 name samples++viz' :: ToJSON a => Int -> Int -> Int -> [String] -> [[a]] -> IO ()+viz' n c cur name samples = do+  putStrLn batch+  when (c+cur < n) $+       viz' n c (cur+c) name (drop c samples)+  where+    total = "total_samples" .= n+    current_sample = "current_sample" .= cur+    chunk = object (zipWith (\ name s -> T.pack name .= s)+                            name+                            (transpose $ take c samples))+    batch = B.unpack $ encode+            (object ["rvars" .= chunk,+                     total,+                     current_sample])
+ hakaru.cabal view
@@ -0,0 +1,25 @@+-- Initial hakaru.cabal generated by cabal init.  For further +-- documentation, see http://haskell.org/cabal/users-guide/++name:                hakaru+version:             0.1+synopsis:            A probabilistic programming embedded DSL+-- description:         +homepage:            http://www.indiana.edu/~ppaml+license:             BSD3+license-file:        LICENSE+author:              The Hakaru Team+maintainer:          ppaml@indiana.edu+-- copyright:           +category:            Language+build-type:          Simple+-- extra-source-files:  +cabal-version:       >=1.10++library+  exposed-modules:     Types, Visual, Syntax, Mixture, Language.Hakaru.Symbolic, Sampler, RandomChoice, Util.Csv, Util.Extras, Util.Coda, Util.Finite, Util.Pretty, Util.HList, Util.FileInterpolater, Examples.Tests, Examples.Examples, Language.Hakaru.ImportanceSampler, Language.Hakaru.Metropolis+  -- other-modules:       +  other-extensions:    RankNTypes, BangPatterns, OverloadedStrings, TypeFamilies, ConstraintKinds, GADTs, FlexibleContexts, TypeOperators, DataKinds, NoMonomorphismRestriction, DeriveDataTypeable, ScopedTypeVariables, ExistentialQuantification, StandaloneDeriving+  build-depends:       base >=4.6 && <4.7, aeson >=0.7 && <0.8, text >=1.1 && <1.2, bytestring >=0.10 && <0.11, pretty >=1.1 && <1.2, logfloat >=0.12 && <0.13, containers >=0.5 && <0.6, random >=1.0 && <1.1, math-functions >=0.1 && <0.2, vector >=0.10 && <0.11, cassava >=0.4 && <0.5, zlib >=0.5 && <0.6, statistics >=0.11 && <0.12, hmatrix >=0.16 && <0.17, parsec >=3.1 && <3.2+  -- hs-source-dirs:      +  default-language:    Haskell2010