packages feed

passage (empty) → 0.1

raw patch · 29 files changed

+3657/−0 lines, 29 filesdep +GraphSCCdep +arraydep +basesetup-changed

Dependencies added: GraphSCC, array, base, containers, directory, filepath, monadLib, mwc-random, pretty, primitive, process, random

Files

@@ -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.
+ LICENSE view
@@ -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+
+ Language/Passage.hs view
@@ -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."+++
+ Language/Passage/AST.hs view
@@ -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.++-}++
+ Language/Passage/Distribution.hs view
@@ -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+    }
+ Language/Passage/Graph.hs view
@@ -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))++
+ Language/Passage/Graph2C.hs view
@@ -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 <math.h>"+     cpp "#include <stdio.h>"+     cpp "#include <omp.h>"+     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)+        ]+
+ Language/Passage/GraphColor.hs view
@@ -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+++
+ Language/Passage/Lang/C.hs view
@@ -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+
+ Language/Passage/Lang/LaTeX.hs view
@@ -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]
+ Language/Passage/SimulatorConf.hs view
@@ -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 })+++
+ Language/Passage/SliceSample.hs view
@@ -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
+ Language/Passage/SliceSampleMWC.hs view
@@ -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)
+ Language/Passage/Term.hs view
@@ -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
+ Language/Passage/UI.hs view
@@ -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
+ Language/Passage/Utils.hs view
@@ -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)
+ README view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,6 @@+module Main(main) where++import Distribution.Simple++main :: IO ()+main = defaultMain
+ cbits/runtime/Makefile view
@@ -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+++
+ cbits/runtime/src/Makefile view
@@ -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+++
+ cbits/runtime/src/generic_slicers.c view
@@ -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;+    }+  }+}+++++
+ cbits/runtime/src/kiss.h view
@@ -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;+}+++++
+ cbits/runtime/src/mt19937ar.h view
@@ -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<N; i++) {+        mt[i] = 1812433253UL * (mt[i-1] ^ (mt[i-1] >> 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<N-M;kk++) {+            y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);+            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];+        }+        for (;kk<N-1;kk++) {+            y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);+            mt[kk] = mt[kk+(M-N)] ^ (y >> 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;+}+++
+ cbits/runtime/src/passage.c view
@@ -0,0 +1,141 @@+#include "passage.h"+#include <math.h>+#include <time.h>+#include <stdlib.h>+#include <fcntl.h>+#include <time.h>+#include <unistd.h>+#include <sys/time.h>++// 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);+}++
+ cbits/runtime/src/passage.h view
@@ -0,0 +1,101 @@+#ifndef __BAYESIAN_DSL_H_INCLUDED+#define __BAYESIAN_DSL_H_INCLUDED++#include <stdlib.h>+#include <stdio.h>+#include <math.h>+#include <omp.h>++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 */
+ cbits/runtime/src/templates/finiteMetropolis.c view
@@ -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;+}++
+ cbits/runtime/src/templates/metropolis_posreal.c view
@@ -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;+    }+  }+
+ cbits/runtime/src/templates/slice.c view
@@ -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
+ passage.cabal view
@@ -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