diff --git a/Bench/Bench.hs b/Bench/Bench.hs new file mode 100644 --- /dev/null +++ b/Bench/Bench.hs @@ -0,0 +1,20 @@ +import Criterion.Main + +import Control.Monad + +import Language.Hakaru.Distribution +import qualified Language.Hakaru.Metropolis as MH +import qualified Language.Hakaru.ImportanceSampler as IS + +giveLast :: IO [(a,b)] -> IO (a,b) +giveLast samples = do s <- samples + return $ last (take 100 s) + +main = defaultMain [ + bcompare [ + bench "is normal 10" $ whnfIO $ giveLast (IS.sample (replicateM 10 (IS.unconditioned (normal 0 10))) []) + , bench "is normal 20" $ whnfIO $ giveLast (IS.sample (replicateM 20 (IS.unconditioned (normal 0 10))) []) + , bench "mh normal 10" $ whnfIO $ giveLast (MH.sample (replicateM 10 (MH.unconditioned (normal 0 10))) []) + , bench "mh normal 20" $ whnfIO $ giveLast (MH.sample (replicateM 20 (MH.unconditioned (normal 0 10))) []) + ] + ] diff --git a/Examples/Tests.hs b/Examples/Tests.hs deleted file mode 100644 --- a/Examples/Tests.hs +++ /dev/null @@ -1,106 +0,0 @@ -{-# LANGUAGE RankNTypes, NoMonomorphismRestriction, BangPatterns #-} - -module Examples.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 --- a/LICENSE +++ 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/Arrow.hs b/Language/Hakaru/Arrow.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Arrow.hs @@ -0,0 +1,6 @@ +{-# LANGUAGE TypeOperators #-} +module Language.Hakaru.Arrow where + +import Language.Hakaru.Types (Dist) + +type a ~~> b = a -> Dist b diff --git a/Language/Hakaru/Distribution.hs b/Language/Hakaru/Distribution.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Distribution.hs @@ -0,0 +1,228 @@ +{-# LANGUAGE RankNTypes, BangPatterns, GADTs #-} +{-# OPTIONS -Wall #-} + +module Language.Hakaru.Distribution where + +import System.Random +import Language.Hakaru.Mixture +import Language.Hakaru.Types +import Data.Ix +import Data.Maybe (fromMaybe) +import Data.List (findIndex, foldl') +import Numeric.SpecFunctions +import qualified Data.Map.Strict as M +import qualified Data.Number.LogFloat as LF + +mapFst :: (t -> s) -> (t, u) -> (s, u) +mapFst f (a,b) = (f a, b) + +dirac :: (Eq a) => a -> Dist a +dirac theta = Dist {logDensity = (\ (Discrete x) -> if x == theta then 0 else log 0), + distSample = (\ g -> (Discrete theta,g))} + +bern :: Double -> Dist Bool +bern p = Dist {logDensity = (\ (Discrete x) -> log (if x then p else 1 - p)), + distSample = (\ g -> case randomR (0, 1) g of + (t, g') -> (Discrete $ t <= p, g'))} + +uniform :: Double -> Double -> Dist Double +uniform lo hi = + let uniformLogDensity lo' hi' x | lo' <= x && x <= hi' = log (recip (hi' - lo')) + uniformLogDensity _ _ _ = log 0 + in Dist {logDensity = (\ (Lebesgue x) -> uniformLogDensity lo hi x), + distSample = (\ g -> mapFst Lebesgue $ randomR (lo, hi) g)} + +uniformD :: (Ix a, Random a) => a -> a -> Dist a +uniformD lo hi = + let uniformLogDensity lo' hi' x | lo' <= x && x <= hi' = log density + uniformLogDensity _ _ _ = log 0 + density = recip (fromInteger (toInteger (rangeSize (lo,hi)))) + in Dist {logDensity = (\ (Discrete x) -> uniformLogDensity lo hi x), + distSample = (\ g -> mapFst Discrete $ randomR (lo, hi) g)} + +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 :: Floating a => a -> a -> a -> a +normalLogDensity mu sd x = (-tau * square (x - mu) + + log (tau / pi / 2)) / 2 + where square y = y * y + tau = 1 / square sd + +normal :: Double -> Double -> Dist Double +normal mu sd = Dist {logDensity = normalLogDensity mu sd . fromLebesgue, + distSample = mapFst Lebesgue . normal_rng mu sd} + +categoricalLogDensity :: (Eq b, Floating a) => [(b, a)] -> b -> a +categoricalLogDensity list x = log $ fromMaybe 0 (lookup x list) +categoricalSample :: (Num b, Ord b, RandomGen g, Random b) => + [(t,b)] -> g -> (t, g) +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 + +categorical :: Eq a => [(a,Double)] -> Dist a +categorical list = Dist {logDensity = categoricalLogDensity list . fromDiscrete, + distSample = mapFst Discrete . categoricalSample list} + +lnFact :: Integer -> Double +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 -> (Integer, 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 -> (Integer,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 () of + () | us >= 0.07 && v <= vr -> (k, g2) + () | k < 0 -> make_poisson g2 + () | us <= 0.013 && v > us -> make_poisson g2 + () | accept_region us v k -> (k, g2) + _ -> make_poisson g2 + + accept_region :: Double -> Double -> Integer -> Bool + accept_region us v k = log (v * arep / (a/(us*us)+b)) <= + -lambda + (fromIntegral k)*lnlam - lnFact k + +poisson :: Double -> Dist Integer +poisson l = + let poissonLogDensity l' x | l' > 0 && x> 0 = (fromIntegral x)*(log l') - lnFact x - l' + poissonLogDensity l' x | x==0 = -l' + poissonLogDensity _ _ = log 0 + in Dist {logDensity = poissonLogDensity l . fromDiscrete, + distSample = mapFst Discrete . poisson_rng l} + +-- 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 _ _ | shape <= 0.0 = error "gamma: got a negative shape paramater" +gamma_rng _ scl _ | scl <= 0.0 = error "gamma: got a negative scale paramater" +gamma_rng shape scl g | shape < 1.0 = (gvar2, g2) + where (gvar1, g1) = gamma_rng (shape + 1) scl g + (w, g2) = randomR (0,1) g1 + gvar2 = scl * gvar1 * (w ** recip shape) +gamma_rng shape scl g = + let d = shape - 1/3 + c = recip $ sqrt $ 9*d + -- Algorithm recommends inlining normal generator + n = normal_rng 1 c + (v, g2) = until (\y -> fst y > 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 -> (scl*d*v3, g3) + False -> gamma_rng shape scl g3 + +gammaLogDensity :: Double -> Double -> Double -> Double +gammaLogDensity shape scl x | x>= 0 && shape > 0 && scl > 0 = + scl * log shape - scl * x + (shape - 1) * log x - logGamma shape +gammaLogDensity _ _ _ = log 0 + +gamma :: Double -> Double -> Dist Double +gamma shape scl = Dist {logDensity = gammaLogDensity shape scl . fromLebesgue, + distSample = mapFst Lebesgue . gamma_rng shape scl} + +beta_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g) +beta_rng a b g | a <= 1.0 && b <= 1.0 = + let (u, g1) = randomR (0.0, 1.0) g + (v, g2) = randomR (0.0, 1.0) g1 + x = u ** (recip a) + y = v ** (recip b) + in case (x+y) <= 1.0 of + True -> (x / (x + y), g2) + False -> beta_rng a b g2 +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 :: Double -> Double -> Double -> Double +betaLogDensity _ _ x | x < 0 || x > 1 = error "beta: value must be between 0 and 1" +betaLogDensity a b _ | a <= 0 || b <= 0 = error "beta: parameters must be positve" +betaLogDensity a b x = (logGamma (a + b) + - logGamma a + - logGamma b + + (a - 1) * log x + + (b - 1) * log (1 - x)) + +beta :: Double -> Double -> Dist Double +beta a b = Dist {logDensity = betaLogDensity a b . fromLebesgue, + distSample = mapFst Lebesgue . beta_rng a b} + +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 :: Floating a => a -> a -> a -> a +laplaceLogDensity mu sd x = - log (2 * sd) - abs (x - mu) / sd + +laplace :: Double -> Double -> Dist Double +laplace mu sd = Dist {logDensity = laplaceLogDensity mu sd . fromLebesgue, + distSample = mapFst Lebesgue . laplace_rng 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 (b, total, h) = (map (/ total) b, h) + +dirichletLogDensity :: [Double] -> [Double] -> Double +dirichletLogDensity a x | all (> 0) x = sum' (zipWith logTerm a x) + logGamma (sum a) + where sum' = foldl' (+) 0 + logTerm b y = (b-1) * log y - logGamma b +dirichletLogDensity _ _ = error "dirichlet: all values must be between 0 and 1" diff --git a/Language/Hakaru/ImportanceSampler.hs b/Language/Hakaru/ImportanceSampler.hs --- a/Language/Hakaru/ImportanceSampler.hs +++ b/Language/Hakaru/ImportanceSampler.hs @@ -9,111 +9,18 @@ -- 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 Language.Hakaru.Types +import Language.Hakaru.Mixture (Prob, empty, point, Mixture(..)) +import Language.Hakaru.Sampler (Sampler, deterministic, smap, sbind) import System.Random import Data.Monoid -import Data.Ix import Data.Dynamic -import Data.List -import Control.Monad +import System.IO.Unsafe 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]) } @@ -121,53 +28,56 @@ bind measure continuation = Measure (\conds -> sbind (unMeasure measure conds) - (\(a,conds) -> unMeasure (continuation a) conds)) + (\(a,cds) -> unMeasure (continuation a) cds)) 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)) +updateMixture :: Typeable a => Cond -> Dist a -> Sampler a +updateMixture (Just cond) dist = + case fromDynamic cond of + Just y -> deterministic (point (fromDensity y) density) + where density = LF.logToLogFloat $ logDensity dist y + Nothing -> error "did not get data from dynamic source" +updateMixture Nothing dist = \g0 -> let (e, g) = distSample dist g0 + in (point (fromDensity e) 1, g) + +conditioned, unconditioned :: Typeable a => Dist a -> Measure a +conditioned dist = Measure (\(cond:conds) -> smap (\a->(a,conds)) + (updateMixture cond dist)) +unconditioned dist = Measure (\ conds -> smap (\a->(a,conds)) + (updateMixture Nothing dist)) + 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 +condition :: Eq b => Measure (a, b) -> b -> Measure a +condition m b' = + Measure (\ conds -> + sbind (unMeasure m conds) + (\ ((a,b), cds) -> + deterministic (if b==b' then point (a,cds) 1 else empty))) -- 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 +empiricalMeasure :: (Ord a) => Int -> Measure a -> [Cond] -> IO (Mixture a) +empiricalMeasure !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) + go k m = once >>= \result -> go (k - 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 +sample :: (Ord a, Show a) => Measure a -> [Cond] -> IO [(a, Prob)] +sample measure conds = do + u <- once + let x = mixToTuple (finish u) + xs <- unsafeInterleaveIO $ sample measure conds + return (x : xs) + where once = getStdRandom (unMeasure measure conds) + mixToTuple = head . M.toList . unMixture logit :: Floating a => a -> a logit !x = 1 / (1 + exp (- x)) diff --git a/Language/Hakaru/Lambda.hs b/Language/Hakaru/Lambda.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Lambda.hs @@ -0,0 +1,22 @@ +-- The lambda-calculus part of the language, which can be shared +module Language.Hakaru.Lambda(lit, dbl, lam, app, fix, ifThenElse) where + +lit :: (Eq a) => a -> a +lit = id + +-- raw lit is a pain to use. These are nicer +dbl :: Double -> Double +dbl = lit + +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 diff --git a/Language/Hakaru/Metropolis.hs b/Language/Hakaru/Metropolis.hs --- a/Language/Hakaru/Metropolis.hs +++ b/Language/Hakaru/Metropolis.hs @@ -1,21 +1,17 @@ {-# LANGUAGE RankNTypes, NoMonomorphismRestriction, BangPatterns, DeriveDataTypeable, GADTs, ScopedTypeVariables, ExistentialQuantification, StandaloneDeriving #-} +{-# OPTIONS -Wall #-} 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 +import Language.Hakaru.Types {- @@ -28,254 +24,160 @@ -} 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) +-- and what does XRP stand for? +data XRP where + XRP :: Typeable e => (Density e, Dist e) -> XRP -unXRP :: Typeable a => XRP -> Maybe (a, Dist a) +unXRP :: Typeable a => XRP -> Maybe (Density 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 LL = LogLikelihood 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) +data DBEntry = DBEntry { + xrp :: XRP, + llhd :: LL, + vis :: Visited, + observed :: Observed } +type Database = M.Map Name DBEntry --- n is structural_name --- d is database --- ll is likelihood of expression --- conds is the observed data --- g is the random seed +data SamplerState g where + S :: { ldb :: Database, -- ldb = local database + -- (total likelihood, total likelihood of XRPs newly introduced) + llh2 :: {-# UNPACK #-} !(LL, LL), + cnds :: [Cond], -- conditions left to process + seed :: g } -> SamplerState g +type Sampler a = forall g. (RandomGen g) => SamplerState g -> (a, SamplerState g) -lit :: (Eq a, Typeable a) => a -> a -lit = id +sreturn :: a -> Sampler a +sreturn x s = (x, s) +sbind :: Sampler a -> (a -> Sampler b) -> Sampler b +sbind s k = \ st -> let (v, s') = s st in k v s' + +smap :: (a -> b) -> Sampler a -> Sampler b +smap f s = sbind s (\a -> sreturn (f a)) + +newtype Measure a = Measure {unMeasure :: Name -> Sampler a } + deriving (Typeable) + return_ :: a -> Measure a -return_ x = Measure (\ (n, d, l, conds, g) -> (x, d, l, conds, g)) +return_ x = Measure $ \ _ -> sreturn x -makeXRP :: (Typeable a, RandomGen g) => Cond -> Dist a - -> Name -> Database -> g - -> (a, Database, Likelihood, Likelihood, g) -makeXRP obs dist' n db g = +updateXRP :: Typeable a => Name -> Cond -> Dist a -> Sampler a +updateXRP n obs dist' s@(S {ldb = db, seed = g}) = case M.lookup n db of - Just (xd, lb, b, ob) -> + Just (DBEntry xd lb _ ob) -> let Just (xb, dist) = unXRP xd - (x,l) = case obs of - Just xd -> - let Just x = fromDynamic xd - in (x, logDensity dist x) + (x,_) = case obs of + Just yd -> + let Just y = fromDynamic yd + in (y, logDensity dist y) Nothing -> (xb, lb) l' = logDensity dist' x - d1 = M.insert n (XRP (x,dist), - l', - True, - ob) db - in (x, d1, l', 0, g) + d1 = M.insert n (DBEntry (XRP (x,dist)) l' True ob) db + in (fromDensity x, + s {ldb = d1, + llh2 = updateLogLikelihood l' 0 s, + seed = g}) Nothing -> - let (xnew, l, g1) = case obs of + let (xnew2, l, g2) = 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 + let (xnew, g1) = distSample dist' g + in (xnew, logDensity dist' xnew, g1) + d1 = M.insert n (DBEntry (XRP (xnew2, dist')) l True (isJust obs)) db + in (fromDensity xnew2, + s {ldb = d1, + llh2 = updateLogLikelihood l l s, + seed = g2}) -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 +updateLogLikelihood :: RandomGen g => + LL -> LL -> SamplerState g -> + (LL, LL) +updateLogLikelihood llTotal llFresh s = + let (l,lf) = llh2 s in (llTotal+l, llFresh+lf) -factor :: Likelihood -> Measure () -factor l = Measure $ \(n, d, (llTotal, llFresh), conds, g) -> - ((), d, (llTotal + l, llFresh), conds, g) +factor :: LL -> Measure () +factor l = Measure $ \ _ -> \ s -> + let (llTotal, llFresh) = llh2 s + in ((), s {llh2 = (llTotal + l, llFresh)}) -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) +condition :: Eq b => Measure (a, b) -> b -> Measure a +condition (Measure m) b' = Measure $ \ n -> + let comp a b s | a /= b = s {llh2 = (log 0, 0)} + comp _ _ s = s + in sbind (m n) (\ (a, b) s -> (a, comp b b' s)) 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) +bind (Measure m) cont = Measure $ \ n -> + sbind (m (0:n)) (\ a -> unMeasure (cont a) (1:n)) -conditioned :: (Cond -> Measure a) -> Measure a -conditioned f = Measure $ \ (n,d,ll,cond:conds,g) -> - unMeasure (f cond) (n, d, ll, conds, g) +conditioned :: Typeable a => Dist a -> Measure a +conditioned dist = Measure $ \ n -> + \s@(S {cnds = cond:conds }) -> + updateXRP n cond dist s{cnds = conds} -unconditioned :: (Cond -> Measure a) -> Measure a -unconditioned f = f Nothing +unconditioned :: Typeable a => Dist a -> Measure a +unconditioned dist = Measure $ \ n -> + updateXRP n Nothing dist 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 + (>>=) = bind -run :: Measure a -> [Cond] -> IO (a, Database, Likelihood) -run (Measure prog) conds = do +run :: Measure a -> [Cond] -> IO (a, Database, LL) +run (Measure prog) cds = do g <- getStdGen - let (v, d, ll, conds1, g') = - prog ([0], M.empty, (0,0), conds, g) + let (v, S d ll [] _) = (prog [0]) (S M.empty (0,0) cds 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 + -> (a, Database, LL, LL, LL, g) +traceUpdate (Measure prog) d cds g = do + -- let d1 = M.map (\ (x, l, _, ob) -> (x, l, False, ob)) d + let d1 = M.map (\ s -> s { vis = False }) d + let (v, S d2 (llTotal, llFresh) [] g1) = (prog [0]) (S d1 (0,0) cds g) + let (d3, stale_d) = M.partition vis d2 + let llStale = M.foldl' (\ llStale' s -> llStale' + llhd s) 0 stale_d (v, d3, llTotal, llFresh, llStale, g1) initialStep :: Measure a -> [Cond] -> IO (a, Database, - Likelihood, Likelihood, Likelihood, StdGen) -initialStep prog conds = do + LL, LL, LL, StdGen) +initialStep prog cds = do g <- getStdGen - return $ traceUpdate prog M.empty conds g + return $ traceUpdate prog M.empty cds 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 +resample :: RandomGen g => Name -> Database -> Observed -> XRP -> g -> + (Database, LL, LL, LL, g) +resample name db ob (XRP (x, dist)) g = + let (x', g1) = distSample dist g + fwd = logDensity dist x' + rvs = logDensity dist x + l' = fwd + newEntry = DBEntry (XRP (x', dist)) l' True ob + db' = M.insert name newEntry db + in (db', l', fwd, rvs, g1) transition :: (Typeable a, RandomGen g) => Measure a -> [Cond] - -> a -> Database -> Likelihood -> g -> [a] -transition prog conds v db ll g = + -> a -> Database -> LL -> g -> [a] +transition prog cds v db ll g = let dbSize = M.size db -- choose an unconditioned choice - (condDb, uncondDb) = M.partition (\ (_, _, _, ob) -> ob) db + (_, uncondDb) = M.partition observed 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 + (name, (DBEntry xd _ _ ob)) = M.elemAt choice uncondDb + (db', _, fwd, rvs, g2) = resample name db ob xd g1 + (v', db2, llTotal, llFresh, llStale, g3) = traceUpdate prog db' cds g2 a = llTotal - ll + rvs - fwd + log (fromIntegral dbSize) - log (fromIntegral $ M.size db2) @@ -283,12 +185,16 @@ (u, g4) = randomR (0 :: Double, 1) g3 in if (log u < a) then - v' : (transition prog conds v' db2 llTotal g4) + v' : (transition prog cds v' db2 llTotal g4) else - v : (transition prog conds v db ll g4) + v : (transition prog cds 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 +mcmc prog cds = do + (v, d, llTotal, _, _, g) <- initialStep prog cds + return $ transition prog cds v d llTotal g +sample :: Typeable a => Measure a -> [Cond] -> IO [(a, Double)] +sample prog cds = do + (v, d, llTotal, _, _, g) <- initialStep prog cds + return $ map (\ x -> (x,1)) (transition prog cds v d llTotal g) diff --git a/Language/Hakaru/Mixture.hs b/Language/Hakaru/Mixture.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Mixture.hs @@ -0,0 +1,61 @@ +{-# LANGUAGE RankNTypes, BangPatterns #-} +{-# OPTIONS -W #-} + +module Language.Hakaru.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/Language/Hakaru/Sampler.hs b/Language/Hakaru/Sampler.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Sampler.hs @@ -0,0 +1,26 @@ +{-# LANGUAGE RankNTypes #-} +{-# OPTIONS -W #-} + +module Language.Hakaru.Sampler (Sampler, deterministic, sbind, smap) where + +import Language.Hakaru.Mixture (Mixture, mnull, empty, scale, point) +import Language.Hakaru.Distribution (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/Language/Hakaru/Symbolic.hs b/Language/Hakaru/Symbolic.hs --- a/Language/Hakaru/Symbolic.hs +++ b/Language/Hakaru/Symbolic.hs @@ -1,78 +1,185 @@ -{-# LANGUAGE GADTs, TypeFamilies #-} -{-# OPTIONS -W #-} +{-# LANGUAGE GADTs, TypeFamilies, ScopedTypeVariables, FlexibleContexts #-} +{-# OPTIONS -Wall #-} module Language.Hakaru.Symbolic where +import Prelude hiding (Real) + +data Real data Prob data Measure a data Dist a +data Exact --- 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) +class IntComp repr where + int :: Integer -> repr Integer + +class BoolComp repr where + bool :: Bool -> repr Bool + +class RealComp repr where + real :: Rational -> repr Real + exp :: repr Real -> repr Real -> repr Real + sqrt, cos, sin :: repr Real -> repr Real + +-- polymorphic operations and integer powering +-- probably should restrict 'a' to be some kind of numeric type? +class SymbComp repr where + add, minus, mul :: repr a -> repr a -> repr a + pow :: repr Real -> repr Integer -> repr Real + scale :: repr Integer -> repr Real -> repr Real + +class MeasMonad repr where + bind :: repr (Measure a) -> (repr a -> repr (Measure b)) + -> repr (Measure b) + ret :: repr a -> repr (Measure a) + +class Distrib repr where + uniform, normal :: repr Real -> repr Real -> repr (Dist Real) + uniformD :: repr Integer -> repr Integer -> repr (Dist Integer) + +class Conditioning repr where conditioned, unconditioned :: repr (Dist a) -> repr (Measure a) -- Printer (to Maple) +data Pos = Front | Back type VarCounter = Int -newtype Maple a = Maple { unMaple :: Bool -> VarCounter -> String } +newtype Maple a = Maple { unMaple :: Pos -> 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 ++ ")" +type family MPL a +type instance MPL Real = Rational +type instance MPL Integer = Integer +type instance MPL Bool = Bool -view e = unMaple e True 0 +-- if it weren't for the constraints, we could/should use Applicative +pure :: Show (MPL a) => MPL a -> Maple a +pure x = Maple $ \_ _ -> show x +liftA1 :: (String -> String) -> Maple a -> Maple a +liftA1 pr x = Maple $ \f h -> pr (unMaple x f h) + +liftA2 :: (String -> String -> String) -> Maple a -> Maple a -> Maple a +liftA2 pr e1 e2 = Maple $ \f h -> pr (unMaple e1 f h) (unMaple e2 f h) + +-- variant for ret +liftA1M :: (String -> String) -> Maple a -> Maple (Measure a) +liftA1M pr x = Maple $ \f h -> pr (unMaple x f h) + +-- variant for pow +liftA2aba :: (String -> String -> String) -> Maple a -> Maple b -> Maple a +liftA2aba pr e1 e2 = Maple $ \f h -> pr (unMaple e1 f h) (unMaple e2 f h) + +-- variant for scale +liftA2baa :: (String -> String -> String) -> Maple b -> Maple a -> Maple a +liftA2baa pr e1 e2 = Maple $ \f h -> pr (unMaple e1 f h) (unMaple e2 f h) + +mkPr :: String -> (String -> String) +mkPr s t = s ++ "(" ++ t ++ ")" + +d2m :: Maple (Dist a) -> Maple (Measure a) +d2m e = Maple $ unMaple e + +infixPr :: String -> (String -> String -> String) +infixPr s a b = a ++ s ++ b + +-- This is quite scary. Probably a mistake. +reify :: forall a. Read a => Pos -> VarCounter -> Maple a -> a +reify f h e = (read (unMaple e f h) :: a) + +name :: String -> VarCounter -> String +name s h = s ++ show h + +var :: String -> VarCounter -> Maple a +var s h = Maple $ \_ _ -> name s h + +binder :: (String -> String -> Maybe String) -> + (String -> String -> VarCounter -> Maybe String) -> + String -> + Maple a -> Maple a -> Maple (Dist a) +binder pr1 pr2 oper e1 e2 = Maple $ pr_ + where + pr_ Front h = let x1 = unMaple e1 Front h + x2 = unMaple e2 Front h in + case (pr1 x1 x2, pr2 x1 x2 h) of + (Just z, Just w) -> z ++ " * " ++ oper ++ " (" ++ w ++ " * " + (Nothing, Just w) -> oper ++ " (" ++ w ++ " * " + (Just z, Nothing) -> z ++ " * " ++ oper ++ " (" + (Nothing, Nothing) -> oper ++ " (" + pr_ Back h = ", " ++ (name "x" h) ++ "=" ++ unMaple e1 Back h ++ + ".." ++ unMaple e2 Back h ++ ")" + +instance RealComp Maple where + -- serious problem here: all exact numbers will be printed as + -- floats, which will really hamper the use of Maple in any + -- serious way. This needs a rethink. + real = pure + exp = liftA2 $ infixPr "^" + sqrt = liftA1 $ mkPr "sqrt" + cos = liftA1 $ mkPr "cos" + sin = liftA1 $ mkPr "sin" + +instance SymbComp Maple where + add = liftA2 $ infixPr "+" + minus = liftA2 $ infixPr "-" + mul = liftA2 $ infixPr "*" + pow = liftA2aba $ infixPr "^" + scale = liftA2baa $ infixPr "*" + +instance IntComp Maple where + int = pure + +instance BoolComp Maple where + bool = pure + +instance MeasMonad Maple where + ret = liftA1M $ mkPr "g" + bind m c = Maple $ \f h -> unMaple m Front h ++ + unMaple (c $ var "x" h) f (succ h) ++ + unMaple m Back h + +instance Distrib Maple where + uniform = binder (\e1 e2 -> Just $ + show (1/((read e2 :: Rational) - (read e1 :: Rational)))) + (\_ _ _ -> Nothing) "Int" + uniformD = binder (\e1 e2 -> + let d = (read e2 :: Integer) - (read e1) in + if d == 1 then Nothing + else Just $ "(1 / " ++ show d ++ ")") + (\_ _ _ -> Nothing) "Sum" + normal = binder (\_ _ -> Nothing) + (\e1 e2 h -> Just $ "PDF (Normal (" ++ e1 ++ ", " ++ e2 ++ "), " ++ (name "x" h) ++ ")") + "Int" + +instance Conditioning Maple where + unconditioned = d2m + conditioned = d2m + +view :: Maple a -> String +view e = unMaple e Front 0 + +lift :: Maple Integer -> Maple Real +lift x = Maple $ \f h -> unMaple x f h + -- TEST CASES -exp1 = unconditioned (uniformC (real 1) (real 3)) `bind` \s -> - ret s +exp1, exp2, exp3, exp4 :: Maple (Measure Real) +exp1 = unconditioned (uniform (real 1) (real 3)) `bind` ret -- 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 +exp2 = unconditioned (uniformD (int 1) (int 3)) `bind` \s -> + unconditioned (uniform (real (-1)) (real 1)) `bind` \x -> + let y = s `scale` x in ret y -- Borel's Paradox -exp3 = unconditioned (uniformD (real 1) (real 2)) `bind` \s -> - unconditioned (uniformC (real (-1)) (real 1)) `bind` \x -> +exp3 = unconditioned (uniformD (int 1) (int 2)) `bind` \s -> + unconditioned (uniform (real (-1)) (real 1)) `bind` \x -> let y = (Language.Hakaru.Symbolic.sqrt ((real 1 ) `minus` - (Language.Hakaru.Symbolic.exp s (real 2)))) `mul` - (Language.Hakaru.Symbolic.sin x) in ret y + (Language.Hakaru.Symbolic.pow (lift s) (int 2)))) `mul` + (Language.Hakaru.Symbolic.sin x) in ret y +exp4 = unconditioned (normal (real 1) (real 4)) `bind` ret + +test, test2, test3, test4 :: String test = view exp1 test2 = view exp2 test3 = view exp3 +test4 = view exp4 diff --git a/Language/Hakaru/Syntax.hs b/Language/Hakaru/Syntax.hs deleted file mode 100644 --- a/Language/Hakaru/Syntax.hs +++ /dev/null @@ -1,185 +0,0 @@ -{-# LANGUAGE TypeFamilies, ConstraintKinds, GADTs, FlexibleContexts #-} - -module Language.Hakaru.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 Language.Hakaru.ImportanceSampler as IS - --- The Metropolis-Hastings semantics - -import qualified Language.Hakaru.Metropolis 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/Language/Hakaru/Types.hs b/Language/Hakaru/Types.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Types.hs @@ -0,0 +1,29 @@ +{-# LANGUAGE RankNTypes, BangPatterns, DeriveDataTypeable, StandaloneDeriving #-} +{-# OPTIONS -W #-} + +module Language.Hakaru.Types where + +import Data.Dynamic +import System.Random + +-- Basic types for conditioning and conditioned sampler +data Density a = Lebesgue !a | Discrete !a deriving Typeable +type Cond = Maybe Dynamic + +fromDiscrete :: Density t -> t +fromDiscrete (Discrete a) = a +fromDiscrete _ = error "got a non-discrete sampler" + +fromLebesgue :: Density t -> t +fromLebesgue (Lebesgue a) = a +fromLebesgue _ = error "got a discrete sampler" + +fromDensity :: Density t -> t +fromDensity (Discrete a) = a +fromDensity (Lebesgue a) = a + +type LogLikelihood = Double +data Dist a = Dist {logDensity :: Density a -> LogLikelihood, + distSample :: forall g. + RandomGen g => g -> (Density a, g)} +deriving instance Typeable1 Dist diff --git a/Language/Hakaru/Util/Coda.hs b/Language/Hakaru/Util/Coda.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Util/Coda.hs @@ -0,0 +1,12 @@ +module Language.Hakaru.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/Language/Hakaru/Util/Csv.hs b/Language/Hakaru/Util/Csv.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Util/Csv.hs @@ -0,0 +1,40 @@ +{-# LANGUAGE TypeOperators #-} + +module Language.Hakaru.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/Language/Hakaru/Util/Extras.hs b/Language/Hakaru/Util/Extras.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Util/Extras.hs @@ -0,0 +1,78 @@ +{-| + 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 Language.Hakaru.Util.Extras where + +import qualified Data.Sequence as S +import System.Random +import Data.Maybe +import qualified Data.Foldable as F + +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 diff --git a/Language/Hakaru/Util/Visual.hs b/Language/Hakaru/Util/Visual.hs new file mode 100644 --- /dev/null +++ b/Language/Hakaru/Util/Visual.hs @@ -0,0 +1,43 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Language.Hakaru.Util.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/Mixture.hs b/Mixture.hs deleted file mode 100644 --- a/Mixture.hs +++ /dev/null @@ -1,61 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/RandomChoice.hs +++ /dev/null @@ -1,145 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/Sampler.hs +++ /dev/null @@ -1,26 +0,0 @@ -{-# 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/Tests/Distribution.hs b/Tests/Distribution.hs new file mode 100644 --- /dev/null +++ b/Tests/Distribution.hs @@ -0,0 +1,106 @@ +{-# LANGUAGE BangPatterns #-} + +module Tests.Distribution where + +import Control.Monad +import qualified System.Random.MWC as MWC + +import Language.Hakaru.Types +import Language.Hakaru.Util.Coda +import Language.Hakaru.Distribution hiding (choose) + +import Test.QuickCheck +import Test.QuickCheck.Monadic as QM +import Test.Framework.Providers.QuickCheck2 (testProperty) + +fromDiscreteToNum = fromIntegral . fromEnum . fromDiscrete +sq x = x * x + +almostEqual :: (Fractional a, Ord a) => a -> a -> a -> Bool +almostEqual tol x y | abs (x - y) < tol = True +almostEqual tol x y = (abs $ (x - y) / (x + y)) < tol + +quickArg :: IO () +quickArg = quickCheckWith stdArgs {maxSuccess = 2000} (\ x -> almostEqual tol x x) + where tol :: Double + tol = 1e-5 + +qtest = [testProperty "checking beta" $ QM.monadicIO betaTest, + testProperty "checking bern" $ QM.monadicIO bernTest, + testProperty "checking gamma" $ QM.monadicIO gammaTest, + testProperty "checking normal" $ QM.monadicIO normalTest, + testProperty "checking laplace" $ QM.monadicIO laplaceTest, + testProperty "checking poisson" $ QM.monadicIO poissonTest] + +betaTest = do + Positive a <- QM.pick arbitrary + Positive b <- QM.pick arbitrary + g <- QM.run $ MWC.create + samples <- QM.run $ replicateM 1000 $ distSample (beta a b) g + let (mean, variance) = meanVariance (map fromLebesgue samples) + QM.assert $ (almostEqual tol mean (mu a b)) && + (almostEqual tol variance (var a b)) + where tol = 1e-1 + mu a b = a / (a + b) + var a b = a*b / ((sq $ a + b) * (a + b + 1)) + +bernTest = do + p <- QM.pick $ choose (0, 1) + g <- QM.run $ MWC.create + samples <- QM.run $ replicateM 1000 $ distSample (bern p) g + let (mean, variance) = meanVariance (map fromDiscreteToNum samples) + QM.assert $ (almostEqual tol mean (mu p)) && + (almostEqual tol variance (var p)) + where tol = 1e-1 + mu p = p + var p = p*(1-p) + +poissonTest = do + lam <- QM.pick $ choose (1, 10) + g <- QM.run $ MWC.create + samples <- QM.run $ replicateM 1000 $ distSample (poisson lam) g + let (mean, variance) = meanVariance (map (fromIntegral . fromDiscrete) samples) + QM.assert $ (almostEqual tol mean (mu lam)) && + (almostEqual tol variance (var lam)) + where tol = 1e-1 + mu lam = lam + var lam = lam + +normalTest = do + mu <- QM.pick arbitrary + sd <- QM.pick $ choose (1, 10) + g <- QM.run $ MWC.create + let nsamples = floor (1000 * sd) -- larger standard deviations need more samples + -- to be shown as correct + samples <- QM.run $ replicateM nsamples $ distSample (normal mu sd) g + let (mean, variance) = meanVariance (map fromLebesgue samples) + QM.assert $ (almostEqual tol mean mu ) && + (almostEqual tol variance (var sd)) + where tol = 1e-1 + var sd = sq sd + +laplaceTest = do + mu <- QM.pick arbitrary + sd <- QM.pick $ choose (1, 10) + g <- QM.run $ MWC.create + let nsamples = floor (1000 * sd) -- larger standard deviations need more samples + -- to be shown as correct + samples <- QM.run $ replicateM nsamples $ distSample (laplace mu sd) g + let (mean, variance) = meanVariance (map fromLebesgue samples) + QM.assert $ (almostEqual tol mean mu ) && + (almostEqual tol variance (var sd)) + where tol = 1e-1 + var sd = 2*(sq sd) + +gammaTest = do + a <- QM.pick $ choose (1, 10) + b <- QM.pick $ choose (1, 10) + g <- QM.run $ MWC.create + samples <- QM.run $ replicateM 1000 $ distSample (gamma a b) g + let (mean, variance) = meanVariance (map fromLebesgue samples) + QM.assert $ (almostEqual tol mean (mu a b)) && + (almostEqual tol variance (var a b)) + where tol = 1e-1 + mu a b = a * b + var a b = a * (b * b) + diff --git a/Tests/ImportanceSampler.hs b/Tests/ImportanceSampler.hs new file mode 100644 --- /dev/null +++ b/Tests/ImportanceSampler.hs @@ -0,0 +1,118 @@ +{-# LANGUAGE BangPatterns #-} + +module Tests.ImportanceSampler where + +import Data.Dynamic +import Language.Hakaru.Types +import Language.Hakaru.Lambda +import Language.Hakaru.Distribution +import Language.Hakaru.ImportanceSampler + +-- import Test.QuickCheck.Monadic +import Tests.Models + +-- Some test programs in our language + +test_mixture :: IO () +test_mixture = sample prog_mixture conds >>= + print . take 10 >> + putChar '\n' >> + empiricalMeasure 1000 prog_mixture conds >>= + print + where conds = [Just (toDyn (Lebesgue 2 :: Density Double))] + +prog_dup :: Measure (Bool, Bool) +prog_dup = do + let c = unconditioned (bern 0.5) + x <- c + y <- c + return (x,y) + +prog_dbn :: Measure Bool +prog_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.90 else bern 0.10) + s2 <- unconditioned (if s1 then bern 0.75 else bern 0.25) + _ <- conditioned (if s2 then bern 0.90 else bern 0.10) + return s2 + +test_dbn :: IO () +test_dbn = sample prog_dbn conds >>= + print . take 10 >> + putChar '\n' >> + empiricalMeasure 1000 prog_dbn conds >>= + print + where conds = [Just (toDyn (Discrete True)), + Just (toDyn (Discrete True))] + +prog_hmm :: Integer -> Measure Bool +prog_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.90 else bern 0.10) + 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_hmm :: IO () +test_hmm = sample (prog_hmm 2) conds >>= + print . take 10 >> + putChar '\n' >> + empiricalMeasure 1000 (prog_hmm 2) conds >>= + print + where conds = [Just (toDyn (Discrete True)), + Just (toDyn (Discrete True))] + +prog_carRoadModel :: Measure (Double, Double) +prog_carRoadModel = do + speed <- unconditioned (uniform 5 15) + let z0 = lit 0 + _ <- conditioned (normal z0 1) + z1 <- unconditioned (normal (z0 + speed) 1) + _ <- conditioned (normal z1 1) + z2 <- unconditioned (normal (z1 + speed) 1) + _ <- conditioned (normal z2 1) + z3 <- unconditioned (normal (z2 + speed) 1) + _ <- conditioned (normal z3 1) + z4 <- unconditioned (normal (z3 + speed) 1) + return (z4, z3) + +test_carRoadModel :: IO () +test_carRoadModel = sample prog_carRoadModel conds >>= + print . take 10 >> + putChar '\n' >> + empiricalMeasure 1000 prog_carRoadModel conds >>= + print + where conds = [Just (toDyn (Lebesgue 0 :: Density Double)), + Just (toDyn (Lebesgue 11 :: Density Double)), + Just (toDyn (Lebesgue 19 :: Density Double)), + Just (toDyn (Lebesgue 33 :: Density Double))] + +prog_categorical :: Measure Bool +prog_categorical = do + rain <- unconditioned (categorical [(True, 0.2), (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.90 else bern 0.1)) + return rain + +test_categorical :: IO () +test_categorical = sample prog_categorical conds >>= + print . take 10 >> + putChar '\n' >> + empiricalMeasure 1000 prog_categorical conds >>= + print + where conds = [Just (toDyn (Discrete True))] + +prog_multiple_conditions :: Measure Double +prog_multiple_conditions = do + b <- unconditioned (beta 1 1) + _ <- conditioned (bern b) + _ <- conditioned (bern b) + return b diff --git a/Tests/Metropolis.hs b/Tests/Metropolis.hs new file mode 100644 --- /dev/null +++ b/Tests/Metropolis.hs @@ -0,0 +1,45 @@ +{-# LANGUAGE BangPatterns #-} + +module Tests.Metropolis where + +import Data.Dynamic + +import Language.Hakaru.Types +import Language.Hakaru.Lambda +import Language.Hakaru.Metropolis +import Language.Hakaru.Distribution (bern, normal, uniform, beta) + +-- import Test.QuickCheck.Monadic +-- import Distribution.TestSuite.QuickCheck +import Tests.Models + +test_mixture :: IO [Bool] +test_mixture = mcmc prog_mixture [Just (toDyn (Lebesgue (-2) :: Density Double))] + +test_multiple_conditions :: IO [Double] +test_multiple_conditions = do + mcmc prog_multiple_conditions [Just (toDyn (Discrete True)), + Just (toDyn (Discrete False))] + +prog_two_normals :: Measure Bool +prog_two_normals = unconditioned (bern 0.5) `bind` \coin -> + ifThenElse coin (conditioned (normal 0 1)) + (conditioned (normal 100 1)) `bind` \_ -> + return coin + +test_two_normals :: IO [Bool] +test_two_normals = mcmc prog_two_normals [Just (toDyn (Lebesgue 1 :: Density Double))] + +test_normal :: IO [Double] +test_normal = mcmc (unconditioned (normal 1 3)) [] + +prog_joint = do + bias <- unconditioned $ beta 1 1 + coin <- unconditioned $ bern bias + return (bias, coin) + +prog_condition = condition prog_joint True + +test_condition :: IO [Double] +test_condition = mcmc prog_condition [] + diff --git a/Tests/Models.hs b/Tests/Models.hs new file mode 100644 --- /dev/null +++ b/Tests/Models.hs @@ -0,0 +1,26 @@ +{-# LANGUAGE NoMonomorphismRestriction, GADTs #-} +module Tests.Models where + +import Prelude hiding (Real) + +import Language.Hakaru.Syntax2 +import Language.Hakaru.Lambda +import qualified Language.Hakaru.Types as T + +-- if we want to forgo the (D m) constraint, need to decorate the +-- program a little more. +prog_mixture :: (Meas m, D m ~ T.Dist ) => m Bool +prog_mixture = do + c <- unconditioned (bern 0.5) + _ <- conditioned (ifThenElse c (normal 1 1) + (uniform 0 3)) + return c + +prog_multiple_conditions :: (Meas m, D m ~ T.Dist) => m Double +prog_multiple_conditions = do + b <- unconditioned (beta 1 1) + _ <- conditioned (bern b) + _ <- conditioned (bern b) + return b + + diff --git a/Tests/Tests.hs b/Tests/Tests.hs new file mode 100644 --- /dev/null +++ b/Tests/Tests.hs @@ -0,0 +1,24 @@ +module Main where + +import Prelude hiding (Real) + +import qualified Tests.ImportanceSampler as IS +import qualified Tests.Metropolis as MH +import qualified Tests.Distribution as D + +import Language.Hakaru.Syntax2 +import Language.Hakaru.Lambda +import qualified Language.Hakaru.Types as T + +import Test.Framework (defaultMain, testGroup) +import Test.Framework.Providers.HUnit + +import Test.HUnit + +tests = [ + testGroup "Distribution checks" D.qtest, + testCase "alwaysPass" (1 @?= 1) + --testCase "alwaysFail" (error "Fail!") + ] + +main = defaultMain tests diff --git a/Types.hs b/Types.hs deleted file mode 100644 --- a/Types.hs +++ /dev/null @@ -1,13 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/Util/Coda.hs +++ /dev/null @@ -1,12 +0,0 @@ -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/Extras.hs b/Util/Extras.hs deleted file mode 100644 --- a/Util/Extras.hs +++ /dev/null @@ -1,84 +0,0 @@ -{-| - 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/Visual.hs b/Visual.hs deleted file mode 100644 --- a/Visual.hs +++ /dev/null @@ -1,43 +0,0 @@ -{-# 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 --- a/hakaru.cabal +++ b/hakaru.cabal @@ -2,10 +2,10 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: hakaru -version: 0.1.2 -synopsis: A probabilistic programming embedded DSL -description: Hakaru is an embedded DSL for performing probabilistic inference. It supports multiple inference backends. -homepage: http://www.indiana.edu/~ppaml +version: 0.1.3 +synopsis: A probabilistic programming embedded DSL +-- description: +homepage: http://indiana.edu/~ppaml/ license: BSD3 license-file: LICENSE author: The Hakaru Team @@ -17,9 +17,83 @@ cabal-version: >=1.10 library - exposed-modules: Sampler, Types, Visual, Language.Hakaru.Syntax, Mixture, Language.Hakaru.Symbolic, RandomChoice, Util.Extras, Util.Coda, Examples.Tests, 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 + exposed-modules: Language.Hakaru.Types, + Language.Hakaru.Symbolic, + Language.Hakaru.Arrow, + Language.Hakaru.Mixture, + Language.Hakaru.Sampler, + Language.Hakaru.ImportanceSampler, + Language.Hakaru.Metropolis, + Language.Hakaru.Lambda, + Language.Hakaru.Distribution, + Language.Hakaru.Util.Csv, + Language.Hakaru.Util.Extras, + Language.Hakaru.Util.Visual, + Language.Hakaru.Util.Coda + other-extensions: RankNTypes, BangPatterns, GADTs, TypeFamilies, TypeOperators, + ConstraintKinds, FlexibleContexts, NoMonomorphismRestriction, + DeriveDataTypeable, ScopedTypeVariables, ExistentialQuantification, + StandaloneDeriving, OverloadedStrings, + FlexibleInstances, RebindableSyntax + build-depends: base >=4.6 && <5.0, + random >=1.0 && <1.1, + transformers >=0.3 && <0.4, + containers >=0.5 && <0.6, + pretty >=1.1 && <1.2, + logfloat >=0.12 && <0.13, + hmatrix >=0.16 && <0.17, + math-functions >=0.1 && <0.2, + vector >=0.10 && <0.11, + cassava >=0.4 && <0.5, + zlib >=0.5 && <0.6, + bytestring >=0.10 && <0.11, + aeson >=0.7 && <0.8, + text >=1.1 && <1.2, + statistics >=0.11 && <0.14, + parsec >=3.1 && <3.2, + array >=0.4, + mwc-random >=0.13 && <0.14, + directory >=1.2 && <1.3, + integration >= 0.2.0 && < 0.3.0, + primitive >= 0.5 && < 0.6, + parallel >=3.2 && <3.3, + monad-loops >= 0.3.0.2 -- hs-source-dirs: default-language: Haskell2010 + ghc-options: -Wall + +test-suite hakaru-test + type: exitcode-stdio-1.0 + main-is: Tests/Tests.hs + other-modules: Tests.ImportanceSampler, + Tests.Models, + Tests.Metropolis, + Tests.Distribution + build-depends: base, Cabal >= 1.10, + QuickCheck, + HUnit, + test-framework, + test-framework-quickcheck2, + test-framework-hunit, + random >=1.0 && <1.1, + pretty >=1.1 && <1.2, + containers >=0.5 && <0.6, + logfloat >=0.12 && <0.13, + math-functions >=0.1 && <0.2, + statistics >=0.11 && <0.14, + hmatrix >=0.16 && <0.17, + vector >=0.10 && <0.11, + hakaru >= 0.1.3, + mwc-random >=0.13 && <0.14, + primitive >= 0.5 && < 0.6, + monad-loops >= 0.3.0.2 + default-language: Haskell2010 + +benchmark bench-all + type: exitcode-stdio-1.0 + hs-source-dirs: Bench + main-is: Bench.hs + build-depends: base, deepseq, ghc-prim, + criterion, hakaru >= 0.1.3 + ghc-options: -O2 + default-language: Haskell2010