diff --git a/Examples/Examples.hs b/Examples/Examples.hs new file mode 100644 --- /dev/null +++ b/Examples/Examples.hs @@ -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 diff --git a/Examples/Tests.hs b/Examples/Tests.hs new file mode 100644 --- /dev/null +++ b/Examples/Tests.hs @@ -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))] diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Language/Hakaru/ImportanceSampler.hs b/Language/Hakaru/ImportanceSampler.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/ImportanceSampler.hs @@ -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)) + diff --git a/Language/Hakaru/Metropolis.hs b/Language/Hakaru/Metropolis.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Metropolis.hs @@ -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 + diff --git a/Language/Hakaru/Symbolic.hs b/Language/Hakaru/Symbolic.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Symbolic.hs @@ -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 diff --git a/Mixture.hs b/Mixture.hs new file mode 100644 --- /dev/null +++ b/Mixture.hs @@ -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) diff --git a/RandomChoice.hs b/RandomChoice.hs new file mode 100644 --- /dev/null +++ b/RandomChoice.hs @@ -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" diff --git a/Sampler.hs b/Sampler.hs new file mode 100644 --- /dev/null +++ b/Sampler.hs @@ -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)) diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/Syntax.hs b/Syntax.hs new file mode 100644 --- /dev/null +++ b/Syntax.hs @@ -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 diff --git a/Types.hs b/Types.hs new file mode 100644 --- /dev/null +++ b/Types.hs @@ -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) diff --git a/Util/Coda.hs b/Util/Coda.hs new file mode 100644 --- /dev/null +++ b/Util/Coda.hs @@ -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 diff --git a/Util/Csv.hs b/Util/Csv.hs new file mode 100644 --- /dev/null +++ b/Util/Csv.hs @@ -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 diff --git a/Util/Extras.hs b/Util/Extras.hs new file mode 100644 --- /dev/null +++ b/Util/Extras.hs @@ -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 diff --git a/Util/FileInterpolater.hs b/Util/FileInterpolater.hs new file mode 100644 --- /dev/null +++ b/Util/FileInterpolater.hs @@ -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 ++ " ") + 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" diff --git a/Util/Finite.hs b/Util/Finite.hs new file mode 100644 --- /dev/null +++ b/Util/Finite.hs @@ -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) + diff --git a/Util/HList.hs b/Util/HList.hs new file mode 100644 --- /dev/null +++ b/Util/HList.hs @@ -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) diff --git a/Util/Pretty.hs b/Util/Pretty.hs new file mode 100644 --- /dev/null +++ b/Util/Pretty.hs @@ -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 < 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] diff --git a/Visual.hs b/Visual.hs new file mode 100644 --- /dev/null +++ b/Visual.hs @@ -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]) diff --git a/hakaru.cabal b/hakaru.cabal new file mode 100644 --- /dev/null +++ b/hakaru.cabal @@ -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