diff --git a/COPYRIGHT b/COPYRIGHT new file mode 100644 --- /dev/null +++ b/COPYRIGHT @@ -0,0 +1,6 @@ +Copyright (c) 2011, Galois Inc. +Copyright (c) 2011, Battelle Memorial Institute +All rights reserved. + +The Passage library is distributed with the BSD3 license. See the LICENSE file +for details. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,43 @@ +Passage: A DSL for describing Bayesian Networks in Haskell + +Copyright (c) 2011, Galois, Inc. +Copyright (c) 2011, Battelle Memorial Institute +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 developer (Galois, Inc.) nor the + names of its 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 Galois Inc. 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. + +This material was prepared as an account of work sponsored by an agency of the +United States Government. Neither the United States Government nor the United +States Department of Energy, nor the Contractor, nor any or their employees, nor +any jurisdiction or organization that has cooperated in the development of these +materials, makes any warranty, express or implied, or assumes any legal +liability or responsibility for the accuracy, completeness, or usefulness or any +information, apparatus, product, software, or process disclosed, or represents +that its use would not infringe privately owned rights. +PACIFIC NORTHWEST NATIONAL LABORATORY +operated by +BATTELLE +for the +UNITED STATES DEPARTMENT OF ENERGY +under Contract DE-AC05-76RL01830 + diff --git a/Language/Passage.hs b/Language/Passage.hs new file mode 100644 --- /dev/null +++ b/Language/Passage.hs @@ -0,0 +1,192 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +module Language.Passage ( + + -- * Types + BayesianNetwork, BayesianGraph(..), StoVar(..), Node, BayesianSimulator + + -- * Constructing models + , logGamma, using, tconst + , vector, matrix, nodeArray, Vector, Matrix, NodeArray + , tcase + , (//) + + -- * Distributions + , module UI + + -- * Extracting graphs + , buildBayesianGraph + + -- * Displayihng graphs + , PP(..), LaTeX(..) + + -- * Simulation + , simulate, genSimulator + , setSampleCount + , setIterationsPerSample + , setWarmupCount + , setThreadNum + , useMersenneTwister + , enableProfiling + , setRandomSeed + , useSpecialSlicers + , splitFiles + , model, observe, monitor, monitorVec, monitorVecs + + -- * LaTeX + , runLatex + ) where + +import Language.Passage.AST +import Language.Passage.UI as UI +import Language.Passage.Graph +import Language.Passage.Lang.LaTeX(LaTeX(..)) +import qualified Language.Passage.Lang.LaTeX as LaTeX +import Language.Passage.Term +import Language.Passage.Utils +import Language.Passage.SimulatorConf + +import qualified Language.Passage.Graph2C as C + +import Control.Exception (finally) +import Control.Monad(when) + +import System.Process(rawSystem, readProcess) +import System.Exit(ExitCode(..)) +import System.Info(os) +import System.FilePath +import System.IO(openFile,hPutStrLn,hClose,IOMode(..)) +import System.Directory(removeDirectoryRecursive, doesDirectoryExist) +import Paths_passage (getDataDir) + +-- | Like monitor, but adds the indexes in the label of the variable. +monitorVec :: String -> Matrix -> [Int] -> BayesianSimulator () +monitorVec name m xs = monitor lab (m (map fromIntegral xs)) + where + lab = name ++ concatMap ix xs + ix x = "[" ++ show x ++ "]" + +monitorVecs :: String -> NodeArray -> [[Int]] -> BayesianSimulator () +monitorVecs name m = mapM_ (monitorVec name m) + +type Node = Expr + +withTempDir :: (FilePath -> IO a) -> IO a +withTempDir f = + do dir <- init `fmap` readProcess "mktemp" ["-d","-t","bayesiandsl.XXXXXX"] "" + -- init drops \n + f dir `finally` removeDirectoryRecursive dir + +runLatex :: BayesianNetwork a -> IO () +runLatex t = withTempDir $ \dir -> + do let file = dir "out" + tex_file = file <.> ".tex" + pdf_file = file <.> ".pdf" + writeFile tex_file (show doc) + runCmd make_pdf ["-output-directory", dir, tex_file] + runCmd show_pdf (pdf_args ++ [pdf_file]) + + where + doc = vcat [ LaTeX.cmd "documentclass" [ text "article" ] + , LaTeX.env "document" [] (latex (snd (buildBayesianGraph t))) + ] + + (show_pdf, pdf_args) + | os == "linux" = ("evince",[]) + | otherwise = ("open",["-W"]) + + make_pdf = "pdflatex" + + + +runCmd :: String -> [String] -> IO () +runCmd f as = + do res <- rawSystem f as + case res of + ExitSuccess -> return () + ExitFailure n -> + fail $ "(error " ++ show n ++ ") Failed to execute " ++ show f + ++ " with arguments " ++ show as + + + +createSimProject :: FilePath -> SimState -> C.SamplerConf -> IO () +createSimProject dir st conf = + do yes <- doesDirectoryExist dir + when yes $ error $ "Directory: " ++ show dir ++ " already exists." + + -- Copy templates + putStrLn $ "Creating directory " ++ show dir + dataDir <- getDataDir + let rt = dataDir "cbits" "runtime" + runCmd "cp" [ "-r", rt, dir ] + + -- Create additional settings + let src_dir = dir "src" + extra_settings = src_dir "extra_settings" + hExtra <- openFile extra_settings WriteMode + hPutStrLn hExtra "# Here one can put additional settings for the build" + when (cfgMersenne st) $ hPutStrLn hExtra "CPPFLAGS+=D__USE_MERSENNE" + when (cfgProfile st) $ hPutStrLn hExtra "CFLAGS+=-pg -g" + + hClose hExtra + + -- Generate simulator + putStrLn "Generating sampler." + -- let c_file = src_dir "sampler" <.> ".c" + mapM_ (\(f,d) -> writeFile (src_dir f) (show d)) $ C.gen_c conf + putStrLn $ "Generated C project: " ++ show src_dir + + -- generate R driver + let rDriver = dir "histogram.R" + writeFile rDriver (genR dir (zip [1..] (map fst (C.monitor conf)))) + putStrLn $ "Generated sample R cmds: " ++ show rDriver + +-- TODO: Generate an R driver that knows what's being observed +genR :: String -> [(Int, String)] -> String +genR name labs = unlines $ [ "library(MASS)" + , "pdf(file='sample.pdf')" + , "table <- read.table('datafile')" + ] ++ map genHist labs + where + genHist (i, s) = + "truehist(table[," ++ show i ++ "], xlab='" ++ s ++ + "', main='" ++ name ++ "')" + + + + + +createSimulator :: FilePath -> SimState -> IO () +createSimulator path st = + case cfgNetwork st of + Nothing -> error $ "No bayesian-network specified; please use \"bayesianNetwork\" to specify one." + Just t -> + case cfgMonitor st of + [] -> error $ "No montitors added; please use \"monitor\" to specify some." + ms -> + let conf = C.SamplerConf { C.graph = t + , C.sampleNum = cfgSampleNum st + , C.itsPerSample = cfgItsPerSample st + , C.warmup = cfgWarmup st + , C.seed = cfgRandomSeed st + , C.observe = cfgObserve st + , C.initialize = cfgInitialize st + , C.monitor = reverse ms + , C.thread_num = cfgThreadNum st + , C.special_slicers = cfgSpecialSlicers st + , C.split_files = cfgSplitFiles st + } + in createSimProject path st conf + +genSimulator :: FilePath -> BayesianSimulator () -> IO () +genSimulator f b = createSimulator f (runSim b) >> return () + +simulate :: FilePath -> BayesianSimulator () -> IO () +simulate f b = do createSimulator f (runSim b) + putStrLn "Running the simulation.." + runCmd "make" [ "--quiet", "-C", f ] + putStrLn "Done." + + + diff --git a/Language/Passage/AST.hs b/Language/Passage/AST.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/AST.hs @@ -0,0 +1,345 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE PatternGuards #-} + +module Language.Passage.AST where + +import qualified Data.Map as M +import qualified Data.IntMap as IM +import qualified Data.IntSet as IS +import qualified Data.Array as A +import MonadLib +import Data.Char(isUpper) + + +import Language.Passage.Utils +import Language.Passage.Lang.LaTeX(LaTeX(..)) +import qualified Language.Passage.Lang.LaTeX as LaTeX +import Language.Passage.Term + +------------------------------------------------------------------------ +-- * AST Nodes in a Bayesian model description +------------------------------------------------------------------------ + +type Expr = Term NodeIdx + +-- | A Bayesian Network is a collection of stochastic nodes. +-- Stochastic nodes may be (optionally) grouped into arrays. +data BayesianGraph = BayesianGraph + { stoNodes :: !(IM.IntMap StoVar) + , stoArryas :: !(IM.IntMap ArrayInfo) + } deriving Show + +-- | A Stochastic variable. +data StoVar = StoVar + { stoVarName :: StoVarName + , stoVarPrior :: PriorInfo + , stoPostDistLL :: !(M.Map Expr Expr) + -- ^ Maps terms that mention the variable to their coefficients, + -- which do not depend on the variable. The term for the + -- distribution is the sum of the products of the map elements + -- (see 'stoPostLL'). + } deriving Show + +-- | The name of a stochastic variable. +data StoVarName + = Unnamed !NodeIdx -- ^ Anonymous stand-alone variable + | InArray !NodeIdx ![Int] -- ^ This sto var belongs to an array + | Named !String -- ^ Standalone variable, with a user-name + deriving Show + +-- | Information about the prior distribution of a stochastic variable. +data PriorInfo = PriorInfo + { priName :: String + , priParams :: [Term NodeIdx] + , priSupport :: DistSupport NodeIdx + , priLL :: Term NodeIdx + } deriving Show + + +-- | The description of an atomic distribution. +data Distribution = Distribution + { distName :: String + , distParams :: [Expr] + , distSupport :: DistSupport NodeIdx + , distLL :: Expr -> Expr + } + + +-- | Support of a distribution +data DistSupport a + = Real + | Discrete (Maybe (Term a)) -- upper bound, from 0 to this number + | PosReal + | Interval (Term a) (Term a) -- lower/upper + deriving (Show,Functor) + + +-- | Information about an array (vecotr/matrix) of stochastic variables. +data ArrayInfo = ArrayInfo + { arrayName :: String + , arrayDimensions :: [(Int,Int)] + , arrayVars :: IS.IntSet + } deriving Show + + +exprToVar :: Expr -> NodeIdx +exprToVar (TVar x) = x +exprToVar e = error $ "Expected a variable expression, got: " ++ + (show (pp e)) + +fvsSupport :: ArrVars -> DistSupport NodeIdx -> IS.IntSet +fvsSupport arr sup = + case sup of + Real -> IS.empty + Discrete t -> maybe IS.empty (leavesOfTerm arr) t + PosReal -> IS.empty + Interval a b -> IS.union (leavesOfTerm arr a) (leavesOfTerm arr b) + + +fvsArray :: BayesianGraph -> NodeIdx -> IS.IntSet +fvsArray bg ix = case IM.lookup ix (stoArryas bg) of + Just ai -> arrayVars ai + Nothing -> IS.empty -- XXX: report error? + + + +latexDist :: LaTeX a => String -> [Term a] -> Doc +latexDist name params = fun <+> commaSep (map latex params) + where fun = case name of + [n] | isUpper n -> LaTeX.mathcal (char n) + _ -> LaTeX.mathrm (text name) + + +------------------------------------------------------------------------ +-- * Constructing the program +------------------------------------------------------------------------ +data ASTState = ASTState + { curIdx :: !Int + , declaredArrays :: !(IM.IntMap ArrayInfo) + , generatedNodes :: !(IM.IntMap StoVar) + } + +newtype BayesianNetwork a = BayesianNetwork (StateT ASTState Id a) + deriving (Functor,Monad) + +updateState :: (ASTState -> (a,ASTState)) -> BayesianNetwork a +updateState f = BayesianNetwork (sets f) + +updateState_ :: (ASTState -> ASTState) -> BayesianNetwork () +updateState_ f = updateState (\s -> let s1 = f s in seq s1 ((), s1)) + +getState :: BayesianNetwork ASTState +getState = BayesianNetwork get + + +using :: Distribution -> BayesianNetwork Expr +using d = updateState $ \st -> + let i = curIdx st + in ( tvar i + , st { curIdx = i + 1 + , generatedNodes = IM.insert i (toStoVar i d) (generatedNodes st) + } + ) + +toStoVar :: NodeIdx -> Distribution -> StoVar +toStoVar i d = StoVar + { stoVarName = Unnamed i + , stoVarPrior = PriorInfo + { priName = distName d + , priParams = distParams d + , priSupport = distSupport d + , priLL = distLL d (tvar i) + } + , stoPostDistLL = M.empty + } + +-- NOTE: The posterior LLs are not yet computed in the graph that's returned. +extractNetwork :: BayesianNetwork a -> (a, BayesianGraph) +extractNetwork (BayesianNetwork m) = + ( a + , BayesianGraph { stoNodes = generatedNodes s + , stoArryas = declaredArrays s + } + ) + where (a, s) = runId (runStateT start m) + start = ASTState { curIdx = 0 + , declaredArrays = IM.empty + , generatedNodes = IM.empty + } + +type Vector = [Expr] -> Expr +type Matrix = [Expr] -> Expr +type NodeArray = [Expr] -> Expr + +-- | Create a 1D vector +vector :: (Int,Int) -- ^ Bounds for the vector indexes + -> (Int -> BayesianNetwork Expr) -- ^ Initializer (should return a node) + -> BayesianNetwork ([Expr] -> Expr) +vector b i = nodeArray [b] (i . head) + +-- | Create a 2D matrix +matrix :: (Int,Int) -> (Int,Int) -- ^ Bounds for 1st and 2nd dimensions. + -> ([Int] -> BayesianNetwork Expr) -- ^ Initializer + -> BayesianNetwork ([Expr] -> Expr) +matrix b1 b2 = nodeArray [b1, b2] + +-- | Create an >= 3D array +nodeArray :: [(Int,Int)] -- ^ Bounds for each dimension. + -> ([Int] -> BayesianNetwork Expr) -- ^ Initializer + -> BayesianNetwork ([Expr] -> Expr) +nodeArray bds initializer = + do (ix,ai,mp) <- newArray bds initializer + return (lkpArrayMap ix ai mp) + + + +data ArrayMap = A !(A.Array Int ArrayMap) + | V !Expr + +lkpArrayMap :: NodeIdx -> ArrayInfo -> ArrayMap -> [Expr] -> Expr +lkpArrayMap x0 _ am = loop (tarr x0) am + where + loop _ (V x) [] = x + loop e (A _) [] = e + loop e (A a) (i : is) + | Just j <- toIx i = loop (tIx e i) (a A.! j) is + loop e _ is = foldl tIx e is -- XXX: could do some more checking! + + toIx (TConst d) = Just (floor d) -- XXX: it'd be better to use proper types... + toIx _ = Nothing + + +newArray :: [(Int,Int)] + -> ([Int] -> BayesianNetwork Expr) + -> BayesianNetwork (NodeIdx, ArrayInfo, ArrayMap) +newArray bds0 initializer = + do aix <- updateState $ \st -> let i = curIdx st in (i, st { curIdx = 1 + i }) + let bds = map dimOK bds0 + (vars,m) <- loop aix IS.empty bds [] + let ai = ArrayInfo { arrayName = "a" ++ show aix + , arrayDimensions = bds + , arrayVars = vars + } + + updateState $ \s -> + ( (aix, ai, m) + , s { declaredArrays = IM.insert aix ai (declaredArrays s) } + ) + + where + dimOK d@(x,y) | x <= y = d + dimOK d = longError [ "Invalid array bounds:" + , " *** Bounds: " ++ show d + ] + + -- The state of "loop" + -- aix: RO, index of the array that we are initializing + -- vars: RW, a set of all the variables in the array (to be stored for later) + -- bds: RW, remaining array dimnesions to process + -- ixes: RW, (reversed) path of indexes to current elent which to initialize + + -- We pass vars explicitly, rather then putting it in the state to + -- avoid constant updates to the array map. + loop aix vars [] ixes0 = + do let ixes = reverse ixes0 + e <- initializer ixes + let v = exprToVar e + updateState $ \st -> + let vars1 = IS.insert v vars + in vars1 `seq` + ( (vars1, V e) + , let upd sv = + case stoVarName sv of + Unnamed _ -> sv { stoVarName = InArray aix ixes } + InArray a' ixes' -> longError + [ "Variable already belongs to an array:" + , " *** Array: " ++ lkpArrayName st a' + , " *** Location: " ++ show ixes' + ] + Named s -> longError + [ "Cannot add explicitly named variables to an array:" + , " *** Variable: " ++ s + ] + in st { generatedNodes = IM.adjust upd v (generatedNodes st) } + ) + + loop aix vars0 (bds@(from,to) : bdss) ixes = + + let loop1 vars as i | i <= to = + do (vars1,a) <- loop aix vars bdss (i:ixes) + loop1 vars1 (a:as) (i+1) + loop1 vars as _ = return (vars, A $ A.array bds + $ zip [ from .. to ] (reverse as)) + in loop1 vars0 [] from + + + + + + + +longError :: [String] -> a +longError = error . unlines + + +infixl 1 // + +-- This operator is used to provide a custom name for a given variable. +(//) :: BayesianNetwork Expr -> String -> BayesianNetwork Expr +m // x = + do e <- m + let v = exprToVar e + updateState $ \s -> (e, newState s v) + + where + newState s v = + case IM.lookup v (generatedNodes s) of + -- XXX: This looks up the name twice. + Just sv -> s { generatedNodes = IM.insert v (setName s sv) + (generatedNodes s) } + Nothing -> + case IM.lookup v (declaredArrays s) of + Just ai -> s { declaredArrays = + IM.insert v ai { arrayName = x } + (declaredArrays s) } + Nothing -> longError + [ "Attempt to rename an unknown node:" + , " *** Node: " ++ show v + ] + + setName s sv = case stoVarName sv of + Unnamed _ -> sv { stoVarName = Named x } + Named n -> longError + [ "Cannot rename a variable multiple times: " + , " *** old name: " ++ n + , " *** new name: " ++ x + ] + InArray a is -> longError + [ "Cannot rename array vairable: " + , " *** array: " ++ lkpArrayName s a ++ show is + , " *** new name" ++ x + ] + + +lkpArrayName :: ASTState -> NodeIdx -> String +lkpArrayName st a = case IM.lookup a (declaredArrays st) of + Nothing -> "(unknown?)" + Just ai -> arrayName ai + + +{- + +Do we need to support arrays of deterministic variables? +Example: + +do x <- bernoulli 0.5 + a <- detVector (1,100) $ \i -> 2 * i + return (a ! x) + +This requires prpoer support for determinsitc nodes, +which is broken at present. + +-} + + diff --git a/Language/Passage/Distribution.hs b/Language/Passage/Distribution.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/Distribution.hs @@ -0,0 +1,158 @@ +module Language.Passage.Distribution where + +import Language.Passage.AST +import Language.Passage.Term(logGamma, tcase) + +logit :: Floating a => a -> a +logit p = log(p/(1-p)) + +logBeta :: Expr -> Expr -> Expr +logBeta x y = logGamma x + logGamma y - logGamma (x + y) + +logFact :: Expr -> Expr +logFact n = logGamma (n + 1) + +logComb :: Expr -> Expr -> Expr +logComb n k = logFact n - logFact k - logFact (n - k) + +-- | A normal distribution with mean 0 and precision 1 +stdNormal :: Distribution +stdNormal = Distribution + { distName = "N(0,1)" + , distParams = [] + , distSupport = Real + , distLL = \x -> -0.5 * x**2 + } + +-- | A normal distribution, with a mean and precision +normal :: Expr -> Expr -> Distribution +normal m t = Distribution + { distName = "N" + , distParams = [m, t] + , distSupport = Real + , distLL = \x -> log t / 2 + - t * (x ** 2) / 2 + + t * x * m + - t * (m ** 2) / 2 + } + +--- | A standard uniform distribution with parameters 0 and 1 +standardUniform :: Distribution +standardUniform = Distribution + { distName = "SU" + , distParams = [0, 1] + , distSupport = Interval 0 1 + , distLL = \_ -> 0 + } + +--- | A uniform distribution with lower and upper bounds +uniform :: Expr -> Expr -> Distribution +uniform lo hi = Distribution + { distName = "U" + , distParams = [lo, hi] + , distSupport = Interval lo hi + -- NB: Uniform distribution, independent of the variable (hence constant function) + , distLL = \_ -> - (log (hi - lo)) + } + +discreteUniform :: Expr -> Distribution +discreteUniform n = Distribution + { distName = "DisreteUniform" + , distParams = [0, n] + , distSupport = Discrete (Just n) + -- NB: Uniform distribution, independent of the variable (hence constant function) + , distLL = \_ -> - (log (n + 1)) + } + +geometric :: Expr -> Distribution +geometric p = Distribution + { distName = "Geometric" + , distParams = [p] + , distSupport = Discrete Nothing + , distLL = \x -> x * log (1 - p) + log p + } + +-- | A categorical distribution with given support size and probabilities +-- | Probabilities are assumed to add to one (not checked here) +categorical :: Expr -> [Expr] -> Distribution +categorical n ps = Distribution + { distName = "Categorical" + , distParams = n:ps + , distSupport = Discrete (Just (n - 1)) + , distLL = \x -> log (tcase x ps) + } + +-- | A Bernoulli distribution with a mean +bernoulli :: Expr -> Distribution +bernoulli p = Distribution + { distName = "B" + , distParams = [p] + , distSupport = Discrete (Just 1) + , distLL = \x -> log (1 - p) + logit p * x + } + +-- | A binomial distribution with given number of samples and probability of success +-- | Number of samples is assumed to be fixed +binomial :: Expr -> Expr -> Distribution +binomial n p = Distribution + { distName = "Binomial" + , distParams = [n, p] + , distSupport = Discrete (Just n) + , distLL = \x -> logComb n x + x * logit p + n * log (1 - p) + } + +negBinomial :: Expr -> Expr -> Distribution +negBinomial r p = Distribution + { distName = "NegativeBinomial" + , distParams = [r, p] + , distSupport = PosReal + , distLL = \x -> logComb (x+r-1) x + r * log (1 - p) + x * log p + } + +poisson :: Expr -> Distribution +poisson lambda = Distribution + { distName = "Poisson" + , distParams = [lambda] + , distSupport = Discrete Nothing + , distLL = \x -> x * log lambda - logFact x - lambda + } + +-- | A beta distribution with the given prior sample sizes. +beta :: Expr -> Expr -> Distribution +beta a b = + Distribution + { distName = "Beta" + , distParams = [a, b] + , distSupport = Interval 0 1 + , distLL = \x -> (a - 1) * log x + (b - 1) * log (1 - x) - logBeta a b + } + +-- | A gamma distribution with the given prior sample sizes. +dgamma :: Expr -> Expr -> Distribution +dgamma a b = + Distribution + { distName = "Gamma" + , distParams = [a, b] + , distSupport = PosReal + , distLL = \x -> a * log b - logGamma a + (a - 1) * log x - b * x + } + +-- | An improper uniform distribution; has no impact on likelihood +improperUniform :: Distribution +improperUniform = + Distribution + { distName = "ImproperUniform" + , distParams = [] + , distSupport = Real + , distLL = const 0 + } + +-- | An improper scale +improperScale :: Distribution +improperScale = + Distribution + { distName = "ImproperScale" + , distParams = [] + , distSupport = PosReal + , distLL = \x -> -log x + } diff --git a/Language/Passage/Graph.hs b/Language/Passage/Graph.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/Graph.hs @@ -0,0 +1,128 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Language.Passage.Graph where + +import qualified Data.IntMap as IM +import qualified Data.Map as M +import qualified Data.IntSet as IS +import Data.List(foldl') + +-- import Debug.Trace + +import Language.Passage.AST +import Language.Passage.Term +import Language.Passage.Utils +import Language.Passage.Lang.LaTeX(LaTeX(..)) +import qualified Language.Passage.Lang.LaTeX as LaTeX + + +stoPostLL :: StoVar -> Term NodeIdx +stoPostLL sv = sum [ b * a | (a,b) <- M.toList (stoPostDistLL sv) ] + +emptyBayesianGraph :: BayesianGraph +emptyBayesianGraph = BayesianGraph { stoNodes = IM.empty + , stoArryas = IM.empty + } + +addToStoLL :: NodeIdx -> Term NodeIdx -> BayesianGraph -> BayesianGraph +addToStoLL ix t bg = bg { stoNodes = IM.alter addLL ix (stoNodes bg) } + where + (x,c) = factorVar (fvsArray bg) ix t + -- addLL sv = sv { stoPostDistLL = M.insertWith plus x c (stoPostDistLL sv) } + addLL (Just sv) = Just $! sv { stoPostDistLL = M.insertWith' plus x c (stoPostDistLL sv) } + + addLL Nothing = Nothing + + plus :: (PP a, Eq a, Show a) => Term a -> Term a -> Term a + plus a b = {- trace ("plus: " ++ "\n a: " ++ show (pp a) + ++ "\n a: " ++ show a + ++ "\n b: " ++ show (pp b) + ++ "\n b: " ++ show b + ++ "\n a+b: " ++ show (pp result)) -} + result + where result = maybe (a+b) id (sAdd a b) +-------------------------------------------------------------------------------- + +buildBayesianGraph :: BayesianNetwork a -> (a, BayesianGraph) +buildBayesianGraph nw = (a, computeLL g) + where (a, g) = extractNetwork nw + +-- | Compute the log-likelihood for a stochastic variable. +computeLL :: BayesianGraph -> BayesianGraph +computeLL bg = foldl' addDef bg (IM.elems (stoNodes bg)) + where addDef m sv = foldl' addSum m (summands (priLL (stoVarPrior sv))) + addSum m t = IS.fold (\i m1 -> addToStoLL i t m1) m + (leavesOfTerm (fvsArray bg) t) + + + + +-------------------------------------------------------------------------------- +-- Pretty printing +-------------------------------------------------------------------------------- + +data PPVar = PPName String + | PPArr String [Int] + deriving Show + +nameToPPName :: BayesianGraph -> StoVar -> PPVar +nameToPPName bg sv = + case stoVarName sv of + Unnamed y -> PPName ("v" ++ show y) + Named y -> PPName y + InArray a b -> + case IM.lookup a (stoArryas bg) of + Just ai -> PPArr (arrayName ai) b + Nothing -> PPArr ("bug_unknown_array_" ++ show a) b + +varName :: BayesianGraph -> NodeIdx -> PPVar +varName bg x = case IM.lookup x (stoNodes bg) of + Just sv -> nameToPPName bg sv + Nothing -> + case IM.lookup x (stoArryas bg) of + Just ai -> PPName (arrayName ai) + Nothing -> PPName ("bug_unknown_variable_" ++ show x) + +namedTerm :: BayesianGraph -> Term NodeIdx -> Term PPVar +namedTerm bg = fmap (varName bg) + +instance PP PPVar where + pp (PPName x) = text x + pp (PPArr x ys) = text x <> hcat (map (brackets . int) ys) + +instance LaTeX PPVar where + latex (PPName x) = LaTeX.var x + latex (PPArr x ys) = LaTeX.var x <> char '_' <> + braces (hcat (punctuate comma (map int ys))) + + +instance PP BayesianGraph where + pp bg = vcat (map ppSto (IM.elems (stoNodes bg))) + + where + ppT t = pp (namedTerm bg t) + ppSto sv = pp (nameToPPName bg sv) <+> text "~~" <+> + ppPri (stoVarPrior sv) <+> + text ":" <+> ppT (stoPostLL sv) + + ppPri i = text (priName i) <+> + commaSep (map (pp . namedTerm bg) (priParams i)) + +instance LaTeX BayesianGraph where + latex bg = + LaTeX.env "tabular" [text "l"] $ vcat $ map (\x -> LaTeX.row [x]) + [ LaTeX.env "tabular" [text "l l"] + (LaTeX.row [ text "Prior distribution" + , text "Posterior log-likelihood" ] $$ + vcat (map ppSto (IM.elems (stoNodes bg)))) + ] + + where ppT t = latex (namedTerm bg t) + row x y z = LaTeX.row (map LaTeX.math [ x <+> LaTeX.sim <+> y, z]) + ppSto sv = + row (latex (nameToPPName bg sv)) + (ppPri (stoVarPrior sv)) + (ppT (stoPostLL sv)) + + ppPri i = latexDist (priName i) (map (namedTerm bg) (priParams i)) + + diff --git a/Language/Passage/Graph2C.hs b/Language/Passage/Graph2C.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/Graph2C.hs @@ -0,0 +1,757 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +module Language.Passage.Graph2C where + +import Language.Passage.Utils hiding (double,int) +import Language.Passage.Term hiding (bin) +import Language.Passage.AST +import Language.Passage.Lang.C +import Language.Passage.GraphColor(groupByColor) + + +import qualified Data.Map as M +import qualified Data.IntMap as IM +import qualified Data.IntSet as IS +import Data.Maybe(maybeToList,fromJust) +import Data.List(sortBy, transpose) +import Data.Function(on) +import MonadLib (ReaderT, StateT, Id + , runId, runStateT, runReaderT + , get, set, asks, mapReader + , forM + , zipWithM + ) + +import Data.Graph(SCC(..)) +import Data.Graph.SCC + + +-------------------------------------------------------------------------------- +-- Compilation of expressions +-------------------------------------------------------------------------------- + +cnameVar :: (NodeIdx, StoVar) -> M CExpr +cnameVar (ix,sv) = + case stoVarName sv of + InArray x is -> + do ai <- lookupArray x + let fixIndex i (from,_) = int_lit (i - from) + name = arrName (x,ai) + return $ foldl arr_ix (var name) + $ zipWith fixIndex is $ arrayDimensions ai + _ -> return $ var $ simpleCName ix + + +simpleCName :: NodeIdx -> CIdent +simpleCName x = ident ("v_" ++ show x) + +variable :: NodeIdx -> M CExpr +variable x = + + -- Is this an observed variable? + do isObs <- isObserved x + case isObs of + Just v -> return (double_lit v) -- Yes, just use the known value. + Nothing -> + + {- Are we compiling within the LL_FUN for this variable? + When we generate the LL_FUN for a stoachastic variable, + we always use it's simple name, even if the variable is stored + in an array in the long run. The reason for this is that in + the LL_FUN, the variable is passed as an argument. -} + do samp <- isSampled x + if samp + then return $ var $ simpleCName x + else + + -- OK, perhaps we have an ordinary stochastic variable? + do mbsv <- lookupVarMb x + case mbsv of + Just sv -> cnameVar (x,sv) + + -- Hmm, we don't know about this variable. + -- The variable must refer to a deterministic node + -- generated to factor repeated compution out + -- of an LL_FUN. + Nothing -> return $ var $ simpleCName x + +term :: Term NodeIdx -> M CExpr +term t = + case t of + TVar x -> variable x + TArr x -> do ai <- lookupArray x + return $ var $ arrName (x,ai) + TConst x -> return (double_lit x) + TApp op ts -> + do ds@(a : bs) <- mapM term ts + let b : _ = bs + bin x = parens a <+> text x <+> parens b + case op of + TExp -> return $ call (ident "exp") ds + TLog -> return $ call (ident "log") ds + TNeg -> return $ char '-' <> parens (head ds) + TAdd -> return $ bin "+" + TMul -> return $ bin "*" + TSub -> return $ bin "-" + TDiv -> return $ bin "/" + TPow -> case ts of + [_ , TConst 2.0] -> (return $ call (ident "square") [a]) + _ -> return $ call (ident "pow") ds + TLogGamma -> return $ call (ident "lgamma") ds + + TCase -> + do i <- newDetVar -- Just used as a new name + let name = ident ("case_fun_" ++ show i) + + -- if we are in the LL function for some variable, + -- we have to pass the sampled variable to the "case" function. + args <- (map simpleCName . maybeToList) `fmap` isSampling + + newLocalFunDecl + (fun_decl double name [ (double,x) | x <- args ]) [ + switch (cast int a) + (zip [ 0 .. ] (map (return . creturn) bs)) + (callS (var (ident "crash_out_of_bounds")) + [ text "__LINE__" ]) + ] + return $ call (var name) (map var args) + + TIx -> + case ts of + [ arr, ix ] -> + do dims <- getArrDimensions arr + case dims of + (from,_) : _ -> + do expr <- term (ix - fromIntegral from) + return (arr_ix a (cast int expr)) + _ -> error $ "Type error: attempt ro index non an array." + _ -> error $ "TIx: Unexpected args: " ++ show ts + +getArrDimensions :: Term NodeIdx -> M [(Int,Int)] +getArrDimensions t = + case t of + TArr x -> arrayDimensions `fmap` lookupArray x + TApp TIx [ a, _ ] -> + do ds <- getArrDimensions a + case ds of + _ : ds1 -> return ds1 + [] -> error $ "Type error: attempt to index a non-array." + _ -> error $ "Type error: not an array" + +-------------------------------------------------------------------------------- + +newtype M a = M (ReaderT R (StateT S Id) a) deriving (Functor, Monad) + +data R = R { config :: SamplerConf + , sampling :: Maybe NodeIdx + } + +data CModule = + CModule + { cpp_stuff :: Doc -- ^ Includes, #define, etc. + , var_decls :: [Doc] -- ^ Variable declarations + , cpp_funs :: Doc -- ^ #included templtes + , fun_decls :: [(Doc,Doc)] -- ^ Function declarations: decl, body + } + +blankMod :: CModule +blankMod = CModule { cpp_stuff = empty + , var_decls = [] + , cpp_funs = empty + , fun_decls = [] + } + +-- XXX: Watch out with the ++ing here... +mergeCModules :: CModule -> CModule -> CModule +mergeCModules m1 m2 = + CModule { cpp_stuff = cpp_stuff m1 $$ cpp_stuff m2 + , var_decls = var_decls m1 ++ var_decls m2 + , cpp_funs = cpp_funs m1 $$ cpp_funs m2 + , fun_decls = fun_decls m1 ++ fun_decls m2 + } + + + +data S = S { main_mod :: CModule + , cur_mod :: Maybe CModule + , helper_mods :: [(NodeIdx, CModule)] + , cnames :: !Int -- ^ Name supply + } + +noHelpers :: S -> S +noHelpers s = s { main_mod = foldr mergeCModules (main_mod s) + $ map snd (helper_mods s) + , helper_mods = [] + } + +getGraph :: M BayesianGraph +getGraph = M (asks (graph . config)) + +lookupArray :: NodeIdx -> M ArrayInfo +lookupArray x = + do mb <- lookupArrayMb x + case mb of + Just a -> return a + Nothing -> error ("Unknown array variable: " ++ show x) + + +lookupArrayMb :: NodeIdx -> M (Maybe ArrayInfo) +lookupArrayMb x = + do g <- getGraph + return (IM.lookup x (stoArryas g)) + +lookupVarMb :: NodeIdx -> M (Maybe StoVar) +lookupVarMb x = + do g <- getGraph + return (IM.lookup x (stoNodes g)) + +lookupVar :: NodeIdx -> M StoVar +lookupVar x = + do mb <- lookupVarMb x + case mb of + Just sv -> return sv + Nothing -> error ("Unknown stochastic variable: " ++ show x) + +nowSampling :: NodeIdx -> M a -> M a +nowSampling x (M a) = M (mapReader (\i -> i { sampling = Just x }) a) + +isSampling :: M (Maybe NodeIdx) +isSampling = M (asks sampling) + +isSampled :: NodeIdx -> M Bool +isSampled x = (Just x ==) `fmap` isSampling + +isObserved :: NodeIdx -> M (Maybe Double) +isObserved ix = IM.lookup ix `fmap` M (asks (observe . config)) + +isInitialized :: NodeIdx -> M (Maybe Double) +isInitialized ix = IM.lookup ix `fmap` M (asks (initialize . config)) + +newDetVar :: M NodeIdx +newDetVar = M $ + do s <- get + let i = cnames s + set s { cnames = i + 1 } + return i + +newHelper :: NodeIdx -> M a -> M a +newHelper i (M m) = M $ + do s <- get + set s { cur_mod = Just blankMod } + a <- m + s1 <- get + set s1 { cur_mod = Nothing + , helper_mods = (i, fromJust (cur_mod s1)) : helper_mods s1 + } + return a + +updHelper :: (CModule -> CModule) -> M () +updHelper f = M $ + do s <- get + case cur_mod s of + Nothing -> error "BUG: updHelper called without a module" + Just m -> + set s { cur_mod = Just (f m) } + +updMain :: (CModule -> CModule) -> M () +updMain f = M $ + do s <- get + set s { main_mod = f (main_mod s) } + + +-- add a new function to the main module. +newFunDecl :: CFunDecl -> [CStmt] -> M () +newFunDecl d body = + updMain $ \m -> m { fun_decls = (d, block body) : fun_decls m } + +-- add a new declaration to the main module. +newDecl :: CDecl -> M () +newDecl d = updMain $ \m -> m { var_decls = d : var_decls m } + + +-- Add "cpp" includes to the main module +cpp :: String -> M () +cpp t = updMain $ \m -> m { cpp_stuff = cpp_stuff m $$ text t } + + + + +-- add a new function to the current helper module +newLocalFunDecl :: CFunDecl -> [CStmt] -> M () +newLocalFunDecl d body = updHelper $ \m -> + m { fun_decls = (d, block body) : fun_decls m } + +-- add a new static variable to the current helper module +newLocalDecl :: CDecl -> M () +newLocalDecl d = updHelper $ \m -> + m { var_decls = static d : var_decls m } + + +cppFun :: String -> M () +cppFun t = cppFun' (text t) + +-- Add "#include function" to the current helper module +cppFun' :: Doc -> M () +cppFun' t = updHelper $ \m -> + m { cpp_funs = cpp_funs m $$ t } + + + +runM :: SamplerConf -> M a -> (a, S) +runM conf (M m) = runId $ runStateT start $ runReaderT info m + where start = S { main_mod = blankMod + , cur_mod = Nothing + , helper_mods = [] + , cnames = maxNode + 1 + } + info = R { config = conf, sampling = Nothing } + (maxNode,_) = IM.findMax $ stoNodes $ graph conf + + +renderMod :: CModule -> Doc +renderMod m = + cpp_stuff m $$ + char ' ' $$ text "/* Variable declarations */" $$ + decls (var_decls m) $$ + char ' ' $$ text "/* Function types */" $$ + decls (map fst (fun_decls m)) $$ + char ' ' $$ text "/* Included templates */" $$ + cpp_funs m $$ + char ' ' $$ text "/* Function definitions */" $$ + vcat [ d $$ b | (d,b) <- fun_decls m ] + + where decls = vcat . map (\d -> d <> semi) + + +renderState :: SamplerConf -> S -> [(FilePath, Doc)] +renderState conf s0 = ("sampler.h", hdr) + : ("sampler.c", main) + : map helper (helper_mods s) + where + s = if split_files conf then s0 else noHelpers s0 + mm = main_mod s + hdr = decls (map extern (var_decls mm)) + main = renderMod mm { + var_decls = concatMap extern_helper (helper_mods s) ++ + + var_decls mm } + helper (i,m) = + ( "slice_" ++ show i ++ ".c" + , renderMod $ m { cpp_stuff = text "#include \"passage.h\"" $$ + text "#include \"sampler.h\"" } + ) + + extern_helper (h,m) + | special_slicers conf + = [ text ("extern double SLICE(" ++ show h ++ ")(double)") + , text ("extern double SLICE_TUNE(" ++ show h ++ ")(double)") + ] + | otherwise = map (extern . fst) (fun_decls m) + + decls = vcat . map (\d -> d <> semi) + + + +-------------------------------------------------------------------------------- + +call_slicer :: (NodeIdx, StoVar) -> M ([CStmt], CExpr,CExpr) +call_slicer x = + do special <- M (asks (special_slicers . config)) + if special then call_special_slicer x else call_generic_slicer x + + + +call_generic_slicer :: (NodeIdx, StoVar) -> M ([CStmt], CExpr,CExpr) +call_generic_slicer (ix,sv) = + do v <- cnameVar (ix,sv) + case priSupport (stoVarPrior sv) of + Real -> + do newDecl $ var_decl double wid + let slice = var (ident "slice_real") + tune = var (ident "tune_slice_real") + return + ( [ assign (var wid) (double_lit 1) ] + , call slice [ llfun, var wid, v ] + , call tune [ llfun, addr_of wid, v ] + ) + + PosReal -> + do newDecl $ var_decl double wid + let z = int_lit 0 + slice = var (ident "slice_pos_real") + tune = var (ident "tune_slice_pos_real") + return + ( [ assign (var wid) (double_lit 1) ] + , call slice [ llfun, var wid, z, v ] + , call tune [ llfun, addr_of wid, z, v ] + ) + + Interval lo hi -> + do e1 <- term lo + e2 <- term hi + let slice = var (ident "slice_real_left_right") + expr = call slice [ llfun, e1, e2, v ] + return ([], expr, expr) + + Discrete (Just t) -> + do e <- term t + let slice = var (ident "slice_discrete_right") + expr = call slice [ llfun, e, v ] + return ([], expr, expr) + + Discrete Nothing -> + do let slice = var (ident "slice_discrete") + expr = call slice [ llfun, v ] + return ([], expr, expr) + + where llfun = var $ ident $ "LL_FUN(" ++ show ix ++ ")" + wid = ident $ "WIDTH(" ++ show ix ++ ")" + + +call_special_slicer :: (NodeIdx, StoVar) -> M ([CStmt], CExpr,CExpr) +call_special_slicer (ix,sv) = + do v <- cnameVar (ix,sv) + cppFun ("#define VAR " ++ show ix) + + let fun = var $ ident $ "SLICE(" ++ show ix ++ ")" + the_tune_fun = var $ ident $ "SLICE_TUNE(" ++ show ix ++ ")" + + tune_fun <- case priSupport (stoVarPrior sv) of + Real -> + do cppFun "#include \"templates/slice.c\"" + return the_tune_fun + + PosReal -> + do cppFun "#define LEFT 0" + cppFun "#include \"templates/slice.c\"" + cppFun "#undef LEFT" + return the_tune_fun + + Interval lo hi -> + do e1 <- term lo + e2 <- term hi + cppFun' (text "#define LEFT " <+> parens e1) + cppFun' (text "#define RIGHT " <+> parens e2) + cppFun "#include \"templates/slice.c\"" + cppFun "#undef LEFT" + cppFun "#undef RIGHT" + return fun + + Discrete (Just t) -> + do e <- term t + cppFun' (text "#define RIGHT" <+> parens e) + cppFun "#include \"templates/finiteMetropolis.c\"" + cppFun "#undef RIGHT" + return fun + + Discrete Nothing -> + do cppFun "#include \"templates/metropolis_posreal.c\"" + return fun + + cppFun "#undef VAR" + return ( [] + , call fun [v] + , call tune_fun [v] + ) + + + +initOrder :: [(NodeIdx, StoVar)] -> M [(NodeIdx, StoVar)] +initOrder ns = + do bg <- getGraph + return $ map check $ stronglyConnComp [ (n,ix,uses bg v) | n@(ix,v) <- ns ] + where + uses bg = IS.toList . fvsSupport (fvsArray bg) . priSupport . stoVarPrior + + check (AcyclicSCC d) = d + check (CyclicSCC _) = error "Cannot initialize: recursive support!" + +init_code :: (NodeIdx, StoVar) -> M CStmt +init_code (x,sv) = + do v <- cnameVar (x,sv) + i <- case priSupport (stoVarPrior sv) of + Real -> return $ double_lit 0 + Discrete _ -> return $ double_lit 0 + PosReal -> return $ double_lit 1 + Interval lo hi -> term (lo + (hi - lo) / 2) + -- duplicates lo but, hopefully, this does not matter too much + + return (assign v i) + +{- If an observed variable is stored in an array, then we need to initialize + the corresponding entry in the array with the observed value. The reason + for this is that there may be expressions of the form: a[i], with "a" + begin an observed array, and "i" which is not statically known. +-} +init_code_initialized :: (NodeIdx, Double) -> M [CStmt] +init_code_initialized (x,d) = + do sv <- lookupVar x + case stoVarName sv of + InArray {} -> + do v <- cnameVar (x,sv) + return [assign v (double_lit d)] + _ -> return [] + +-------------------------------------------------------------------------------- + +ll_summand :: (Term NodeIdx, Term NodeIdx) + -> M ([CStmt], Term NodeIdx, IS.IntSet) +ll_summand (x,c) = + do bg <- getGraph + if isSimpleTerm c + then return ([], x * c, varsOf bg (x * c)) + else do c1 <- newDetVar + let c2 = simpleCName c1 + newLocalDecl $ var_decl double c2 + expr <- term c + return ( [assign (var c2) expr] + , x * tvar c1 + , varsOf bg c `IS.union` varsOf bg x + ) + where + varsOf bg = leavesOfTerm (fvsArray bg) + +data StoVarCode = + StoVarCode + { tuneCode :: [CStmt] + , sliceCode :: [CStmt] + , locality :: (Int,[Int]) -- array number, indexes. + -- clobals would be a 1-dim array + -- call "-1". (XXX) + } + +sto_var :: (NodeIdx, StoVar) -> M ( [CStmt] -- init code + , (NodeIdx, StoVarCode, IS.IntSet) + ) +sto_var (ix,sv) = newHelper ix $ + do let xParam = simpleCName ix + iname = ident ("INIT_DET_VARS(" ++ show ix ++ ")") + llname = ident $ "LL_FUN(" ++ show ix ++ ")" + + -- Here we compute a "locality" for the variable. + -- This is useful when we group work by thread because + -- we prefer to put close updates together. + loc <- case stoVarName sv of + InArray aix ixes -> return (aix,ixes) + _ -> do newDecl $ var_decl double xParam + return (-1,[0]) -- XXX: count which vars are close + + (is,ts,vs) <- unzip3 `fmap` mapM ll_summand (M.toList (stoPostDistLL sv)) + expr <- nowSampling ix (term (sum ts)) + init_dets <- case concat is of + [] -> return [] + have_dets -> + do newLocalFunDecl (fun_decl void iname []) have_dets + return [ callS (var iname) [] ] + newLocalFunDecl (fun_decl double llname [(double,xParam)]) [ creturn expr ] + x <- cnameVar (ix,sv) + (initW, sliceExpr,sliceTuneExpr) <- call_slicer (ix,sv) + ic <- init_code (ix,sv) + return ( initW ++ [ic] -- initializaztion code + , ( ix + , StoVarCode + { tuneCode = init_dets ++ [ assign x sliceTuneExpr ] + , sliceCode = init_dets ++ [ assign x sliceExpr ] + , locality = loc + } + , IS.unions vs -- sto. var. deps. + ) + ) + + +-------------------------------------------------------------------------------- + +data SamplerConf = SamplerConf + { graph :: BayesianGraph + , sampleNum :: Int + , itsPerSample :: Int + , warmup :: Int + , thread_num :: Int + , seed :: [Int] + , monitor :: [(String, Term NodeIdx)] + , observe :: IM.IntMap Double -- Map node indexes to observed values. + , initialize :: IM.IntMap Double -- Map node indexes to initialized values. + , special_slicers :: Bool -- Generate a custom slicer per variable? + , split_files :: Bool -- Make one file per variable? + } + +declareArr :: (NodeIdx, ArrayInfo) -> M () +declareArr (ix,i) = + newDecl $ array_decl double (arrName (ix,i)) (map size (arrayDimensions i)) + where size (x,y) = y - x + 1 + +arrName :: (NodeIdx,ArrayInfo) -> CIdent +arrName (x,_) = ident ("a_" ++ show x) + +{- +genParGroups :: Int -> (StoVarCode -> [CStmt]) -> [[StoVarCode]] -> [CStmt] +genParGroups cpus which xs = concatMap makeSections xs + where + entries_per_thread len = (len + cpus - 1) `div` cpus + + makeSections vs = [ pragma "omp sections" + , block $ concatMap makeSection + $ chunks (entries_per_thread len) + $ sortBy (compare `on` locality) vs + ] + where len = length vs + + makeSection vs = [ pragma "omp section" + , block (concatMap which vs) + ] + +-} + +genParGroups :: Int -> (StoVarCode -> [CStmt]) -> [[StoVarCode]] -> [[CStmt]] +genParGroups cpus which = map concat . transpose . map threadBlocks + + where + entries_per_thread len = (len + cpus - 1) `div` cpus + + -- Allocate a list of indipendent statements to different threads. + -- Each blocks start with a barrier + threadBlocks vs = map makeBlock + $ addBlanks cpus + $ chunks (entries_per_thread len) + $ sortBy (compare `on` locality) vs + where len = length vs + + makeBlock vs = pragma "omp barrier" : concatMap which vs + + -- If there is not enough work for all threads, we insert + -- an empty list, so that we still get a barrier, otherwise + -- things get our of sync. + addBlanks n [] = replicate n [] + addBlanks n (x : xs) = x : seq m (addBlanks m xs) + where m = n - 1 + + + +-- Split a list into chunks of the given length. +-- If we run out of elements we make empty lists. +chunks :: Int -> [a] -> [[a]] +chunks n xs = case splitAt n xs of + (as,bs) -> as : case bs of + [] -> [] + _ -> chunks n bs + + + +genThread :: [(CStmt,CStmt)] -> Int -> ([CStmt], [CStmt]) -> M [CStmt] +genThread monitor_code n (tune_code,sample_code) = + do newFunDecl (fun_decl void name []) $ + [ var_decl unsigned_long (ident "i") <> semi + , var_decl unsigned_long (ident "j") <> semi + ] + + ++ + ifMaster ( map fst monitor_code ++ + [ nl, toStdErr [ string_lit "Tuning width parameters.\n" ] ] + ) + + ++ + [ text ("for (i = 0; i < warm_up_steps; ++i)") + $$ nest 2 (block tune_code) + , barrier + ] + + ++ + ifMaster [ toStdErr [string_lit "Sampling.\n"] ] + + ++ + [ text ("for (i = 0; i < number_of_samples; ++i)") $$ + nest 2 (block $ + ifMaster [ ppProg ] + + ++ + [ text ("for (j = 0; j < steps_per_sample; ++j)") + $$ nest 2 (block sample_code) + , barrier + ] + + ++ + ifMaster (printRowLabel : map snd monitor_code ++ [ nl ]) + ) + ] + return [ pragma "omp section" + , callS name [] + ] + + where + name = ident ("thread_" ++ show n) + ifMaster xs = if n == 0 then xs else [] + barrier = pragma "omp barrier" + toStdErr xs = callS (var (ident "fprintf")) (var (ident "stderr") : xs) + toStdOut xs = callS (var (ident "printf")) xs + ppProg = callS (ident "progress") [ var (ident "i") ] + + printRowLabel = toStdOut [ string_lit "%lu", ident "i" ] + nl = toStdOut [ string_lit "\n" ] + + + + + + + +gen_c :: SamplerConf -> [(FilePath, Doc)] +gen_c conf = renderState conf $ snd $ runM conf $ + do let bg = graph conf + cpp "#include " + cpp "#include " + cpp "#include " + cpp "#include \"passage.h\"" + + mapM_ declareArr $ IM.toList $ stoArryas bg + + -- We generate sampling code only for stochastic variables that + -- are not observed: + let observedVars = IM.keysSet (observe conf) + sampledNodes = filter (not . (`IS.member` observedVars) . fst) + $ IM.toList $ stoNodes bg + (ins,deps) <- unzip `fmap` (mapM sto_var =<< initOrder sampledNodes) + let dropObserved (x,y,zs) = + (x,y, IS.filter (not . (`IS.member` observedVars)) zs) + par_groups = map (map snd) $ groupByColor $ map dropObserved deps + + -- Observed stochastic variables just get initialized once. + -- Variables that are not in an array don't even need to be initialized + -- but it is important the we initialize arrays, because of expressions + -- a[i], where "a" is observed but "i" is not. + obs_ins <- mapM init_code_initialized $ IM.toList $ observe conf + + init_ins <- mapM init_code_initialized $ IM.toList $ initialize conf + + let cpus = thread_num conf + + newFunDecl (fun_decl void (ident "set_defaults") []) $ + [ assign (var (ident "number_of_samples")) $ int_lit $ sampleNum conf + , assign (var (ident "steps_per_sample")) $ int_lit $ itsPerSample conf + , assign (var (ident "warm_up_steps")) $ int_lit $ warmup conf + , assign (var (ident "num_threads")) $ int_lit cpus + , assign (var (ident "have_seed")) $ int_lit $ length $ seed conf + ] ++ [ assign (arr_ix (var (ident "seeds")) (int_lit n)) (int_lit v) + | (n,v) <- zip [ 0 .. ] (reverse (seed conf)) ] + + newFunDecl (fun_decl void (ident "init_vars") []) + $ concat $ obs_ins ++ ins + + -- Code to print the values of monitored expressions. + monitor_code <- forM (monitor conf) $ \(lab,x) -> + do expr <- term x + return ( callS (ident "printf") [ string_lit ("\t" ++ lab) ] + , callS (ident "printf") [ string_lit ("\t%f"), expr ] + ) + + let tune_codes = genParGroups cpus tuneCode par_groups + slice_codes = genParGroups cpus sliceCode par_groups + + threads <- zipWithM (genThread monitor_code) + [ 0 .. ] (zip tune_codes slice_codes) + + + -- The main sampling function. + newFunDecl (fun_decl void (ident "sampler") []) $ + [ pragma "omp sections" + , block (concat threads) + ] + diff --git a/Language/Passage/GraphColor.hs b/Language/Passage/GraphColor.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/GraphColor.hs @@ -0,0 +1,53 @@ +-- | A simple graph coloring algorithm. +module Language.Passage.GraphColor (groupByColor) where + +import qualified Data.IntMap as IM +import qualified Data.IntSet as IS +import Data.List (sortBy, groupBy) +import Data.Foldable(foldl') +import Data.Maybe (mapMaybe) +import Data.Function (on) + +type Color = Int +type Coloring a = IM.IntMap (a,Color) + +choose :: Coloring a -> (Int, (a, IS.IntSet)) -> Coloring a +choose coloring (key,(a,ns)) = IM.insert key (a, head candidates) coloring + where used = map snd $ mapMaybe (`IM.lookup` coloring ) $ IS.toList ns + candidates = [ x | x <- [ 0 .. ], not (x `elem` used) ] + +addNode :: IM.IntMap (a,IS.IntSet) + -> (Int, a, IS.IntSet) + -> IM.IntMap (a, IS.IntSet) +addNode g0 (x,a,xs) = IS.fold addBack (IM.insertWith jnUseNew x node g0) xs + where + addBack y g = IM.insertWith jnUseOld y (missing y,IS.singleton x) g + node = (a,xs) + jnUseOld (_,vs1) (b,vs2) = (b, IS.union vs1 vs2) + jnUseNew (_,vs1) (_,vs2) = (a, IS.union vs1 vs2) + missing y = error ("BUG: Variable " ++ show y ++ " is missing?") + +colorG :: [(Int,a,IS.IntSet)] -> Coloring a +colorG = foldl' choose IM.empty + . sortBy earlier + . IM.toList + . foldl' addNode IM.empty + + where earlier (_,(_,x)) (_,(_,y)) = compare (IS.size y) (IS.size x) + + +colorGroups :: Coloring a -> [[(Int,a)]] +colorGroups = map (map snd) + . groupBy ((==) `on` fst) + . sortBy (compare `on` fst) + . map swap + . IM.toList + + where swap (x,(a,c)) = (c,(x,a)) + + +groupByColor :: [(Int,a, IS.IntSet)] -> [[(Int,a)]] +groupByColor = colorGroups . colorG + + + diff --git a/Language/Passage/Lang/C.hs b/Language/Passage/Lang/C.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/Lang/C.hs @@ -0,0 +1,121 @@ +module Language.Passage.Lang.C where + +import Language.Passage.Utils hiding (double,int) +import Data.Word + +-------------------------------------------------------------------------------- +-- Generating C +-------------------------------------------------------------------------------- + +type CExpr = Doc +type CStmt = Doc +type CType = Doc +type CDecl = Doc +type CFunDecl = Doc +type CIdent = Doc + + +ident :: String -> CIdent +ident = text + +call :: CExpr -> [CExpr] -> CExpr +call f xs = f <> commaSep xs + + +cast :: CType -> CExpr -> CExpr +cast t e = parens t <+> parens e + +arr_ix :: CExpr -> CExpr -> CExpr +arr_ix a i = a <> brackets i + +var :: CIdent -> CExpr +var x = x + +double_lit :: Double -> CExpr +double_lit = text . show + +int_lit :: Int -> CExpr +int_lit = pp + +unsigned_lit :: Word -> CExpr +unsigned_lit = pp + +string_lit :: String -> CExpr +string_lit = text . show + +-- could be an array, in general +addr_of :: CIdent -> CExpr +addr_of x = char '&' <> x + +main_name :: CIdent +main_name = ident "main" + + +double :: CType +double = text "double" + +void :: CType +void = text "void" + +int :: CType +int = text "int" + +unsigned :: CType +unsigned = text "unsigned" + +unsigned_long :: CType +unsigned_long = text "unsigned long" + +array_decl :: CType -> CIdent -> [Int] -> CDecl +array_decl t a ds = t <+> a <> hcat (map (brackets . pp) ds) + +extern :: CDecl -> CDecl +extern d = text "extern" <+> d + +static :: CDecl -> CDecl +static d = text "static" <+> d + +static_fun :: CFunDecl -> CDecl +static_fun d = text "static" $$ d + +fun_decl :: CType -> CIdent -> [(CType, CIdent)] -> CFunDecl +fun_decl res f args = res <+> f <> args_decl + where + args_decl = case args of + [] -> parens (text "void") + _ -> commaSep [ t <+> x | (t, x) <- args ] + +var_decl :: CType -> CIdent -> CDecl +var_decl t x = t <+> x + +-- In C this is an expr, but we don't need this "flexibility". +assign :: CExpr -> CExpr -> CStmt +assign x y = x <+> text "=" <+> y <> semi + +creturn :: CExpr -> CStmt +creturn x = text "return" <+> x <> semi + +switch :: CExpr -> [(Int,[CStmt])] -> CStmt -> CStmt +switch e as d = + text "switch" <+> parens e <+> char '{' + $$ nest 2 ( vcat (map ppCase as) + $$ text "default:" <+> d) + $$ char '}' + where ppCase (i,s) = text "case" <+> int_lit i <> char ':' <+> vcat s + +cbreak :: CStmt +cbreak = text "break" <> semi + + +callS :: CExpr -> [CExpr] -> CStmt +callS f xs = call f xs <> semi + +block :: [CStmt] -> CStmt +block ss = char '{' $+$ nest 2 (vcat ss) $+$ char '}' + +declBlock :: [CDecl] -> [CStmt] -> CStmt +declBlock ds ss = char '{' $+$ nest 2 (vcat ds $+$ vcat ss) $+$ char '}' + +pragma :: String -> CStmt +pragma x = text "#pragma" <+> text x + diff --git a/Language/Passage/Lang/LaTeX.hs b/Language/Passage/Lang/LaTeX.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/Lang/LaTeX.hs @@ -0,0 +1,70 @@ +module Language.Passage.Lang.LaTeX where + +import Language.Passage.Utils +import qualified Data.Set as S + + +class LaTeX t where + latex :: t -> Doc + +instance LaTeX Double where + latex = double + +instance LaTeX Int where + latex = int + +cmd :: String -> [Doc] -> Doc +cmd a as = char '\\' <> text a <> case as of + [] -> empty + _ -> hcat (map braces as) + +frac, pow :: Doc -> Doc -> Doc +frac a b = cmd "frac" [ a, b ] +pow a b = braces a <> char '^' <> braces b + +lg :: Doc +lg = cmd "log" [] + +expon = cmd "exp" [] + +logGamma :: Doc -> Doc +logGamma a = cmd "log" [cmd "Gamma" [a]] + +literal :: String -> Doc -> Doc +literal s a = cmd s [a] + +knownVars :: S.Set String +knownVars = S.fromList + [ "alpha", "beta", "gamma", "delta", "epsilon", "varepsilon" + , "zeta", "eta", "theta", "vartheta", "kappa", "lambda", "mu", "nu", "xi" + , "pi", "varpi", "rho", "varrho", "sigma", "varsigma" + , "tau", "upsilon", "phi", "varphi", "chi", "psi", "imega" + ] + +var :: String -> Doc +var x | x `S.member` knownVars = cmd x [] + | short x = text x + | True = cmd "ensuremath" [ cmd "mathit" [ text x ]] + where short [_] = True + short _ = False + +math :: Doc -> Doc +math x = char '$' <> x <> char '$' + +env :: String -> [Doc] -> Doc -> Doc +env x opts body = cmd "begin" (text x : opts) $$ body $$ cmd "end" [text x] + +sim :: Doc +sim = cmd "sim" [] + +row :: [Doc] -> Doc +row xs = hsep (punctuate (text "&") xs) <+> text "\\\\" + +hline :: Doc +hline = cmd "hline" [] + +mathcal :: Doc -> Doc +mathcal d = cmd "mathcal" [d] + +mathrm :: Doc -> Doc +mathrm d = cmd "mathrm" [d] diff --git a/Language/Passage/SimulatorConf.hs b/Language/Passage/SimulatorConf.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/SimulatorConf.hs @@ -0,0 +1,143 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +module Language.Passage.SimulatorConf where + +import MonadLib +import qualified Data.IntMap as IM +import qualified Data.IntSet as IS +import Language.Passage.Term +import Language.Passage.AST +import Language.Passage.Utils +import Language.Passage.Graph + +-- The "simulator" infrastructure +data SimState = SimState + { cfgSampleNum :: Int + , cfgItsPerSample :: Int + , cfgWarmup :: Int + , cfgMersenne :: Bool + , cfgProfile :: Bool + , cfgMonitor :: [(String, Term NodeIdx)] + , cfgObserve :: IM.IntMap Double -- NodeIdx |-> observed values + , cfgInitialize :: IM.IntMap Double + , cfgRandomSeed :: [Int] -- reversed list of seeds for threads + , cfgNetwork :: Maybe BayesianGraph + , cfgThreadNum :: Int + , cfgSpecialSlicers :: Bool + , cfgSplitFiles :: Bool + } + +initSimState :: SimState +initSimState = SimState { cfgSampleNum = 100000 + , cfgItsPerSample = 10 + , cfgWarmup = 1000 + , cfgMersenne = False + , cfgProfile = False + , cfgMonitor = [] + , cfgObserve = IM.empty + , cfgInitialize = IM.empty + , cfgRandomSeed = [] + , cfgNetwork = Nothing + , cfgThreadNum = 2 + , cfgSpecialSlicers = False + , cfgSplitFiles = False + } + + +newtype BayesianSimulator a = B (StateT SimState Id a) + deriving (Functor, Monad) + +upd :: (SimState -> SimState) -> BayesianSimulator () +upd f = B $ sets_ f + +getField :: (SimState -> a) -> BayesianSimulator a +getField f = B (f `fmap` get) + +runSim :: BayesianSimulator a -> SimState +runSim (B m) = snd $ runId $ runStateT initSimState m + +setWarmupCount :: Int -> BayesianSimulator () +setWarmupCount i = upd $ \s -> s { cfgWarmup = i } + +setSampleCount :: Int -> BayesianSimulator () +setSampleCount i = upd $ \s -> s { cfgSampleNum = i } + +setIterationsPerSample :: Int -> BayesianSimulator () +setIterationsPerSample i = upd $ \s -> s { cfgItsPerSample = i } + +-- | Set the random seed for a thread. This function may be calledd +-- multiple times to set the seeds for multiple threads. +-- The seeds are used in order: first call is for thread 0, next for thread 1, +-- etc. +setRandomSeed :: Int -> BayesianSimulator () +setRandomSeed d = upd $ \s -> s { cfgRandomSeed = d : cfgRandomSeed s } + +setThreadNum :: Int -> BayesianSimulator () +setThreadNum n = upd $ \s -> s { cfgThreadNum = n } + +useMersenneTwister :: Bool -> BayesianSimulator () +useMersenneTwister b = upd (\s -> s { cfgMersenne = b }) + +-- | When using a specialized slizer, we generate a custom slicer +-- for each stochastic variable. The benefit of this is that, in principle, +-- this may result in more efficient code, at the cost of longer compilation +-- time, and larger binary. The alternative is to use a generic slicing +-- function which is parameterized by the log-likelihood function for +-- a variable. +useSpecialSlicers :: Bool -> BayesianSimulator () +useSpecialSlicers b = upd (\s -> s { cfgSpecialSlicers = b }) + +-- | Generate a separate file for each stochastic variable. +-- The benefit of this flag is that it makes it possible to compile +-- multiple files in parallel. The drawback is that some optimizations +-- may be lost because the files are compiled separately. Also, there +-- is some overhead for processing multiple files. +splitFiles :: Bool -> BayesianSimulator () +splitFiles b = upd (\s -> s { cfgSplitFiles = b }) + +enableProfiling :: Bool -> BayesianSimulator () +enableProfiling b = upd $ \s -> s { cfgProfile = b } + +model :: BayesianNetwork a -> BayesianSimulator a +model t = do upd $ \s -> s { cfgNetwork = Just nw } + return a + where (a, nw) = buildBayesianGraph t + +observe :: Term NodeIdx -> Double -> BayesianSimulator () +observe t d = + case splitArray t of + -- (TVar idx, its) | Just is <- mapM isConst its + (TVar idx, []) -> + do obs <- getField cfgObserve + case IM.insertLookupWithKey unused idx d obs of + (Nothing, m1) -> upd (\s -> s { cfgObserve = m1 }) + (Just _, _) -> + error $ "observe: Model error. Node was observed before: " + ++ show t + + _ -> error $ "observe: Model error. Only nodes can be observed, received: " + ++ show t + + where unused = error "BUG: observe--not used" + +initialize :: Term NodeIdx -> Double -> BayesianSimulator () +initialize t d = + case splitArray t of + -- (TVar idx, its) | Just is <- mapM isConst its + (TVar idx, []) -> + do init <- getField cfgInitialize + case IM.insertLookupWithKey unused idx d init of + (Nothing, m1) -> upd (\s -> s { cfgInitialize = m1 }) + (Just _, _) -> + error $ "initialize: Model error. Node was initialized before: " + ++ show t + + _ -> error $ "initialize: Model error. Only nodes can be initialized, received: " + ++ show t + + where unused = error "BUG: initialize--not used" + +monitor :: String -> Expr -> BayesianSimulator () +monitor nm e = upd (\s -> s { cfgMonitor = (nm, e) : cfgMonitor s }) + + + diff --git a/Language/Passage/SliceSample.hs b/Language/Passage/SliceSample.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/SliceSample.hs @@ -0,0 +1,110 @@ +module Language.Passage.SliceSample where + +import System.Random + +{- GNUPLOT commands to see the distribution: +binwidth=5 +bin(x,width)=width*floor(x/width) + +plot 'datafile' using (bin($1,binwidth)):(1.0) smooth freq with boxes +-} + +{- +Some R code for testing: + +db <- read.table('data.beta', header=F)[,1] +dg <- read.table('data.gamma', header=F)[,1] +dn <- read.table('data.norm', header=F)[,1] + +library(MASS) + +s <- seq(0,5,length=1000) + +truehist(db) +lines(s, dbeta(s,2,5),col='red') + +truehist(dg) +lines(s, dgamma(s,2,5),col='red') + +qqnorm(dn) + +-} + +genAll :: IO () +genAll = mapM_ (genTest 100000) [0..3] + +genTest :: Integer -> Int -> IO () +genTest cnt i = do + putStr $ "Generating " ++ show nm ++ ". " + xs <- genPts cnt slicer 1 5 source + writeFile nm $ unlines (map show xs) + putStrLn "Done." + where gammaLL a b x = (a - 1) * log x - b * x + normLL = (* 0.5) . negate . (**2) + betaLL a b x = (a - 1) * log x + (b-1) * log (1 - x) + (nm, slicer, source) = + [ ("data.gamma", slicePos, gammaLL 2 5) + , ("data.norm", slice, normLL) + , ("data.beta", sliceUnit, betaLL 2 5) + , ("data.truncNormal",sliceUnit, normLL) + ] !! i + +genPts :: Integer -> (StdGen -> Double -> Double -> (Double -> Double) -> (Double, Double, StdGen)) + -> Double -> Double -> (Double -> Double) -> IO [Double] +genPts cnt slicer w initX ll = do + g <- newStdGen + return $ go g initX cnt + where go _ _ 0 = [] + go g0 x0 c = let (x1, _, g1) = slicer g0 w x0 ll + in x1 : go g1 x1 (c-1) + +sliceUnit :: StdGen -- ^ Source of randomness + -> Double -- ^ Width of initial sampling interval + -> Double -- ^ variable + -> (Double -> Double) -- ^ Log likelyhood, in terms of varibale + -> (Double, Double, StdGen) -- ^ New value for variable, together with its log-likelihood +sliceUnit g w x0 ll = genericSlice g w x0 ll (\_ _ -> 0) (\_ _ -> 1) + +slicePos :: StdGen -- ^ Source of randomness + -> Double -- ^ Width of initial sampling interval + -> Double -- ^ variable + -> (Double -> Double) -- ^ Log likelyhood, in terms of varibale + -> (Double, Double, StdGen) -- ^ New value for variable, together with its log-likelihood +slicePos g w x0 ll = genericSlice g w x0 ll (\_ _ -> 0) (\y lo -> search y (right_pts y lo)) + where search y = head . dropWhile (\p -> ll p > y) + right_pts _ lo = [ right, right + w .. ] + where right = x0 + lo + +slice :: StdGen -- ^ Source of randomness + -> Double -- ^ Width of initial sampling interval + -> Double -- ^ variable + -> (Double -> Double) -- ^ Log likelyhood, in terms of varibale + -> (Double, Double, StdGen) -- ^ New value for variable, together with its log-likelihood +slice g w x0 ll = genericSlice g w x0 ll (\y lo -> search y (left_pts y lo)) (\y lo -> search y (right_pts y lo)) + where search y = head . dropWhile (\p -> ll p > y) + left_pts _ lo = [ left, left - w .. ] + where left = x0 - lo + right_pts _ lo = [ right, right + w .. ] + where right = x0 - lo + w + +genericSlice :: StdGen -- ^ Source of randomness + -> Double -- ^ Width of initial sampling interval + -> Double -- ^ variable + -> (Double -> Double) -- ^ Log likelyhood, in terms of varibale + -> (Double -> Double -> Double) -- ^ left bound + -> (Double -> Double -> Double) -- ^ right bound + -> (Double, Double, StdGen) -- ^ New value for variable, together with its log-likelihood +genericSlice g w x0 ll left right = + let (r, g1) = randomR (0, 1) g + y = log r + ll x0 + (lo, g2) = randomR (0, w) g1 + in pickRandom g2 x0 y ll (left y lo) (right y lo) + + +pickRandom :: StdGen -> Double -> Double -> (Double -> Double) -> Double -> Double -> (Double, Double, StdGen) +pickRandom g0 x0 y ll = go g0 + where go g l r = if ll_x1 < y then if x1 < x0 then go g1 x1 r + else go g1 l x1 + else (x1, ll_x1, g1) + where (x1, g1) = randomR (l, r) g + ll_x1 = ll x1 diff --git a/Language/Passage/SliceSampleMWC.hs b/Language/Passage/SliceSampleMWC.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/SliceSampleMWC.hs @@ -0,0 +1,119 @@ +module Language.Passage.SliceSampleMWC where + +import System.Random.MWC +import Control.Monad.Primitive + +{- GNUPLOT commands to see the distribution: +binwidth=5 +bin(x,width)=width*floor(x/width) + +plot 'datafile' using (bin($1,binwidth)):(1.0) smooth freq with boxes +-} + +{- +Some R code for testing: + +db <- read.table('data.beta', header=F)[,1] +dg <- read.table('data.gamma', header=F)[,1] +dn <- read.table('data.norm', header=F)[,1] + +library(MASS) + +s <- seq(0,5,length=1000) + +truehist(db) +lines(s, dbeta(s,2,5),col='red') + +truehist(dg) +lines(s, dgamma(s,2,5),col='red') + +qqnorm(dn) + +-} + +genAll :: IO () +genAll = mapM_ (genTest 100000) [0..3] + +genTest :: Integer -> Int -> IO () +genTest cnt i = do + putStr $ "Generating " ++ show nm ++ ". " + xs <- genPts cnt slicer 1 5 source + writeFile nm $ unlines (map show xs) + putStrLn "Done." + where gammaLL a b x = (a - 1) * log x - b * x + normLL = (* 0.5) . negate . (**2) + betaLL a b x = (a - 1) * log x + (b-1) * log (1 - x) + (nm, slicer, source) = + [ ("data.gamma", slicePos, gammaLL 2 5) + , ("data.norm", slice, normLL) + , ("data.beta", sliceUnit, betaLL 2 5) + , ("data.truncNormal",sliceUnit, normLL) + ] !! i +type Slicer = Gen RealWorld + -> Double + -> Double + -> (Double -> Double) + -> IO (Double, Double) + +genPts :: Integer -> Slicer -> Double -> Double -> (Double -> Double) -> IO [Double] +genPts cnt slicer w initX ll = do + g <- create + let + go _ 0 = return [] + go x0 c = do + (x1,_) <- slicer g w x0 ll + xs <- go x1 (c-1) + return (x1:xs) + go initX cnt + +sliceUnit :: Slicer +sliceUnit g w x0 ll = genericSlice g w x0 ll (\_ _ -> 0) (\_ _ -> 1) + +slicePos :: Slicer +slicePos g w x0 ll = genericSlice g w x0 ll (\_ _ -> 0) (\y lo -> search y (right_pts y lo)) + where + search y = head . dropWhile (\p -> ll p > y) + right_pts _ lo = [ right, right + w .. ] + where right = x0 + lo + +slice :: Slicer +slice g w x0 ll = genericSlice g w x0 ll (\y lo -> search y (left_pts y lo)) (\y lo -> search y (right_pts y lo)) + where search y = head . dropWhile (\p -> ll p > y) + left_pts _ lo = [ left, left - w .. ] + where left = x0 - lo + right_pts _ lo = [ right, right + w .. ] + where right = x0 - lo + w + +genericSlice :: (PrimMonad m) => + Gen (PrimState m) + -> Double + -> Double + -> (Double -> Double) + -> (Double -> Double -> Double) + -> (Double -> Double -> Double) + -> m (Double, Double) +genericSlice g w x0 ll left right = do + r <- uniform g + let y = log r + ll x0 + lo <- uniformR (0,w) g + pickRandom g x0 y ll (left y lo) (right y lo) + + +pickRandom + :: (Control.Monad.Primitive.PrimMonad m) => + Gen (Control.Monad.Primitive.PrimState m) + -> Double + -> Double + -> (Double -> Double) + -> Double + -> Double + -> m (Double, Double) + +pickRandom g x0 y ll = go + where + go l r = do + x1 <- uniformR (l, r) g + let ll_x1 = ll x1 + if ll_x1 < y then + if x1 < x0 then go x1 r else go l x1 + else return (x1, ll_x1) diff --git a/Language/Passage/Term.hs b/Language/Passage/Term.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/Term.hs @@ -0,0 +1,443 @@ +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE DeriveFunctor #-} + +module Language.Passage.Term where + +import Control.Monad(mplus) +import qualified Data.IntSet as IS +import Data.Ratio(numerator,denominator) +import Data.Maybe(fromMaybe) + +import Language.Passage.Utils +import qualified Language.Passage.Lang.LaTeX as LaTeX +import Language.Passage.Lang.LaTeX(LaTeX(..)) + + +data Op = TLog | TNeg | TAdd | TMul | TSub | TDiv | TPow | TLogGamma | TExp + | TCase -- for "arrays" of det. nodes. 1st arg index, rest array. + | TIx + | TLit String + deriving (Eq, Show, Ord) + +-- | A term in a stochastic context +data Term a = TVar a + | TArr a -- A node corresponding to an array + | TConst Double + | TApp Op [Term a] + deriving (Eq, Ord, Show, Functor) + +{- +sizeOf :: Term a -> Int +sizeOf (TVar{}) = 1 +sizeOf (TConst{}) = 0 +sizeOf (TApp _ xs) = 1 + sum (map sizeOf xs) +-} + +type ArrVars = NodeIdx -> IS.IntSet + +-- | Nodes hanging off of a term +leavesOfTerm :: ArrVars -> Term NodeIdx -> IS.IntSet +leavesOfTerm _ (TVar b) = IS.singleton b +leavesOfTerm arr (TArr b) = arr b +leavesOfTerm _ (TConst{}) = IS.empty +leavesOfTerm arr (TApp _ ts) = IS.unions (map (leavesOfTerm arr) ts) + +-- | Returns 'True' for terms that are not applications. +isSimpleTerm :: Term a -> Bool +isSimpleTerm t = + case t of + TApp _ _ -> False + TArr _ -> True + TVar _ -> True + TConst _ -> True + +tcase :: Term a -> [Term a] -> Term a +tcase e es = TApp TCase (e:es) + +tvar :: a -> Term a +tvar = TVar + +tarr :: a -> Term a +tarr = TArr + +tconst :: Double -> Term a +tconst = TConst + +isConst :: Term a -> Maybe Double +isConst (TConst d) = Just d +isConst _ = Nothing + +un :: Op -> Term a -> Term a +un op x = TApp op [x] + +bin :: Op -> Term a -> Term a -> Term a +bin op x y = TApp op [x,y] + + +termOp :: Term a -> Maybe Op +termOp (TApp op _) = Just op +termOp _ = Nothing + +logGamma :: Term a -> Term a +logGamma t = case t of + TConst a -> TConst (lgg a) + _ -> un TLogGamma t + where lgg :: Double -> Double + lgg z = 0.5 * (log (2*pi) - log z) + z*(log (z+1/(12*z-0.1/z)) - 1) + + + +tIx :: Term a -> Term a -> Term a +tIx a i = TApp TIx [a,i] + +-- Split a term into an arrya and indexes. +-- Example: a[1][2][3] ---> (a,[1,2,3]) +splitArray :: Term a -> (Term a, [Term a]) +splitArray t0 = loop t0 [] + where loop (TApp TIx [a,i]) is = loop a (i:is) + loop t is = (t,is) + + + +-------------------------------------------------------------------------------- +-- Smarter printing infrastructure + +precedence :: Op -> (Fixity, Rational) +precedence op = + case op of + TExp -> (Prefix, 100) + TLog -> (Prefix, 100) + TLogGamma -> (Prefix, 100) + TCase -> (Prefix, 100) + TNeg -> (Prefix, 100) + TAdd -> (Infix ToLeft, 6) + TSub -> (Infix ToLeft, 6) + TMul -> (Infix ToLeft, 7) + TDiv -> (Infix ToLeft, 7) + TPow -> (Infix ToRight, 8) + TIx -> (Infix ToLeft, 9) + TLit s -> (Prefix, 100) + +instance PP Op where + pp op = + case op of + TCase -> text "choose" + TExp -> text "exp" + TLogGamma -> text "logGamma" + TLog -> text "log" + TNeg -> text "-" + TAdd -> text "+" + TSub -> text "-" + TMul -> text "*" + TDiv -> text "/" + TPow -> text "^" + TIx -> text "!" + TLit s -> text s + +wrapLatex :: Op -> Posn -> Maybe Op -> Doc -> Doc +wrapLatex _ _ Nothing doc = doc +wrapLatex op pos (Just op1) doc = if shouldWrap then wrap else doc + where + shouldWrap = + case (pos,op) of + (ToLeft, TPow) -> True + (_, TPow) -> False + + (_, TAdd) -> op1 == TIx + (ToRight, TIx) -> op1 == TIx + (_, TIx) -> False + + (ToLeft, TSub) -> False + (_, TSub) -> op1 == TSub || op1 == TAdd || op1 == TNeg + || op1 == TIx + + (ToLeft, TMul) -> op1 == TSub || op1 == TAdd + (_, TMul) -> op1 == TSub || op1 == TAdd || op1 == TNeg + || op1 == TIx + + (_, TDiv) -> False + + (_, TLog) -> op1 == TAdd || op1 == TSub || op == TMul + (_, TExp) -> op1 == TAdd || op1 == TSub || op == TMul + (_, TLogGamma) -> op1 == TAdd || op1 == TSub || op == TMul + (_, TNeg ) -> op1 == TAdd || op1 == TSub || op == TMul + || op1 == TIx + (_, TCase) -> False + + wrap = text "\\left(" <> doc <> text "\\right)" + +instance LaTeX a => LaTeX (Term a) where + + latex (TApp op ts) = + case op of + TAdd -> dL <+> char '+' <+> dR + TSub -> dL <+> char '-' <+> dR + TMul -> dL <+> dR + TDiv -> LaTeX.frac dL dR + TPow -> LaTeX.pow dL dR + TLog -> LaTeX.lg <+> dL + TExp -> LaTeX.expon <+> dL + TLogGamma -> LaTeX.logGamma dL + TNeg -> char '-' <> dL + TIx -> dL <+> char '!' <+> dR -- XXX: This could be nicer + TCase -> let a : as = map latex ts + in commaSep as <> char '_' <> braces a + + TLit s -> (LaTeX.literal s) dL + where ds = zipWith pr (ToLeft : ToRight : repeat None) ts + pr pos t = wrapLatex op pos (termOp t) (latex t) + dL : ds1 = ds + dR : _ = ds1 + + latex (TVar x) = latex x + latex (TArr x) = latex x + latex (TConst a) = double a + + + +ppTerm :: (PP a) => (Posn,Rational) -> Term a -> Doc +ppTerm (pos,prec) (TApp op [l,r]) + | (Infix dir, myprec) <- precedence op = + let this = ppTerm (ToLeft, myprec) l <+> pp op <+> + ppTerm (ToRight, myprec) r + in if myprec > prec || (myprec == prec && pos == dir) + then this + else parens this + +ppTerm (_,p) (TApp op ts) = + let this = pp op <+> commaSep [ ppTerm (None,0) t | t <- ts ] + in if snd (precedence op) > p then this else parens this + +ppTerm (_,n) (TArr x) = text "!" <> ppPrec n x +ppTerm (_,n) (TVar x) = text "?" <> ppPrec n x +ppTerm _ (TConst a) = double a + +instance PP a => PP (Term a) where + ppPrec n = ppTerm (None,n) + +liftTerm1 :: (Term a -> Term a) -> (Double -> Double) -> Term a -> Term a +liftTerm1 _ c (TConst a) = TConst (c a) +liftTerm1 s _ a = s a + +liftTerm2 :: (Term a -> Term a -> Term a) -> (Double -> Double -> Double) -> Term a -> Term a -> Term a +liftTerm2 _ c (TConst a) (TConst b) = TConst (a `c` b) +liftTerm2 s _ a b = a `s` b + +tbd1 :: Show a => String -> Term a -> b +tbd1 w x = tbd ("Term." ++ w ++ " " ++ show x) + +tbd2 :: Show a => String -> Term a -> Term a -> b +tbd2 w x y = tbd ("Term." ++ w ++ show (x, y)) + +instance (Eq a, Show a) => Num (Term a) where + -- addition + TConst 0 + b = b + a + TConst 0 = a + TConst x + TConst y = TConst (x + y) + a + (TApp TAdd [b,c]) = (a + b) + c + a + TApp TNeg [b] = a - b + TApp TNeg [a] + b = b - a + TApp TDiv [a, x] + TApp TDiv [b, y] + | x == y = (a+b) / x + a + b = liftTerm2 (bin TAdd) (+) a b +{- + a + b + | sizeOf a < sizeOf b = walkAdd a b + | True = walkAdd b a +-} + + -- multiplication + TConst 0 * _ = TConst 0 + TConst 1 * b = b + TConst (-1) * b = negate b + TConst x * TConst y = TConst (x * y) + TConst x * TApp TDiv [TConst y, z] + = TConst (x*y) / z + TConst x * TApp TDiv [z, TConst y] + = TConst (x/y) * z + a * b@(TConst _) = b * a -- constants float left + TApp TNeg [a] * b = negate (a * b) + a * TApp TNeg [b] = negate (a * b) + a * b = liftTerm2 (bin TMul) (*) a b + + -- subtraction + a - TConst 0 = a + a - TApp TNeg [b] = a + b + TApp TNeg [a] - b = negate (a + b) + TApp TSub [b, c] - d = TApp TSub [b, c+d] + a - b = liftTerm2 (bin TSub) (-) a b + + -- negation + negate (TApp TNeg [x]) = x + negate x = liftTerm1 (un TNeg) negate x + + -- others + abs = liftTerm1 (tbd1 "abs") abs + signum = liftTerm1 (tbd1 "signum") signum + fromInteger = TConst . fromInteger + +instance (Eq a, Show a) => Fractional (Term a) where + -- division + a / TConst 1 = a + TConst x / TConst y | y /= 0 = TConst (x / y) + (TApp TDiv [TConst c1, x]) / TConst c2 = TConst (c1/c2) / x + a / b = liftTerm2 (bin TDiv) (/) a b + + recip x = 1 / x + fromRational x = fromInteger (numerator x) / fromInteger (denominator x) + +instance (Eq a, Show a) => Floating (Term a) where + pi = TConst pi + exp = liftTerm1 (un TExp) exp + sqrt = liftTerm1 (tbd1 "sqrt") sqrt + log = liftTerm1 (un TLog) log + + -- power + _ ** TConst 0 = TConst 1 + TConst 0 ** _ = TConst 0 + a ** b = liftTerm2 (bin TPow) (**) a b + + -- TBD: Add support for these as needed + logBase = liftTerm2 (tbd2 "logBase") logBase + sin = liftTerm1 (tbd1 "sin") sin + tan = liftTerm1 (tbd1 "tan") tan + cos = liftTerm1 (tbd1 "cos") cos + asin = liftTerm1 (tbd1 "asin") asin + atan = liftTerm1 (tbd1 "atan") atan + acos = liftTerm1 (tbd1 "acos") acos + sinh = liftTerm1 (tbd1 "sinh") sinh + tanh = liftTerm1 (tbd1 "tanh") tanh + cosh = liftTerm1 (tbd1 "cosh") cosh + asinh = liftTerm1 (tbd1 "asinh") asinh + atanh = liftTerm1 (tbd1 "atanh") atanh + acosh = liftTerm1 (tbd1 "acosh") acosh + +-- Add two terms symbolically, trying to simplify as much as possible. +-- TODO: This process is currently quite ad-hoc, due to the lack of a +-- well defined normal form for terms. Also, developers will likely to +-- be able to add their own rules if necessary.. +-- +-- sAdd returns Nothing if it did Nothing, and Just t, if it was able +-- to do something interesting. +sAdd :: (Show t, Eq t) => Term t -> Term t -> Maybe (Term t) + +-- x + x == 2x +sAdd x y + | x == y = Just (2 * x) + +-- x - x == 0 +sAdd x y + | x == -y || -x == y = Just 0 + +-- a + ab = a (b+1) +-- b + ab = b (a+1) +sAdd x (TApp TMul [a,b]) + | x == a = Just (a * (b + 1)) + | x == b = Just (b * (a + 1)) + +-- x + (q + x) = q + 2*x +sAdd x (TApp TAdd [q, y]) + | x == y = Just (q + 2*x) + | x == q = Just (y + 2*x) + +-- -x + (q - x) = q - 2*x +sAdd (TApp TNeg [x]) (TApp TSub [q, y]) + | x == y = Just (q - 2*x) + +-- ay + ab = a (y+b) +-- xb + ab = (x+a) b +sAdd (TApp TMul [x,y]) (TApp TMul [a,b]) + | x == a = Just (x * (y `add` b)) + | y == b = Just ((x `add` a) * y) + where add p q = fromMaybe (p + q) (sAdd p q) + +-- Try associative rule to see if it simplifies things +-- We arbitrarily prefer the first rule below if both apply +-- x + (a+b) = (x+a) + b +-- x + (a+b) = a + (x+b) +sAdd x (TApp TAdd [a,b]) = case fmap (+ b) (sAdd x a) `mplus` fmap (a +) (sAdd x b) of + Nothing -> Nothing + r@(Just (TApp TAdd [t1, t2])) -> maybe r Just (sAdd t1 t2) + r -> r + +-- x + (a - b) = (x+a) - b +sAdd x (TApp TSub [a,b]) = fmap (subtract b) (sAdd x a) + +-- Otherwise we weren't able to do anything smart; so just report Nothing +sAdd _ _ = Nothing + + +-- Walk a deep-tree and add a term; taking care of collapsing if possible. +walkAdd :: (Eq a, Show a) => Term a -> Term a -> Term a +walkAdd a b = case walk b of + Just b' -> b' + Nothing -> bin TAdd a b + where walk t | a == t = Just (2*t) + walk (TApp TMul [c, t]) | a == t = Just $ (c+1) * t + walk t = case t of + TApp TAdd [x, y] -> case walk x of + Just x' -> Just $ bin TAdd x' y + Nothing -> case walk y of + Just y' -> Just $ bin TAdd x y' + Nothing -> Nothing + _ -> Nothing + +-- Split a term into its summands. +-- We use distributivity laws to try to get as smaller terms as possible, +-- in the hope that they might contain fewer varaibles. +-- (i.e., we convert the term to a sum-of-products). +summands :: (Eq a, Show a) => Term a -> [Term a] +summands te = loop [te] + where loop [] = [] + loop (t : ts) = + case t of + TApp op [x,y] + -- x+y; split them + | op == TAdd -> loop (x : y : ts) + -- x-y; split them + | op == TSub -> loop (x : negate y : ts) + -- x*(y1+y2+..+yN); distribute + | op == TMul, composite ys -> loop (map (x *) ys ++ ts) + -- (x1+x2+..+xN)*y; distribute + | op == TMul, composite xs -> loop (map (* y) xs ++ ts) + -- (x1+x2+..+xN)/y; distribute + | op == TDiv, composite xs -> loop (map (/ y) xs ++ ts) + where xs = summands x + ys = summands y + -- Does it have at least two terms? + composite (_ : _ : _) = True + composite _ = False + -- Otherwise keep t; and factor the rest + _ -> t : loop ts + + +-- Split a term into a product of two terms, with the property that +-- only the first term contains the given variable. +factorVar :: ArrVars -> NodeIdx -> Term NodeIdx -> (Term NodeIdx, Term NodeIdx) +factorVar arr x t = + case t of + TApp op ts + | op == TNeg, [t1] <- ts -> let (a,b) = factorVar arr x t1 in (a, negate b) + | op == TMul -> let ([a,b],bs) = unzip $ map (factorVar arr x) ts + in (opt_mul a b, product bs) + | op == TDiv -> let ([a,b],[c,d]) = unzip $ map (factorVar arr x) ts + in (a / b, c / d) + _ | x `IS.member` leavesOfTerm arr t -> (t, 1) + _ -> (1,t) + + where + + opt_mul a b = fromMaybe (a * b) (mul a b) + + mul a b | a == b = Just (a ** 2) + mul (TApp TPow [a,b]) (TApp TPow [a1, c]) | a == a1 = Just (a ** (b + c)) + mul a (TApp TPow [b, c]) | a == b = Just (a ** (1 + c)) + mul (TApp TPow [b,c]) a | a == b = Just (a ** (1 + c)) + mul a (TApp TMul [b,c]) = + case mul a b of + Just b1 -> Just (opt_mul b1 c) + Nothing -> case mul a c of + Just c1 -> Just (opt_mul b c1) + Nothing -> Nothing + mul a (TApp TDiv [b,c]) = Just (opt_mul a b / c) + mul _ _ = Nothing diff --git a/Language/Passage/UI.hs b/Language/Passage/UI.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/UI.hs @@ -0,0 +1,58 @@ +module Language.Passage.UI where + +import qualified Language.Passage.Distribution as D +import Language.Passage.AST +import Control.Monad + +-- | A normal distribution, with a mean and precision +normal :: Expr -> Expr -> BayesianNetwork Expr +normal m t = using (D.normal m t) + +--- | A standard uniform distribution with parameters 0 and 1 +standardUniform :: BayesianNetwork Expr +standardUniform = using D.standardUniform + +--- | A uniform distribution with lower and upper bounds +uniform :: Expr -> Expr -> BayesianNetwork Expr +uniform lo hi = using (D.uniform lo hi) + +-- | A Bernoulli distribution with a mean +bernoulli :: Expr -> BayesianNetwork Expr +bernoulli t = using (D.bernoulli t) + +categorical n ps = using (D.categorical n ps) + +-- | A beta distribution with the given prior sample sizes. +beta :: Expr -> Expr -> BayesianNetwork Expr +beta a b = using (D.beta a b) + +-- | A gamma distribution with the given prior sample sizes. +dgamma :: Expr -> Expr -> BayesianNetwork Expr +dgamma a b = using (D.dgamma a b) + +-- | A chi-square distribution with the given degrees of freedom. +chiSquare :: Expr -> BayesianNetwork Expr +chiSquare df = dgamma (0.5*df) 0.5 + +-- | An exponential distribution with the given rate (inverse scale) +dexp :: Expr -> BayesianNetwork Expr +dexp lambda = dgamma 1 lambda + +-- | A Student's T distribution, given the degrees of freedom. +studentT :: Expr -> BayesianNetwork Expr +studentT df = do + v <- chiSquare df + normal 0 v + +symDirichlet :: Int -> Expr -> BayesianNetwork [Expr] +symDirichlet n alpha = do + gs <- replicateM n $ dgamma alpha 1 + let s = sum gs + return [g/s | g <- gs] + +-- | An improper uniform distribution; has no impact on likelihood +improperUniform :: BayesianNetwork Expr +improperUniform = using D.improperUniform + +improperScale :: BayesianNetwork Expr +improperScale = using D.improperScale diff --git a/Language/Passage/Utils.hs b/Language/Passage/Utils.hs new file mode 100644 --- /dev/null +++ b/Language/Passage/Utils.hs @@ -0,0 +1,63 @@ +module Language.Passage.Utils + ( module Language.Passage.Utils + , module Text.PrettyPrint + ) where + +import Text.PrettyPrint hiding (float,double,rational) +import Data.Ratio(numerator,denominator) +import Data.Word + +-- | A node index +type NodeIdx = Int + +tbd :: String -> a +tbd msg = error $ "TBD: " ++ msg + + + +data Fixity = Prefix | Infix Posn +data Posn = ToLeft | ToRight | None + deriving Eq + + +-- | Pretty printing +class Show t => PP t where + ppPrec :: Rational -> t -> Doc + pp :: t -> Doc + + pp = ppPrec 0 + ppPrec _ = pp + +instance PP Double where + pp = double + +instance PP Int where + pp = int + +instance PP Word where + pp = text . show + +instance PP Integer where + pp = integer + +instance PP Char where + pp = char + +commaSep :: [Doc] -> Doc +commaSep = parens . hsep . punctuate comma + + +ppFrac :: Show a => a -> Doc +ppFrac x = case break (== '.') candidate of + (as,_:bs) | all (== '0') bs -> text as + _ -> text candidate + where candidate = show x + +float :: Float -> Doc +float = ppFrac + +double :: Double -> Doc +double = ppFrac + +rational :: Rational -> Doc +rational r = integer (numerator r) <> text "/" <> integer (denominator r) diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,1 @@ +Passage (Parallel Sampler Generator) is an Bayesian modeling EDSL. Given a model specification and data, Passage generates low-level code for sampling the posterior distribution in parallel using OpenMP threads. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,6 @@ +module Main(main) where + +import Distribution.Simple + +main :: IO () +main = defaultMain diff --git a/cbits/runtime/Makefile b/cbits/runtime/Makefile new file mode 100644 --- /dev/null +++ b/cbits/runtime/Makefile @@ -0,0 +1,21 @@ +.PHONY: all clean sampler + + +all: sample.pdf + evince $< || open $< + +sampler: + make -C src + cp src/sampler $@ + +sample.pdf: histogram.R datafile + R --no-save -f $< + +datafile: sampler + ./sampler > $@ + +clean: + -rm datafile datafile.png + + + diff --git a/cbits/runtime/src/Makefile b/cbits/runtime/src/Makefile new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/Makefile @@ -0,0 +1,21 @@ +-include extra_settings + +CC=gcc + +CPPFLAGS +=-DNDEBUG +CFLAGS +=-Wall -O3 -fomit-frame-pointer -msse2 -fopenmp + +TGT=sampler + +HEADER_FILES=passage.h mt19937ar.h +SRC_FILES=$(wildcard *.c) +OBJ_FILES=$(patsubst %.c,%.o,$(SRC_FILES)) + +$(TGT): $(HEADER_FILES) $(OBJ_FILES) + $(CC) $(OBJ_FILES) -o $(TGT) -lm -fopenmp + +clean: + -rm $(TGT) *.o + + + diff --git a/cbits/runtime/src/generic_slicers.c b/cbits/runtime/src/generic_slicers.c new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/generic_slicers.c @@ -0,0 +1,205 @@ +#include "passage.h" + +/* This duplicates code from the templates. + It'd be nice to eliminate this, at a later stage. +*/ + +double slice_real(double (*ll)(double), double width, double x0) +{ + double x1; + double r; + double y; + + double lo; + + double left; + double right; + + r = genrand_real3(); + y = log(r) + ll(x0); + + lo= getRandomRange(0, width); + + left = x0 - lo; + while(ll(left) > y) left -= width; + + right = x0 - lo + width; + while(ll(right) > y) right += width; + + + for ( x1 = getRandomRange(left, right) + ; ll(x1) < y + ; x1 = getRandomRange(left, right) + ) { + if (x1 < x0) left = x1; else right = x1; + } + + return x1; +} + + +double tune_slice_real(double (*ll)(double), double *width, double x0) +{ + double x1; + double r; + double y; + + double lo; + + double left; + double right; + + int steps_out = 0; + int steps_in = 0; + + r = genrand_real3(); + y = log(r) + ll(x0); + + lo= getRandomRange(0, *width); + + left = x0 - lo; + while(ll(left) > y) { left -= *width; ++steps_out; } + + right = x0 - lo + *width; + while(ll(right) > y) { right += *width; ++steps_out; } + + + for ( x1 = getRandomRange(left, right) + ; ll(x1) < y + ; x1 = getRandomRange(left, right) + ) { + if (x1 < x0) left = x1; else right = x1; + ++steps_in; + } + + *width *= 0.75 + steps_out / ((double) (1 << (steps_in+2))); + + + return x1; +} + +double slice_pos_real + (double (*ll)(double), double width, double left, double x0) +{ + double x1; + double r; + double y; + double lo; + double right; + + r = genrand_real3(); + y = log(r) + ll(x0); + lo= getRandomRange(0, width); + right = x0 - lo + width; + while(ll(right) > y) right += width; + + for ( x1 = getRandomRange(left, right) + ; ll(x1) < y + ; x1 = getRandomRange(left, right) + ) { + if (x1 < x0) left = x1; else right = x1; + } + + return x1; +} + +double tune_slice_pos_real + (double (*ll)(double), double *w, double left, double x0) +{ + double x1; + double r; + double y; + double lo; + double right; + + int steps_out = 0; + int steps_in = 0; + + r = genrand_real3(); + y = log(r) + ll(x0); + lo= getRandomRange(0, *w); + + right = x0 - lo + *w; + while(ll(right) > y) { right += *w; ++steps_out; } + + for ( x1 = getRandomRange(left, right) + ; ll(x1) < y + ; x1 = getRandomRange(left, right) + ) { + if (x1 < x0) left = x1; else right = x1; + ++steps_in; + } + + /* Trying to optimize the width, + assuming an average "step in" cuts the interval in half. + Current rule of thumb is to take the weighted average of the current width + with the estimated ideal width */ + *w *= 0.75 + steps_out / ((double) (1 << (steps_in+2))); + + /* fprintf(stderr, "(%d, [(%f,%d)])\n", VAR, WIDTH(VAR), error); */ + return x1; +} + + + +double slice_real_left_right + ( double (*ll)(double) + , double left + , double right + , double x0 + ) +{ + double x1; + double r; + double y; + + r = genrand_real3(); + y = log(r) + ll(x0); + + for ( x1 = getRandomRange(left, right) + ; ll(x1) < y + ; x1 = getRandomRange(left, right) + ) + { + if (x1 < x0) left = x1; else right = x1; + } + + return x1; +} + +double slice_discrete_right(double (*ll)(double), double n, double x) { + double y; + y = floor(getRandomRangeOpenRight(0,n+1)); + return (log(genrand_real3()) < ll(y) - ll(x)) ? y : x; +} + +double slice_discrete(double (*ll)(double), double x0) +{ + double y0; + y0 = ll(x0); + + // x1 is the proposal for the next value + double x1 = -x0 * log(genrand_real3()); + double y1 = ll(x1); + + double rat = x0 / x1; + // logP is the log of the probability of transitioning to x1 + double logP = y1 - y0 + log(rat) + 1/rat - rat; + + if(logP > 0) { + return x1; + } + else { + if(log(genrand_real3()) < logP) { + return x1; + } + else { + return x0; + } + } +} + + + + + diff --git a/cbits/runtime/src/kiss.h b/cbits/runtime/src/kiss.h new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/kiss.h @@ -0,0 +1,34 @@ +/*****************************************************************************/ +/* Implementation of a 32-bit KISS generator which uses no multiply instructions */ + +extern unsigned int rx, ry, rz, rw, rc; +#pragma omp threadprivate (rx,ry,rz,rw,rc) + +static inline +void init_genrand(unsigned long seed) { + rx = seed; + ry = seed + 1; + rz = seed + 2; + rw = seed + 3; + rc = seed & 1; +} + +static inline +unsigned long genrand_int32(void) { + + int t; + ry ^= ry << 5; + ry ^= ry >> 7; + ry ^= ry << 22; + t = rz + rw + rc; + rz = rw; + rc = t < 0; + rw = t & 2147483647; + rx += 1411392427; + return rx + ry + rw; +} + + + + + diff --git a/cbits/runtime/src/mt19937ar.h b/cbits/runtime/src/mt19937ar.h new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/mt19937ar.h @@ -0,0 +1,125 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + Copyright (C) 2005, Mutsuo Saito, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of its contributors may not 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. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +/* This code has been modified, so that it works in a parallel setting. + The changes are as follows: + - The functions are parmetrized by a state paramemter + - We removed the check to see if things have been initialized +*/ + + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ + + +extern unsigned long mt[N]; /* the array for the state vector */ +extern int mti; /* mti==N+1 means mt[N] is not initialized */ +extern unsigned long mag01[2]; +#pragma omp threadprivate(mt,mti,mag01) + +/* initializes mt[N] with a seed */ +static inline +void init_genrand(unsigned long seed) +{ + unsigned long i; + + // Initialize mag01 here, because global initializers do not work with OMP + /* mag01[x] = x * MATRIX_A for x=0,1 */ + mag01[0] = 0x0UL; + mag01[1] = MATRIX_A; + + mt[0]= seed & 0xffffffffUL; + for (i=1; i> 30)) + i; + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[i] &= 0xffffffffUL; + /* for >32 bit machines */ + } + mti = N; +} + +/* generates a random number on [0,0xffffffff]-interval */ +static inline +unsigned long genrand_int32(void) +{ + unsigned long y; + + if (mti >= N) { /* generate N words at one time */ + int kk; + + for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; + } + for (;kk> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + + + diff --git a/cbits/runtime/src/passage.c b/cbits/runtime/src/passage.c new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/passage.c @@ -0,0 +1,141 @@ +#include "passage.h" +#include +#include +#include +#include +#include +#include +#include + +// An arbitrary limit to the number of seeds we might have +#define MAX_SEEDS 1024 + +unsigned long number_of_samples; +unsigned long steps_per_sample; +unsigned long warm_up_steps; +static unsigned long seeds[MAX_SEEDS]; +unsigned long num_threads; +int have_seed; + + +/* State for the random number generator */ + +#ifdef __USE_MERSENNE +unsigned long mag01[2]; +unsigned long mt[N]; /* the array for the state vector */ +int mti; /* mti==N+1 means mt[N] is not initialized */ +#else +unsigned int rx, ry, rz, rw, rc; +#endif + + + + + + + +/* Make up a random seed */ +static +unsigned long getSeed() { + struct timeval tv; + unsigned long seed; + int problem = 1; + int fd; + + if ((fd = open("/dev/urandom", O_RDONLY)) >= 0 || + (fd = open("/dev/random", O_RDONLY)) >= 0) { + ssize_t n = read(fd, &seed, sizeof(seed)); + if (n == sizeof(seed)) problem = 0; + close(fd); + } + if (problem) { + gettimeofday(&tv, NULL); + seed = (getpid() << 16) ^ tv.tv_sec ^ tv.tv_usec; + } + + return seed; +} + +static +int get_num(unsigned long *out) { + char *end; + unsigned long r; + r = strtoul(optarg, &end, 10); + if (*end == '\0' && optarg != '\0') { *out = r; return 1; } + return 0; +} + +/* Generated by DSL */ +extern void sampler(void); +extern void set_defaults(void); +extern void init_vars(void); + +int main(int argc, char *argv[]) { + int r, i; + + set_defaults(); + + do { + r = getopt(argc, argv, "n:i:w:s:h"); + if (r == -1) break; + + switch (r) { + case 'n': if (!get_num(&number_of_samples)) goto err; break; + case 'i': if (!get_num(&steps_per_sample)) goto err; break; + case 'w': if (!get_num(&warm_up_steps)) goto err; break; + case 's': + if (have_seed >= MAX_SEEDS) goto err; + if (!get_num(&seeds[have_seed])) goto err; ++have_seed; + break; + + default: goto err; + } + + } while (1); + if (optind != argc) goto err; + + if (num_threads >= MAX_SEEDS) { + fprintf(stderr, "Too many threads (limit %d)\n", MAX_SEEDS); + return 2; + } + + for (; have_seed < num_threads; ++have_seed) + seeds[have_seed] = getSeed(); + + for (i = 0; i < have_seed; ++i) + fprintf(stderr, "Seed[%d]: %lu\n", i, seeds[i]); + fprintf(stderr, "Samples: %lu\n", number_of_samples); + fprintf(stderr, "Steps/sample: %lu\n", steps_per_sample); + fprintf(stderr, "Warm-up: %lu\n", warm_up_steps); + fprintf(stderr, "Threads: %lu\n", num_threads); + + #pragma omp parallel num_threads(num_threads) + { + init_genrand(seeds[omp_get_thread_num()]); + + #pragma omp single + init_vars(); + + sampler(); + } + + return 0; + +err: + fprintf(stderr, "usage:\n"); + fprintf(stderr, " -n NUM\tNumber of samples to generate.\n"); + fprintf(stderr, " -i NUM\tNumber of iterations per sample.\n"); + fprintf(stderr, " -w NUM\tNumber of iterations for calibrating sampling window.\n"); + fprintf(stderr, " -s NUM\tFixed seed(s) of randomness.\n"); + fprintf(stderr, " -h \tThis help.\n"); + return 1; +} + + + +void crash_out_of_bounds(int line) { + fprintf(stderr, "Array index out of bounds on line %d\n", line); + exit(1); +} + + diff --git a/cbits/runtime/src/passage.h b/cbits/runtime/src/passage.h new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/passage.h @@ -0,0 +1,101 @@ +#ifndef __BAYESIAN_DSL_H_INCLUDED +#define __BAYESIAN_DSL_H_INCLUDED + +#include +#include +#include +#include + +extern unsigned long number_of_samples; +extern unsigned long steps_per_sample; +extern unsigned long warm_up_steps; +extern unsigned long num_threads; +extern int have_seed; + +#ifdef __USE_MERSENNE +#include "mt19937ar.h" +#else +#include "kiss.h" +#endif + +/* generates a random number on [0,1]-real-interval */ +static inline +double genrand_real1(void) { + return genrand_int32()*(1.0/4294967295.0); + /* divided by 2^32-1 */ +} + +/* generates a random number on [0,1)-real-interval */ +static inline +double genrand_real2() { + return genrand_int32()*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on (0,1)-real-interval */ +static inline +double genrand_real3() { + return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* Return a double-random value in [lo, hi] */ +static inline +double getRandomRange(double lo, double hi) { + return lo + (hi - lo) * genrand_real1(); +} + +/* Return a double-random value in [lo, hi) */ +static inline +double getRandomRangeOpenRight(double lo, double hi) { + return lo + (hi - lo) * genrand_real2(); +} + +static inline +double square (double x) { return x * x; } + + +static inline +void progress(unsigned long n) { + int l = fprintf(stderr, "%lu", n); + l += fprintf(stderr," of %lu (%2lu%%)", number_of_samples, + 100 * n / number_of_samples); + for (; l > 0; --l) putc('\b', stderr); +} + + +void crash_out_of_bounds(int line) __attribute__((noreturn)); + +/* Generic samplers */ + + +double slice_real(double (*ll)(double), double width, double x0); +double tune_slice_real(double (*ll)(double), double *width, double x0); +double slice_pos_real + (double (*ll)(double), double width, double left, double x0); +double tune_slice_pos_real + (double (*ll)(double), double *w, double left, double x0); +double slice_real_left_right + ( double (*ll)(double) + , double left + , double right + , double x0 + ); + +double slice_discrete_right(double (*ll)(double), double n, double x); +double slice_discrete(double (*ll)(double), double x0); + +#define SLICE_NAME(var) slice_##var +#define SLICE_TUNE_NAME(var) slice_tune_##var +#define INIT_DET_VARS_NAME(var) init_det_vars_##var +#define LL_FUN_NAME(var) ll_##var +#define WIDTH_NAME(var) w_##var + +#define SLICE(x) SLICE_NAME(x) +#define SLICE_TUNE(x) SLICE_TUNE_NAME(x) +#define INIT_DET_VARS(x) INIT_DET_VARS_NAME(x) +#define LL_FUN(x) LL_FUN_NAME(x) +#define WIDTH(x) WIDTH_NAME(x) + + +#endif /* __BAYESIAN_DSL_H_INCLUDED */ diff --git a/cbits/runtime/src/templates/finiteMetropolis.c b/cbits/runtime/src/templates/finiteMetropolis.c new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/templates/finiteMetropolis.c @@ -0,0 +1,10 @@ +/* Independent Metropolis sampling over [0..n-1] */ +/* static inline */ +double SLICE(VAR)(double x) { + double y; + double n = RIGHT; + y = floor(getRandomRangeOpenRight(0,n+1)); + return (log(genrand_real3()) < LL_FUN(VAR)(y) - LL_FUN(VAR)(x)) ? y : x; +} + + diff --git a/cbits/runtime/src/templates/metropolis_posreal.c b/cbits/runtime/src/templates/metropolis_posreal.c new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/templates/metropolis_posreal.c @@ -0,0 +1,28 @@ +/* Metropolis-Hastings over the positive reals */ +/* static inline */ +double SLICE(VAR)(double x0) +{ + double y0; + + y0 = LL_FUN(VAR)(x0); + + // x1 is the proposal for the next value + double x1 = -x0 * log(genrand_real3()); + double y1 = LL_FUN(VAR)(x1); + + double rat = x0 / x1; + // logP is the log of the probability of transitioning to x1 + double logP = y1 - y0 + log(rat) + 1/rat - rat; + + if(logP > 0) { + return x1; + } + else { + if(log(genrand_real3()) < logP) { + return x1; + } + else { + return x0; + } + } + diff --git a/cbits/runtime/src/templates/slice.c b/cbits/runtime/src/templates/slice.c new file mode 100644 --- /dev/null +++ b/cbits/runtime/src/templates/slice.c @@ -0,0 +1,101 @@ +#if !(defined(LEFT) && defined(RIGHT)) +static double WIDTH(VAR) = 1; +#endif + +/* Slice over the reals, or a sub-range of the reals. */ +/* static inline */ +double SLICE(VAR)(double x0) +{ + double x1; + double r; + double y; +#if !(defined(LEFT) && defined(RIGHT)) + double lo; +#endif + double left; + double right; + + r = genrand_real3(); + y = log(r) + LL_FUN(VAR)(x0); +#if !(defined(LEFT) && defined(RIGHT)) + lo= getRandomRange(0, WIDTH(VAR)); +#endif + +#if defined(LEFT) + left = LEFT; +#else + left = x0 - lo; + while(LL_FUN(VAR)(left) > y) left -= WIDTH(VAR); +#endif + +#if defined(RIGHT) + right = RIGHT; +#else + right = x0 - lo + WIDTH(VAR); + while(LL_FUN(VAR)(right) > y) right += WIDTH(VAR); +#endif + + for ( x1 = getRandomRange(left, right) + ; LL_FUN(VAR)(x1) < y + ; x1 = getRandomRange(left, right) + ) { + if (x1 < x0) left = x1; else right = x1; + } + + return x1; +} + + +#if !(defined(LEFT) && defined(RIGHT)) + +/* Slice over the reals */ +/* static inline */ +double SLICE_TUNE(VAR)(double x0) +{ + double x1; + double r; + double y; + double lo; + double left; + double right; + + int steps_out = 0; + int steps_in = 0; + + r = genrand_real3(); + y = log(r) + LL_FUN(VAR)(x0); + lo= getRandomRange(0, WIDTH(VAR)); + +#if defined(LEFT) + left = LEFT; +#else + left = x0 - lo; + while(LL_FUN(VAR)(left) > y) { left -= WIDTH(VAR); ++steps_out; } +#endif + +#if defined(RIGHT) + right = RIGHT; +#else + right = x0 - lo + WIDTH(VAR); + while(LL_FUN(VAR)(right) > y) { right += WIDTH(VAR); ++steps_out; } +#endif + + for ( x1 = getRandomRange(left, right) + ; LL_FUN(VAR)(x1) < y + ; x1 = getRandomRange(left, right) + ) { + if (x1 < x0) left = x1; else right = x1; + ++steps_in; + } + + /* Trying to optimize the width, + assuming an average "step in" cuts the interval in half. + Current rule of thumb is to take the weighted average of the current width + with the estimated ideal width */ + WIDTH(VAR) *= 0.75 + steps_out / ((double) (1 << (steps_in+2))); + + /* fprintf(stderr, "(%d, [(%f,%d)])\n", VAR, WIDTH(VAR), error); */ + return x1; +} + +#endif diff --git a/passage.cabal b/passage.cabal new file mode 100644 --- /dev/null +++ b/passage.cabal @@ -0,0 +1,54 @@ +Name: passage +Version: 0.1 +Category: Statistical Modeling, Code Generation +Synopsis: Parallel code generation for hierarchical Bayesian modeling. +Description: + Passage is a PArallel SAmpler GEnerator. The user specifies a hierarchical + Bayesian model and data using the Passage EDSL, and Passage generates code + to sample the posterior distribution in parallel. + . + Currently Passage targets C with OpenMP threads. +Copyright: 2011 Galois, Inc. and Battelle Memorial Institute +License: BSD3 +License-file: LICENSE +Stability: Experimental +Author: Chad Scherrer (Pacific Northwest National Laboratory), + Levent Erkok (Galois, Inc), + Iavor Diatchki (Galois, Inc), + Matthew Sottile (Galois, Inc) +Maintainer: Chad Scherrer +Build-Type: Simple +Cabal-Version: >= 1.6 +Extra-Source-Files: README, COPYRIGHT, LICENSE +Data-Files: cbits/runtime/src/mt19937ar.h + , cbits/runtime/src/kiss.h + , cbits/runtime/src/passage.h + , cbits/runtime/src/passage.c + , cbits/runtime/src/generic_slicers.c + , cbits/runtime/src/templates/slice.c + , cbits/runtime/src/templates/finiteMetropolis.c + , cbits/runtime/src/templates/metropolis_posreal.c + , cbits/runtime/src/Makefile + , cbits/runtime/Makefile + +Library + ghc-options : -Wall + ghc-prof-options: -auto-all -caf-all + Build-Depends : base >= 3 && < 5, containers, monadLib, pretty, + process, primitive, filepath, directory, random, GraphSCC, + mwc-random, array + Exposed-modules : Language.Passage + , Language.Passage.Distribution + , Language.Passage.UI + Other-modules : Language.Passage.AST + , Language.Passage.Graph + , Language.Passage.Graph2C + , Language.Passage.Lang.LaTeX + , Language.Passage.Lang.C + , Language.Passage.SliceSample + , Language.Passage.SliceSampleMWC + , Language.Passage.Term + , Language.Passage.Utils + , Language.Passage.SimulatorConf + , Language.Passage.GraphColor + , Paths_passage