hbayes 0.1 → 0.2

raw patch · 13 files changed

+1200/−683 lines, 13 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Bayes.hs +50/−4
• Bayes/Examples.hs +6/−1
• Bayes/Examples/Tutorial.hs +19/−10
• Bayes/Factor.hs +130/−257
• Bayes/FactorElimination.hs +166/−402
• Bayes/FactorElimination/JTree.hs +486/−0
• Bayes/ImportExport/HuginNet.hs +2/−1
• Bayes/PrivateTypes.hs +184/−0
• Bayes/Test.hs +12/−0
• Bayes/Test/CompareEliminations.hs +2/−2
• Bayes/Test/ReferencePatterns.hs +122/−0
• Bayes/VariableElimination.hs +2/−2
• hbayes.cabal +19/−4

@@ -35,6 +35,7 @@   , newEdge   , edgeEndPoints   , connectedGraph+  , dag   -- * SimpleGraph implementation   -- ** The SimpleGraph type   , DirectedSG@@ -63,7 +64,7 @@ import Control.Monad.State.Strict import Control.Monad.Writer.Strict import Control.Applicative((<$>))-import Bayes.Factor+import Bayes.Factor hiding(isEmpty) import Data.Maybe import qualified Data.Map as Map import qualified Data.Foldable as F@@ -74,6 +75,7 @@ import Test.QuickCheck import Test.QuickCheck.Arbitrary import Data.List(sort,intercalate,nub)+import Bayes.PrivateTypes hiding(isEmpty)  --import Debug.Trace --debug a = trace (show a) a@@ -120,6 +122,29 @@      foldM createEdge g edges    +-- | Warning : the generated graph is not at all a bayesian network+-- The variables in the CPT have no reason to correspond to the edges+-- connected to that CPT.+-- Only the main variable (first variable) is linked to the right vertex+instance Factor f => Arbitrary (DirectedSG () f) where+  arbitrary = do +    let createVertex g i = do +          let value = fromJust $ factorWithVariables [DV (Vertex i) 2] [0.1,0.9]+          return $ addVertex (Vertex i) value g+        createEdge g (va,vb) = do +          return $ addEdge (edge va vb) () g ++    nbVertex <- choose (1,8) :: Gen Int+    +    g <- foldM createVertex emptyGraph [1..nbVertex]++    let allPairs = [(Vertex x,Vertex y) | x <- [1..nbVertex], y <- [1..nbVertex], x /= y]+        anEdge (x,y) = arbitrary :: Gen Bool++    edges <- filterM anEdge allPairs++    foldM createEdge g edges+ testEdgeRemoval_prop :: DirectedSG String String -> Property testEdgeRemoval_prop g = (not . hasNoEdges) g ==>    let Just e = someEdge g@@ -248,7 +273,26 @@     ingoing :: g a b -> Vertex -> Maybe [Edge]     outgoing :: g a b -> Vertex -> Maybe [Edge] +-- | Get the root node for the graph+rootNode :: DirectedGraph g => g a b -> Maybe Vertex+rootNode g = +  let someRoots = filter (isRoot g) . allVertices $ g+  in +  case someRoots of +    (h:l) -> Just h +    _ -> Nothing+  where +    isRoot g v =+      case ingoing g v of +        Just [] -> True +        _ -> False +-- | Check if the graph is a directed Acyclic graph+dag :: DirectedGraph g => g a b -> Bool +dag g = case rootNode g of +  Nothing -> isEmpty g +  Just r -> dag (removeVertex r g)+ -- | Check if the graph is connected connectedGraph :: Graph g => g a b -> Bool  connectedGraph g = @@ -384,12 +428,14 @@ class FoldableWithVertex g where   -- | Fold with vertex    foldrWithVertex :: (Vertex -> a -> b -> b) -> b -> g c a -> b +  foldlWithVertex' :: (b -> Vertex -> a -> b) -> b -> g c a -> b   instance FoldableWithVertex (SimpleGraph local) where   foldrWithVertex f s (SP _ vm _) = IM.foldrWithKey (\k (_,v) y -> f (Vertex k) v y) s vm+  foldlWithVertex' f s (SP _ vm _) = IM.foldlWithKey' (\y k (_,v)  -> f y (Vertex k) v) s vm  _addLabeledVertex vertexName vert@(Vertex v) value (SP em vm name) =-  let vm' = IM.insertWith noRedundancy v (emptyNeighborhood,value) vm+  let vm' = IM.insertWith' noRedundancy v (emptyNeighborhood,value) vm       name' = IM.insert v vertexName name    in   SP em vm' name'@@ -498,7 +544,7 @@     else        Just . head . M.keys $ em -_addVertex vert@(Vertex v) value (SP em vm nm) = SP em (IM.insertWith noRedundancy v (emptyNeighborhood,value) vm) nm+_addVertex vert@(Vertex v) value (SP em vm nm) = SP em (IM.insertWith' noRedundancy v (emptyNeighborhood,value) vm) nm  _removeVertex v@(Vertex vertex) g@(SP _ vm _)  = maybe g removeVertexWithValue (IM.lookup vertex vm)   where@@ -587,7 +633,7 @@   tell "\n"   let r = IM.lookup k nm   when (isJust r) $ do-     tell $ fromJust r+     tell $ "Node " ++ fromJust r   tell "\n"   tell $ show v   tell "\n"

@@ -69,7 +69,9 @@ module Bayes.Examples(    example  , exampleJunction+#ifndef LOCAL  , exampleImport+#endif  , exampleDiabete  , exampleAsia  , examplePoker@@ -86,6 +88,8 @@ import qualified Data.Map as Map import System.Directory(getHomeDirectory) import System.FilePath((</>))++#ifndef LOCAL import Paths_hbayes  -- | Example showing how to import a graph described into@@ -95,6 +99,7 @@     path <- getDataFileName "cancer.net"     r <- importBayesianGraph path     return (runBN $ fromJust r)+#endif  -- | Genereic loading functions to load some other -- examples from the author's dropbox.@@ -136,7 +141,7 @@   -- | Standard example found in many books about Bayesian Networks.-example :: (DVSet,SBN CPT)+example :: ([DV],SBN CPT) example = runBN $ do      winter <- variable "winter" (t :: Bool)     sprinkler <- variable "sprinkler" (t :: Bool) 

@@ -124,16 +124,23 @@ module Bayes.Examples.Tutorial(     -- * Tests with the standard network        inferencesOnStandardNetwork+#ifndef LOCAL     -- * Tests with the cancer network     , inferencesOnCancerNetwork+#endif     , Coma(..)     , miscTest+--    , miscDiabete 	) where   import Bayes.Factor import Bayes import Bayes.VariableElimination+#ifndef LOCAL import Bayes.Examples(example, exampleJunction,exampleImport,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)+#else +import Bayes.Examples(example, exampleJunction,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)+#endif import Bayes.FactorElimination import Data.Function(on) import qualified Data.Map as Map@@ -141,11 +148,11 @@ import System.Exit(exitSuccess) import qualified Data.List as L((\\)) -miscDiabete = do -  (varmap,perso) <- exampleDiabete-  let jtperso = createJunctionTree nodeComparisonForTriangulation perso-      cho0 = fromJust . Map.lookup "cho_0" $ varmap-  print $ posterior jtperso cho0+--miscDiabete = do +--  (varmap,perso) <- exampleDiabete+--  let jtperso = createJunctionTree nodeComparisonForTriangulation perso+--      cho0 = fromJust . Map.lookup "cho_0" $ varmap+--  print $ posterior jtperso cho0  miscTest s = do    (varmap,perso) <- anyExample s@@ -171,6 +178,7 @@ -- from the cancer network. data Coma = Present | Absent deriving(Eq,Enum,Bounded) +#ifndef LOCAL -- | Inferences with the cancer network inferencesOnCancerNetwork = do      print "CANCER NETWORK"@@ -182,12 +190,13 @@     mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer $ x)) [varA,varB,varC,varD,varE]      print "UPDATED EVIDENCE : Coma present"-    let jtcancer' = updateEvidence [varD =: Present] jtcancer +    let jtcancer' = changeEvidence [varD =: Present] jtcancer      mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varD,varE]      print "UPDATED EVIDENCE : Coma absent"-    let jtcancer' = updateEvidence [varD =: Absent] jtcancer +    let jtcancer' = changeEvidence [varD =: Absent] jtcancer      mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varD,varE]+#endif   -- | Inferences with the standard network inferencesOnStandardNetwork = do@@ -221,20 +230,20 @@     print m     putStrLn "" -    let jt' = updateEvidence [wet =: True] jt +    let jt' = changeEvidence [wet =: True] jt       print "Posterior Marginal : probability of rain if grass wet"     let m = posterior jt' rain     print m     putStrLn "" -    let jt'' = clearEvidence jt'+    let jt'' = changeEvidence [] jt'     print "Prior Marginal : probability of rain"     let m = posterior jt rain     print m     putStrLn "" -    let jt3 = updateEvidence [wet =: True, sprinkler =: True] jt'+    let jt3 = changeEvidence [wet =: True, sprinkler =: True] jt'      print "Posterior Marginal : probability of rain if grass wet and sprinkler used"     let m = posterior jt3 rain

@@ -16,10 +16,9 @@  -- * Implementation  , Vertex(..)  -- ** Discrete variables and instantiations- , DV(..)- , DVSet(..)+ , DV+ --, DVSet(..)  , DVI- , DVISet(..)  , setDVValue  , instantiationValue  , instantiationVariable@@ -40,7 +39,7 @@  import qualified Data.Vector.Unboxed as V import Data.Vector.Unboxed((!))-import Data.Maybe(fromJust,mapMaybe)+import Data.Maybe(fromJust,mapMaybe,isJust) import qualified Data.List as L import Text.PrettyPrint.Boxes hiding((//)) import Test.QuickCheck@@ -48,79 +47,19 @@ import qualified Data.IntMap as IM import Control.Monad import System.Random(Random)+import Data.List(partition)+import Bayes.PrivateTypes  --import Debug.Trace  --debug a = trace ("\nDEBUG\n" ++ show a ++ "\n") a --- | Vertex type used to identify a vertex in a graph-newtype Vertex = Vertex {vertexId :: Int} deriving(Eq,Ord) -instance Show Vertex where -    show (Vertex v) = "v" ++ show v---- | A Set of variables used in a factor. s is the set and a the variable-class Set s where-    -- | Empty set-    emptySet :: s a-    -- | Union of two sets-    union :: Eq a => s a -> s a -> s a-    -- | Intersection of two sets-    intersection :: Eq a => s a -> s a -> s a-    -- | Difference of two sets-    difference :: Eq a => s a -> s a -> s a-    -- | Check if the set is empty-    isEmpty :: s a -> Bool-    -- | Check if an element is member of the set-    isElem :: Eq a => a -> s a -> Bool-    -- | Add an element to the set-    addElem :: Eq a => a -> s a -> s a-    -- | Number of elements in the set-    nbElements :: s a -> Int--    -- | Check if a set is subset of another one-    subset :: Eq a => s a -> s a -> Bool--    -- | Check set equality-    equal :: Eq a => s a -> s a -> Bool-    equal sa sb = (sa `subset` sb) && (sb `subset` sa)--instance Set [] where-    emptySet = []-    union = L.union-    intersection = L.intersect-    difference a b = a L.\\ b-    isEmpty [] = True -    isEmpty _ = False-    isElem = L.elem -    addElem a l = if a `elem` l then l else a:l-    nbElements = length-    subset sa sb = all (`elem` sb) sa---- | A discrete variable has a number of levels which is required to size the factors-class BayesianDiscreteVariable v where-    dimension :: v -> Int -- -- | A vertex associated to another value (variable dimension, variable value ...) class LabeledVertex l where     variableVertex :: l -> Vertex --- | A discrete variable-data DV = DV !Vertex !Int deriving(Eq,Ord) --- | A set of discrete variables-type DVSet = [DV]--instance Show DV where-    show (DV v d) = show v ++ "(" ++ show d ++ ")"---- | Discrete Variable instantiation. A variable and its value-data DVI a = DVI DV !a deriving(Eq)--instance Show a => Show (DVI a) where -   show (DVI (DV v _) i) = show v ++ "=" ++ show i- -- | Convert a variable instantation to a factor -- Useful to create evidence factors factorFromInstantiation :: Factor f => DVI Int -> f@@ -129,88 +68,16 @@     in     fromJust . factorWithVariables [dv] . map (setValue) $ [0..dimension dv-1] --- | A set of variable instantiations-type DVISet a = [DVI a] -instance BayesianDiscreteVariable DV where-    dimension (DV _ d) = d --- | Create a discrete variable instantiation for a given discrete variable-setDVValue :: DV -> a -> DVI a-setDVValue v a = DVI v a -getMinBound :: Bounded a => a -> a -getMinBound _ = minBound---- | Create a variable instantiation using values from--- an enumeration-(=:) :: (Bounded b, Enum b) => DV -> b -> DVI Int -(=:) a b = setDVValue a (fromEnum b - fromEnum (getMinBound b))---- | Extract value of the instantiation-instantiationValue (DVI _ v) = v---- | Discrete variable from the instantiation-instantiationVariable (DVI dv _) = dv- instance LabeledVertex (DVI a) where     variableVertex (DVI v _) = variableVertex v  instance LabeledVertex DV where     variableVertex (DV v _) = v --- | Extend indexing to full variable set using a bool--- list and a default value--- For instance [True, False, True, False] 5 [2,3] ---> [2,5,3,5]-extend :: [Bool] -> a -> [a] -> [a]-extend [] _ l = l-extend (h:t) d [] = d:extend t d []-extend (False:t) d l = d:extend t d l-extend (True:t) d (h:l') = h:extend t d l' --- | Inner loop function using full indices for full variables-type InnerLoop a = [Int] -> a---- | Outer loop function using result from inner loop--- and outer vars indices-type OuterLoop a b = [Int] -> [a] -> b---- | Iter on outer var and inner var--- Inner body is called with indiced for full set--- Outer body is called with indices for outer set-forSubA :: DVSet -- ^ All variables-        -> DVSet -- ^ Outer variables-        -> (DVSet -> [Int] -> [a]) -- ^ Inner loop body-        -> OuterLoop a b -- ^ Outer loop function-        -> [b]-forSubA allvars outervars inner outer = -    let sCode s e = if (e `isElem` s) then True else False-        selection s = map (sCode s) allvars-        computeOuter iouter =-            let outerIdx =  extend (selection outervars) 0 iouter-                innerValues = inner allvars outerIdx-            in -            outer iouter innerValues-    in-    map computeOuter (forAllIndices outervars)---- | Use indices with full variable set-forSubB :: DVSet -- ^ Inner vars -        -> InnerLoop a -- ^ Inner loop function-        -> DVSet -- ^ All vars-        -> [Int] -- ^ Outer indices-        -> [a]-forSubB innervars f allvars  outerIdx  = -        let sCode s e = if (e `isElem` s) then True else False-            selection s = map (sCode s) allvars-            computeInner iinner =-                let innerIdx = extend (selection innervars) 0 iinner-                    idx = zipWith (+) outerIdx innerIdx-                    in -                    f idx-        in-        map computeInner (forAllIndices innervars)- -- | Norm the factor normedFactor :: Factor f => f -> f  normedFactor f = factorDivide f (factorNorm f)@@ -218,7 +85,7 @@ -- | A factor as used in graphical model -- It may or not be a probability distribution. So it has no reason to be -- normalized to 1-class FactorPrivate f => Factor f where+class Factor f where     -- | When all variables of a factor have been summed out, we have a scalar     isScalarFactor :: f -> Bool      -- | An empty factor with no variable and no values@@ -226,20 +93,24 @@     -- | Check if a given discrete variable is contained in a factor     containsVariable :: f -> DV  -> Bool     -- | Give the set of discrete variables used by the factor-    factorVariables :: f -> DVSet    +    factorVariables :: f -> [DV]         -- | Return A in P(A | C D ...). It is making sense only if the factor is a conditional propbability     -- table. It must always be in the vertex corresponding to A in the bayesian graph     factorMainVariable :: f -> DV-    factorMainVariable = head . factorVariables+    factorMainVariable f = let vars = factorVariables f +      in+      case vars of +        [] -> error "Can't get the main variable of a scalar factor"+        (h:_) -> h      -- | Create a new factors with given set of variables and a list of value     -- for initialization. The creation may fail if the number of values is not     -- coherent with the variables and their levels.     -- For boolean variables ABC, the value must be given in order     -- FFF, FFT, FTF, FTT ...-    factorWithVariables :: DVSet -> [Double] -> Maybe f+    factorWithVariables :: [DV] -> [Double] -> Maybe f     -- | Value of factor for a given set of variable instantitation.     -- The variable instantion is like a multi-dimensional index.-    factorValue :: f -> DVISet Int -> Double+    factorValue :: f -> [DVI Int] -> Double     -- | Position of a discrete variable in te factor (p(AB) is differennt from p(BA) since values     -- are not organized in same order in memory)     variablePosition :: f -> DV -> Maybe Int@@ -258,107 +129,42 @@      -- | Create an evidence factor from an instantiation.     -- If the instantiation is empty then we get nothing-    evidenceFrom :: DVISet Int -> Maybe f+    evidenceFrom :: [DVI Int] -> Maybe f           -- | Divide all the factor values     factorDivide :: f -> Double -> f     factorDivide f d = (1.0 / d) `factorScale` f  +    factorToList :: f -> [Double]+     -- | Multiply factors.      factorProduct :: [f] -> f-    factorProduct [] = factorFromScalar 1.0-    factorProduct l = -        let allVars = L.foldl1' union . map factorVariables $ l-        in -        if L.null allVars -            then -                factorFromScalar (product . map factorNorm $ l) -            else-                let getFactorValueAtIndex i factor = factorValuePrivate factor (reorder i factor)-                    instantiationProduct instantiation = product . map (getFactorValueAtIndex instantiation) $ l-                    values = [instantiationProduct x | x <- forAllInstantiations allVars]-                in -                fromJust $ factorWithVariables allVars values      -- | Project out a factor. The variable in the DVSet are summed out-    factorProjectOut :: DVSet -> f -> f-    factorProjectOut s f = -        let alls = factorVariables f-            s' = alls `difference` s-        in -        if null s'-            then -                factorFromScalar (factorNorm f)-            else-                let dstValues = forSubA alls s' -                                   (forSubB s $ factorValuePrivate f)-                                   (\i c -> sum c)-                in -                fromJust $ factorWithVariables s' dstValues+    factorProjectOut :: [DV] -> f -> f+     -- | Project to. The variable are kept and other variables are removed-    factorProjectTo :: DVSet -> f -> f +    factorProjectTo :: [DV] -> f -> f      factorProjectTo s f =          let alls = factorVariables f              s' = alls `difference` s          in          factorProjectOut s' f --- | Used internaly when we know the position of a variable in the factor--- then we can identify the variable with an int. May be a bit faster for some--- algorithms-class FactorPrivate f where-    factorValuePrivate :: f -> [Int] -> Double---- | Return all the index (position in the factor) for a DV-allValues :: DV -> [Int]-allValues (DV _ i) = [0..i-1]---- | Generate all indexes for a set of variables-forAllIndices :: DVSet -> [[Int]]-forAllIndices = mapM allValues---- | Generate all instantiations of variables-forAllInstantiations :: DVSet -> [DVISet Int]-forAllInstantiations = mapM oneInstantiation- where-    oneInstantiation v@(DV vertex _) = map (setDVValue v) . allValues $ v- -- | Change the layout of values in the -- factor to correspond to a new variable order-changeVariableOrder :: DVSet -- ^ Old order-                    -> DVSet -- ^ New order +-- Used to import external files+changeVariableOrder :: DVSet s -- ^ Old order+                    -> DVSet s' -- ^ New order                      -> [Double] -- ^ Old values                     -> [Double] -- ^ New values-changeVariableOrder oldOrder newOrder oldValues =+changeVariableOrder (DVSet oldOrder) newOrder oldValues =     let oldFactor = fromJust $ factorWithVariables oldOrder oldValues :: CPT     in     [factorValue oldFactor i | i <- forAllInstantiations newOrder]  --- | Order the variable to get a multiindex which is--- making sense in the CPT. Only the subset in CPT is selectionned and reordered-reorder :: Factor f => DVISet Int -> f  -> [Int]-reorder i f = -    let nbDestVars = nbElements . factorVariables $ f-        v = V.replicate nbDestVars 0-        asDV v = DV v 0-        vectorPair bdvi = do -            pos <- variablePosition f . asDV . variableVertex $ bdvi-            let value = instantiationValue bdvi-            return (pos, value)-        allPos = mapMaybe vectorPair i-    in-    let testError = maybe False (const True) $ do -        guard $ length allPos == nbDestVars-        guard $ and . map ( (< nbDestVars) . fst)  $ allPos-        return ()-    in-    case testError of-      False -> error ("reorder has not set all destination indexes ! allpos = " ++ show allPos ++ " nbDestVars = " ++ show nbDestVars ++ "\n" ) -      True -> V.toList $ v V.// allPos-- -- | Mainly used for conditional probability table like p(A B | C D E) but the normalization to 1 -- is not imposed. And the conditionned variables are not different from the conditionning ones. -- The dimensions for each variables are listed.@@ -368,11 +174,11 @@ -- the knowledge of the dependencies is. -- So, this same structure is used for a probability too (conditioned on nothing) data CPT = CPT {-           dimensions :: DVSet -- ^ Dimensions for all variables-         , mapping :: IM.IntMap Int -- ^ Mapping from vertex number to position in dimensions-         , values :: V.Vector Double -- ^ Table of values+           dimensions :: ![DV] -- ^ Dimensions for all variables+         , mapping :: !(IM.IntMap Int) -- ^ Mapping from vertex number to position in dimensions+         , values :: !(V.Vector Double) -- ^ Table of values          }-         | Scalar Double+         | Scalar !Double  debugCPT (Scalar d) = do     putStrLn "SCALAR CPT"@@ -496,9 +302,11 @@ -- meaning of the variables according to their position. isomorphicFactor :: Factor f => f -> f -> Bool isomorphicFactor fa fb = maybe False (const True) $ do -    let va = factorVariables fa -        vb = factorVariables fb -    guard (va `equal` vb)+    let sa = factorVariables fa +        sb = factorVariables fb +        va = DVSet sa +        vb = DVSet sb+    guard (sa `equal` sb)     guard (factorDimension fa == factorDimension fb)     guard $ and [factorValue fa ia `nearlyEqual` factorValue fb ia | ia <- forAllInstantiations va]     return ()@@ -508,22 +316,25 @@ Following functions are used to typeset the factor when displaying it  -}-vname :: Int -> Int -> Box-vname vc i = text $ "v" ++ show vc ++ "=" ++ show i+-- | Display a variable name and its size+vname :: Int -> DVI Int -> Box+vname vc i = text $ "v" ++ show vc ++ "=" ++ show (instantiationValue i) -dispFactor :: FactorPrivate f => f -> DV -> [Int] -> DVSet -> Box+dispFactor :: Factor f => f -> DV -> [DVI Int] -> [DV] -> Box dispFactor cpt h c [] = -    let dstIndexes = allValues h+    let dstIndexes = allInstantiationsForOneVariable h         dependentIndexes =  reverse c         factorValueAtPosition p = -            let v = factorValuePrivate cpt p+            let v = factorValue cpt p             in             text . show  $ v     in     vsep 0 center1 . map (factorValueAtPosition . (:dependentIndexes)) $ dstIndexes  dispFactor cpt dst c (h@(DV (Vertex vc) i):l) = -    hsep 1 top . map (\i -> vcat center1 [vname vc i,dispFactor cpt dst (i:c) l])  $ (allValues h)+    let allInst = allInstantiationsForOneVariable h+    in+    hsep 1 top . map (\i -> vcat center1 [vname vc i,dispFactor cpt dst (i:c) l])  $ allInst  instance Show CPT where     show (Scalar v) = "\nScalar Factor:\n" ++ show v@@ -532,11 +343,17 @@     show c@(CPT d _ v) =          let h@(DV (Vertex vc) _) = head d             table = dispFactor c h [] (tail d)-            dstColumn = vcat center1 $ replicate (length d - 1) (text "") ++ map (vname vc) (allValues h)+            dstIndexes = map head (forAllInstantiations . DVSet $ [h])+            -- In P(A | B ...), the dst column is containing the possible values for the+            -- variables A with a header made of space to be aligned with the other part of the table.+            -- In the other part of the table, this header is containing the variable values for the other varibles+            dstColumn = vcat center1 $ replicate (length d - 1) (text "") ++ map (vname vc) dstIndexes         in         "\n" ++ show d ++ "\n" ++ render (hsep 1 top [dstColumn,table])  instance Factor CPT where+    factorToList (Scalar v) = [v]+    factorToList (CPT _ _ v) = V.toList v     emptyFactor = emptyCPT     isScalarFactor (Scalar _) = True     isScalarFactor _ = False@@ -548,55 +365,111 @@     factorWithVariables = createCPTWithDims     factorVariables (CPT v _ _) = v     factorVariables (Scalar _) = []-    factorNorm f@(CPT _ _ _) = sum [ factorValuePrivate f x | x <- forAllIndices (factorVariables f)]+    factorNorm f@(CPT d _ vals) = +        let vars = DVSet d+            strides = indexStrides vars+        in+        sum [ vals!(indexPosition strides x) | x <- indicesForDomain vars]     factorNorm (Scalar v) = v     variablePosition (CPT _ m _) (DV (Vertex i) _) = IM.lookup i m     variablePosition (Scalar _) _ = Nothing     factorScale s (Scalar v) = Scalar (s*v)-    factorScale s f = -        let newValues = map (s *) [ factorValuePrivate f x | x <- forAllIndices (factorVariables f)]+    factorScale s f@(CPT d _ vals) = +        let vars = DVSet d+            strides = indexStrides vars+            newValues = map (s *) [ vals!(indexPosition strides x) | x <- indicesForDomain vars]         in          fromJust $ factorWithVariables (factorVariables f) newValues     factorValue (Scalar v) _ = v -    factorValue f i = -        let multiIndex = reorder i f+    factorValue f@(CPT d _ v) i = +        let vars = DVSet d+            (dv,pos) = instantiationDetails i+            strides = indexStridesFor vars dv         in -        factorValuePrivate f multiIndex+        v!(indexPosition strides pos)     evidenceFrom [] = Nothing      evidenceFrom l = -        let index = map instantiationValue l -            variables = map instantiationVariable l+        let (variables,index) = instantiationDetails l+            DVSet nakedVars = variables             setValueForIndex i = if i == index then 1.0 else 0.0          in-        factorWithVariables variables . map setValueForIndex $ forAllIndices variables--instance FactorPrivate CPT where-    factorValuePrivate = getCPTValue+        factorWithVariables nakedVars . map setValueForIndex $ indicesForDomain variables+    factorProduct [] = factorFromScalar 1.0+    factorProduct l = +        let allVars = DVSet $ L.foldl1' union . map factorVariables $ l+            DVSet nakedVars = allVars+            (scalars,cpts) = partition isScalarFactor l+            stridesFromCPT (CPT d _ _) = indexStridesFor (DVSet d) allVars+            ps = product . map (flip factorValue undefined) $ scalars+            cptsStrides = map stridesFromCPT cpts+        in +        if L.null cpts +            then +                factorFromScalar ps+            else+                let getFactorValueAtIndex i (strides,factor@(CPT _ _ vals)) = vals!(indexPosition strides i)+                    instantiationProduct instantiation = product . map (getFactorValueAtIndex instantiation) $ (zip cptsStrides cpts)+                    values = [ps * instantiationProduct x | x <- indicesForDomain allVars]+                in +                values `seq` fromJust $ factorWithVariables nakedVars values+    factorProjectOut _ f@(Scalar v) = f+    factorProjectOut s f@(CPT d _ v) = +        let variablesToSum = s+            variablesToKeep = d `difference` s +            keepSet = DVSet variablesToKeep+            sumSet = DVSet variablesToSum +            strides = indexStridesFor (DVSet d) (DVSet $ variablesToKeep ++ variablesToSum) +            values = do +                  keepIndex <- indicesForDomain keepSet +                  let l = do+                        sumIndex <- indicesForDomain sumSet +                        return $ v!(indexPosition strides $ combineIndex strides keepIndex sumIndex)+                  return (sum l)+        in+        values `seq` fromJust $ factorWithVariables variablesToKeep values+        +-- | Used to combined the keep and sum indices in the factor projection+combineIndex :: Strides s'' -> [Index s] -> [Index s'] -> [Index s''] +combineIndex _ la lb = map (Index . fromIndex) la ++ map (Index .fromIndex) lb  -- | An empty CPT emptyCPT :: CPT emptyCPT = CPT [] IM.empty V.empty --- | Convertion of a multiindex to its--- position inside of the data vector of a 'CPT'-indexPosition :: DVSet -> [Int] -> Int-indexPosition [] _ = 0-indexPosition d pos = -    let dim = map dimension d+newtype Strides s = Strides [Int] deriving(Eq,Show)++-- | Generate strides to read the first CPT using an index having meaning in the second CPT+indexStridesFor :: DVSet s -- ^ DVSet to be read+                -> DVSet s' -- ^ DVSet to interpret the index+                -> Strides s'+indexStridesFor dr@(DVSet drvars) di@(DVSet divars) =+    let Strides originStrides = indexStrides dr+        reference = zip drvars originStrides +        getNewStrides dv = maybe 0 id (lookup dv reference)+    in +    Strides $ map getNewStrides divars+    ++-- | Generate the strides to read a given factor using a multiindex+-- using the same order as the factor variables+indexStrides :: DVSet s -> Strides s+indexStrides d@(DVSet dvars)  = +    let dim = map dimension dvars         pos' = scanr (*) (1::Int) (tail dim)-        c = sum . map (\(x,y) -> x * y) $ (zip pos' pos)     in -    c+    Strides pos'+-- | Convertion of a multiindex to its+-- position inside of the data vector of a 'CPT'+indexPosition :: Strides s -> [Index s] -> Int+{-# INLINE indexPosition #-}+indexPositions _ []  = 0+indexPosition (Strides pos') pos = sum . map (\(x,y) -> x * fromIndex y) $ (zip pos' pos) --- | Get the value at a given position. Positions are starting at zero-getCPTValue :: CPT -> [Int] -> Double-getCPTValue (Scalar v) _ = v-getCPTValue cpt@(CPT d _ v) pos = v!(indexPosition d pos)  -- | Create a CPT given some dimensions and a list of Doubles. -- Returns nothing is the length are not coherents.-createCPTWithDims :: DVSet -> [Double] -> Maybe CPT+createCPTWithDims :: [DV] -> [Double] -> Maybe CPT createCPTWithDims dims values =      let createDVIndex i (DV (Vertex v) _)  = (v,i)         m = IM.fromList . zipWith createDVIndex ([0,1..]::[Int]) $ dims

@@ -9,37 +9,43 @@     -- * Triangulation     , nodeComparisonForTriangulation     , numberOfAddedEdges+    , weight+    , weightedEdges     , triangulate     -- * Junction tree-    , minimumSpanningTree     , createClusterGraph     , Cluster     , createJunctionTree+    , createUninitializedJunctionTree     , JunctionTree     -- * Shenoy-Shafer message passing     , collect      , distribute     , posterior      -- * Evidence-    , clearEvidence-    , updateEvidence+    , changeEvidence     -- * Test      , junctionTreeProperty_prop-    , createVerticesJunctionTree     , VertexCluster+    -- * For debug +    , junctionTreeProperty+    , maximumSpanningTree+    , fromVertexCluster     ) where  import Bayes import qualified Data.Foldable as F import Data.Maybe(fromJust,mapMaybe,isJust)-import Control.Monad(mapM)+import Control.Monad(mapM,guard) import Bayes.Factor hiding (isEmpty) import Data.Function(on)-import Data.List(minimumBy,maximumBy,inits)+import Data.List(minimumBy,maximumBy,inits,foldl') import qualified Data.Set as Set import qualified Data.Map as Map import qualified Data.Functor as Functor import qualified Data.Tree as T +import Bayes.FactorElimination.JTree+import Control.Applicative((<$>))  import Test.QuickCheck hiding ((.||.), collect) import Test.QuickCheck.Arbitrary@@ -63,6 +69,15 @@     in      length [edge x y | x <- nodes, y <- nodes, x /= y, not (isLinkedWithAnEdge g x y)] +weightedEdges :: (UndirectedGraph g, Factor f) +              => g a f +              -> Vertex +              -> Int +weightedEdges g v = +    let nodes = fromJust $ neighbors g v+    in +    sum [weight g x * weight g y | x <- nodes, y <- nodes, x /= y, not (isLinkedWithAnEdge g x y)]+ -- | Weight of a node weight :: (UndirectedGraph g, Factor f)        => g a f @@ -101,6 +116,13 @@ instance Show VertexCluster where      show (VertexCluster s) = show . Set.toList $ s +++instance IsCluster Cluster where +  overlappingEvidence (Cluster c) e = filter (\x -> Set.member (instantiationVariable x) c) e+  clusterVariables (Cluster s) = Set.toList s+  mkSeparator (Cluster sa) (Cluster sb) = Cluster $ Set.intersection sa sb+ -- | Triangulate a graph using a cost function -- The result is the triangulated graph and the list of clusters -- which may not be maximal.@@ -153,155 +175,6 @@       (Nothing,_) -> checkIfMaximal (current:reversedPrefix) (head suffix) (tail suffix)       (Just r,l) -> checkIfMaximal (r:reverse l) (head suffix) (tail suffix) ---- | Create the cluster graph-createClusterGraph :: UndirectedGraph g-                   => [VertexCluster] -                   -> g Int VertexCluster-createClusterGraph c =-  let numberedClusters = zip c (map Vertex [0..])-      addCluster (c,v) g = addVertex v c g-      graphWithoutEdges = foldr addCluster emptyGraph numberedClusters-      separatorSize ca cb = Set.size $ Set.intersection (fromVertexCluster ca) (fromVertexCluster cb)-      allEdges = [(cx,cy) | cx <- numberedClusters, cy <- numberedClusters, cx /= cy]-      addClusterEdge ((ca,va),(cb,vb)) g = addEdge (edge va vb) (separatorSize ca cb) g-  in -  foldr addClusterEdge graphWithoutEdges allEdges---{---Minimum spanning tree using Prim's algorithm-  --}---- | Tree with values on edges-data Tree b a = Node a [(b,Tree b a)] deriving(Eq)--{---Implementation of show for the tree- --}-standardHaskellTree :: (Show f, Show b) => Tree b (JTNodeValue f) -> T.Tree String -standardHaskellTree n@(Node a []) = T.Node (show $ nodeCluster n) []-standardHaskellTree n@(Node a l) = T.Node (show $ nodeCluster n) (map (standardHaskellTree  . snd) l)--standardVertexTree :: Tree () VertexCluster -> T.Tree String -standardVertexTree n@(Node a []) = T.Node (show a) []-standardVertexTree n@(Node a l) = T.Node (show a) (map (standardVertexTree  . snd) l)-  -showFactorsAndEdges :: (Show f, Show b) => Tree b (JTNodeValue f) -> (String -> String) -showFactorsAndEdges  n@(Node a []) = (++ show (nodeValueFactor a))-showFactorsAndEdges  n@(Node a l) = foldl1 (.) (map factorAndEdge l) . (++ show (nodeValueFactor a)) -  where -    factorAndEdge (s,t) = showFactorsAndEdges t . (++ show s) --instance (Show f ,Show b)=> Show (Tree b (JTNodeValue f)) where -  show t = "JUNCTION TREE\n" ++ T.drawTree (standardHaskellTree t) ++ "\n" ++ showFactorsAndEdges t "" ++ "\n------\n"--instance Show (Tree () VertexCluster) where -  show t = "JUNCTION TREE\n" ++ T.drawTree (standardVertexTree t) ++ "\n"--instance Functor.Functor (Tree b) where -  fmap f (Node a l) = Node (f a) (map (mapEdge f) l)-    where -      mapEdge f (e,c) = (e, fmap f c)---- | Expand a tree (encoded as a list of edges)--- by adding vertices and keeping track of the vertices which have--- already been added.--- The selection of where to connect the new vertices is based upon cost of the new edges-expand :: UndirectedGraph g -       => g Int f -       -> [Edge] -- ^ List of edges-       -> [Vertex] -- ^ Vertices in Tree-       -> [Vertex] -- ^ Vertices to add-       -> [Edge] -- ^ Updated sets and edge list-expand g theEdges inTree remaining | null remaining = theEdges-                                   | otherwise = -                                        let (treeVertex,outVertex) = maximumBy (compare `on` (edgeCost g)) $ [(vin,vout) | vin <- inTree, vout <-remaining,isLinkedWithAnEdge g vin vout]-                                        in -                                        expand g (edge treeVertex outVertex : theEdges) (outVertex : inTree)-                                          (filter (/= outVertex) remaining)--  where -    edgeCost g (va,vb) = fromJust $ edgeValue g (edge va vb)--leaf x = Node x []-treeEdge c b = (c,b)---- | Create a tree based upon a description with edges-createTreeFromMap :: Vertex -- ^ Root vertex-                  -> Map.Map Vertex [Vertex] -- ^ Tree branches-                  -> Tree () Vertex -createTreeFromMap root m = -  let growTree m t@(Node a _) | Map.null m = t-                              | otherwise = -                                    case Map.lookup a m of -                                      Nothing -> t -                                      Just l -> Node a . map (treeEdge () . growTree m . leaf) $ l-  in-  growTree m (leaf root)-                   --- | Implementing the Prim's algorithm for minimum spanning tree-minimumSpanningTree :: UndirectedGraph g -                    => g Int f -                    -> Tree () f -minimumSpanningTree g = -  let startRoot = fromJust $ someVertex g -      remainingVertices = filter (/= startRoot) (allVertices g)-      foundEdges = expand g [] [startRoot] remainingVertices-      m = Map.fromListWith (++) . map ((\(a,b) -> (a,[b])) . edgeEndPoints) $ foundEdges-      theTree = createTreeFromMap startRoot m-  in -  Functor.fmap (fromJust . vertexValue g) theTree-      -   -{---Junction tree algorithm---}---- | Check if all variables of a factor are included in a cluster-vertexClusterIsContainingFactor :: Factor f => VertexCluster -> f -> Bool -vertexClusterIsContainingFactor c f = -  let factorVars = Set.fromList . map variableVertex . factorVariables $ f-  in -  Set.isSubsetOf factorVars (fromVertexCluster c)---- | Check if all variables of a factor are included in a cluster-clusterIsContainingVariable :: DV -> Cluster  -> Bool -clusterIsContainingVariable v c  =  -  Set.member v (Set.fromList $ fromCluster c)---- | Separator which can be in 3 state depending how many messages have passed through it-data Separator f = NoMessage !Cluster-                 | Collect !Cluster !f -                 | Distribute !Cluster !f !f -- Upward and downward message-                 deriving(Eq)--instance Show f => Show (Separator f) where -  show (NoMessage c) = "NoMessage: " ++ show c -  show (Collect c u) = "Collect: " ++ show c ++ "\n" ++ "\n <----- \n" ++ show u ++ "\n"-  show (Distribute c u d) = "Distribute: " ++ show c ++ "\n <----- \n" ++ show u ++ "\n" ++ " -----> \n" ++ show d ++ "\n"----- | Evidence if some is used for the node-type Evidence f = f---- | Evidence for cluster, factor for cluster-data JTNodeValue f = JTNodeValue !Cluster !(Evidence f) !f deriving(Eq,Show)---- | Cluster of discrete variables.--- Discrete variables instead of vertices are needed because the--- factor are using 'DV' and we need to find--- which factors must be contained in a given cluster.-newtype Cluster = Cluster (Set.Set DV) deriving(Eq,Show)--fromCluster (Cluster s) = Set.toList s - -- | Convert the clusters from vertex to 'DV' clusters vertexClusterToCluster :: (Factor f , Graph g)                        => g e f @@ -313,150 +186,115 @@   in    Cluster . Set.fromList $ variables --- | Vertices contained in a cluster-clusterVertices :: VertexCluster -> [Vertex]-clusterVertices = Set.toList . fromVertexCluster --- | Find all factors contained in a cluster-findFactorsForCluster :: (Factor f , Graph g)-                      => BayesianNetwork g f-                      -> VertexCluster-                      -> [f]-findFactorsForCluster g c = -  filter (vertexClusterIsContainingFactor c) . mapMaybe (vertexValue g) . clusterVertices $ c---- | The junction tree-type JunctionTree f = Tree (Separator f) (JTNodeValue f)---- | Get the potential for a cluster-mkNodePotential :: (Graph g, Factor f, Show f)-                => BayesianNetwork g f -                -> VertexCluster -                -> Set.Set Vertex-                -> (JTNodeValue f, Set.Set Vertex)-mkNodePotential g c set =  -  let -- Factor found in a cluster but they may already be used in another cluster-      foundFactors = findFactorsForCluster g c-      -- Get the vertices for the factor-      vertexForFactors = map (variableVertex . factorMainVariable) foundFactors -      -- Keep only the factors which are not already used-      isNotUsed (v,f) = Set.member v set-      factorsNotYetUsed = filter isNotUsed (zip vertexForFactors foundFactors)-      set' = Set.difference set (Set.fromList $ map fst factorsNotYetUsed)-      factorsToUse = map snd factorsNotYetUsed-    -      potential = factorProduct factorsToUse-  in -  (JTNodeValue (vertexClusterToCluster g c) (factorFromScalar 1.0) potential, set')---- | Generate the evidence potential for a given cluster-evidenceForCluster :: Factor f => DVISet Int -> Cluster -> Maybe (Evidence f)-evidenceForCluster assignments cluster@(Cluster c) = -  let c' = Set.fromList (map instantiationVariable assignments) -      common = Set.intersection c' c -      selectedVariables = filter (\c -> Set.member (instantiationVariable c) common) assignments+-- | Create the cluster graph+createClusterGraph :: (UndirectedGraph g, Factor f, Graph g')+                   => g' e f+                   -> [VertexCluster] +                   -> g Int Cluster+createClusterGraph bn c =+  let numberedClusters = zip c (map Vertex [0..])+      addCluster g (c,v)  = addVertex v (vertexClusterToCluster bn c) g+      graphWithoutEdges = foldl' addCluster emptyGraph numberedClusters+      separatorSize ca cb = Set.size $ Set.intersection (fromVertexCluster ca) (fromVertexCluster cb)+      allEdges = [(cx,cy) | cx <- numberedClusters, cy <- numberedClusters, cx /= cy]+      addClusterEdge g ((ca,va),(cb,vb)) = addEdge (edge va vb) (separatorSize ca cb) g   in -  evidenceFrom selectedVariables+  foldl' addClusterEdge graphWithoutEdges allEdges  --- | Get the cluster for a node-nodeCluster :: Tree a (JTNodeValue f) -> Cluster -nodeCluster (Node (JTNodeValue c _ _ ) _) = c --emptyCluster :: Cluster -emptyCluster = Cluster Set.empty--nodeValueFactor (JTNodeValue _ _ f ) = f-nodeValueEvidence (JTNodeValue _ e _) = e--nodeValueWithNewEvidence (JTNodeValue a e b) e' = JTNodeValue a e' b-clearNodeValueEvidence (JTNodeValue a _ b)  = JTNodeValue a (factorFromScalar 1.0) b---- | Get the cluster for a separator-separatorCluster :: Separator f -> Cluster -separatorCluster (NoMessage c) = c-separatorCluster (Collect c _) = c -separatorCluster (Distribute c _ _) = c -+{- -upMessage (Distribute _ u _) = Just u -upMessage (Collect _ u ) = Just u -upMessage _ = Nothing +Maximum spanning tree using Prim's algorithm+  +-} -downMessage (Distribute _ _ d) = Just d -downMessage _ = Nothing +-- | Get all possible edges between the leaves and the remaining nodes+possibilities :: (Ord c , UndirectedGraph g) +              => g Int c -- ^ Original graph to get the edge value +              -> JTree c (Vertex,f) -- ^ Tree to get the vertex for a leaf+              -> [Vertex] -- ^ Vertices to add to the tree+              -> [c] -- ^ List of leaves+              -> [(Vertex,c,Int)] -- ^ Found edge to add+possibilities g currentT remaining leavesClusters = do +  rv <- remaining+  lv <- leavesClusters+  let NodeValue (lvVertex,lvCluster) _ = nodeValue currentT lv+  guard (isLinkedWithAnEdge g rv lvVertex)+  let ev = fromJust $ edgeValue g (edge rv lvVertex)+  return $ (rv,lv,ev) -computeSeparatorCluster :: (Factor f, Graph g) -                        => BayesianNetwork g f -                        -> VertexCluster -                        -> VertexCluster-                        -> Cluster-computeSeparatorCluster g parent child = -  let theNodeCluster (Node c _) = c -      childVertices = fromVertexCluster child -      parentVertices = fromVertexCluster parent-      separatorVertices = VertexCluster $ Set.intersection childVertices parentVertices+-- | Find the max edge to add to the tree+findMax :: (UndirectedGraph g, Ord c, Factor f)+        => g Int c -- ^ Graph+        -> [Vertex] -- ^ Nodes to add +        -> JTree c (Vertex,f)+        -> ([Vertex],(Vertex,c),c) +findMax g remaining currentT = +  let leavesClusters = treeNodes currentT+      edgeValue (_,_,e) = e+      (rf,lf,ef) = maximumBy (compare `on` edgeValue) (possibilities g currentT remaining leavesClusters)+      remaining' = filter (/= rf) remaining +      foundCluster = fromJust $ vertexValue g rf   in-  vertexClusterToCluster g  separatorVertices+  (remaining', (rf, foundCluster), lf) -dfs :: (n -> n -> e -> e') -- Parent, child node and their egde-    -> (n -> a -> (n', a)) -- Node and current value -> new value and new nod-    -> Tree e n  -- Tree to traverse-    -> a -- Start value-    -> (Tree e' n', a) -- New tree and new value-dfs edgef nodef n@(Node nodevalue []) current = -  let (newnodevalue, newval) = nodef nodevalue current -  in -  (Node newnodevalue [],newval) -dfs edgef nodef n@(Node nodevalue children) current =-  let (newnodevalue, newval) = nodef nodevalue current -      applyEdgeFunction (e,Node childvalue _) = edgef nodevalue childvalue e-      applyToChildren childrenNode val = dfs edgef nodef childrenNode val-      edges' = map applyEdgeFunction children-      recurseOnChildren s r [] = (s,reverse r)-      recurseOnChildren s r (a:l) = -        let (a',s') = applyToChildren a s-        in -        recurseOnChildren s' (a':r) l -      (lastval,newSubTrees) = recurseOnChildren newval [] (map snd children)-  in -  (Node newnodevalue (zip edges' newSubTrees),lastval)+removeVertices :: JTree c (Vertex,f) -> JTree c f+removeVertices t = t { nodeValueMap = Map.map removeVertexFromNode (nodeValueMap t)+                     , separatorValueMap = Map.map removeVertexFromSeparator (separatorValueMap t)+                     }+ where +  removeVertexFromNode (NodeValue (_,f) (_,e)) = NodeValue f e +  removeVertexFromSeparator (SeparatorValue (_,u) (Just (_,d))) = SeparatorValue u (Just d)+  removeVertexFromSeparator (SeparatorValue (_,u) Nothing) = SeparatorValue u Nothing +  removeVertexFromSeparator EmptySeparator = EmptySeparator -setFactorEdgeUpdate :: (Graph g, Factor f) -                    => BayesianNetwork g f -                    -> VertexCluster -                    -> VertexCluster-                    -> () -                    -> Separator f-setFactorEdgeUpdate g parentvalue childvalue _ = NoMessage $ computeSeparatorCluster g parentvalue childvalue +-- | Implementing the Prim's algorithm for minimum spanning tree+maximumSpanningTree :: (UndirectedGraph g, IsCluster c, Factor f, Ord c) +                    => g Int c +                    -> JTree c f+maximumSpanningTree g = +    let rootNodeVertex = fromJust $ someVertex g +        rootNodeValue = fromJust $ vertexValue g rootNodeVertex+        unitFactor = factorFromScalar 1.0 +        startTree = singletonTree rootNodeValue (rootNodeVertex,unitFactor) (rootNodeVertex,unitFactor) +        remainingVertices = filter (/= rootNodeVertex) (allVertices g) +    in +    removeVertices $ buildTree g remainingVertices startTree  -setFactorNodeUpdate :: (Graph g, Factor f, Show f) -                    => BayesianNetwork g f -                    -> VertexCluster-                    -> Set.Set Vertex -                    -> (JTNodeValue f, Set.Set Vertex)-setFactorNodeUpdate g nodeValue set = mkNodePotential g nodeValue set+buildTree :: (UndirectedGraph g , IsCluster c, Factor f, Ord c)+          => g Int c +          -> [Vertex]+          -> JTree c (Vertex,f) +          -> JTree c (Vertex,f) +buildTree g [] currentT = currentT +buildTree g l@(h:t) currentT = +    let unitFactor = factorFromScalar 1.0+        (l',(foundElemVertex,foundElemValue),leaf) = findMax g l currentT+        sep = mkSeparator foundElemValue leaf+        newTree = addSeparator leaf sep foundElemValue . +                  addNode foundElemValue (foundElemVertex,unitFactor) (foundElemVertex,unitFactor) $ currentT+    in +    buildTree g l' newTree+   +{- --- | Set a factor for a node-setFactors :: (Graph g, Factor f, Show f)-           => BayesianNetwork g f -- ^ Bayesian graph-           -> Tree () VertexCluster  -- ^ Cluster tree with no factors-           -> Set.Set Vertex-           -> (JunctionTree f,Set.Set Vertex) -- ^ Initialized junction tree-setFactors g = dfs (setFactorEdgeUpdate g) (setFactorNodeUpdate g) +Junction tree algorithm +-} + -- | Create a junction tree with only the clusters and no factors-createVerticesJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g)-                           => (UndirectedSG () b -> Vertex -> Vertex -> Ordering) -- ^ Weight function on the moral graph-                           -> g () b -- ^ Input directed graph-                           -> Tree () VertexCluster -- ^ Junction tree-createVerticesJunctionTree cmp g =  +createUninitializedJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f)+                                => (UndirectedSG () f -> Vertex -> Vertex -> Ordering) -- ^ Weight function on the moral graph+                                -> g () f -- ^ Input directed graph+                                -> JunctionTree f -- ^ Junction tree+createUninitializedJunctionTree cmp g =     let theMoralGraph = moralGraph g       (clusters,_) = triangulate (cmp theMoralGraph) theMoralGraph-      g'' = createClusterGraph clusters :: UndirectedSG Int VertexCluster+      g'' = createClusterGraph g clusters :: UndirectedSG Int Cluster   in -  minimumSpanningTree g''+  maximumSpanningTree g''  -- | Create a function tree createJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f, Show f)@@ -464,137 +302,63 @@                   -> BayesianNetwork g f -- ^ Input directed graph                   -> JunctionTree f -- ^ Junction tree createJunctionTree cmp g = -  let cTree = createVerticesJunctionTree cmp g -      factorSet = Set.fromList (allVertices g) -- Tracking of factors which have not yet been put in the junction tree+  let cTree = createUninitializedJunctionTree cmp g        -- A vertex is linked with a factor so vertex is used as the identifier-      (newTree, _) = setFactors g cTree factorSet-  in -  distribute Nothing . collect $ newTree---collectMessages :: Factor f => (Separator f , JunctionTree f) -> (Separator f , JunctionTree f)-collectMessages (separator, Node nc []) = -  let sc = separatorCluster separator-      newPotential = factorProduct [nodeValueFactor nc,nodeValueEvidence nc] -      newMessage = factorProjectTo (fromCluster sc) newPotential-  in-  (Collect sc newMessage, Node nc []) -- Copy node factor to node current potential-collectMessages (separator,(Node nc l)) = -  let sc = separatorCluster separator-      messagesFromSubTrees = map collectMessages l -      newPotential = factorProduct (nodeValueEvidence nc:nodeValueFactor nc:(mapMaybe (upMessage . fst) messagesFromSubTrees))-      newMessage = factorProjectTo (fromCluster sc) newPotential -  in -  (Collect sc newMessage, Node nc messagesFromSubTrees)---- | Collect phase of the junction tree-collect :: Factor f => JunctionTree f -> JunctionTree f -collect t = let (_,t') = collectMessages (NoMessage emptyCluster, t) in t'--notSameCluster a b = nodeCluster a /= nodeCluster b ---- | Distribute phase of the junction tree-distribute :: Factor f => Maybe (Separator f) -> JunctionTree f -> JunctionTree f -distribute down n@(Node nc []) = n-distribute down (Node nc l) = -  let receivedDownMessage = if isJust down then fromJust . downMessage . fromJust $ down else factorFromScalar 1.0-      getUpMessage (edge,c) = upMessage edge -      upMessagesForSendingTo i = fromJust . mapM getUpMessage . filter ((i `notSameCluster`) . snd) $ l-      newPotential i = factorProduct (nodeValueFactor nc:nodeValueEvidence nc:receivedDownMessage:upMessagesForSendingTo i)-      newMessage sc i = factorProjectTo (fromCluster sc) (newPotential i)-      distributeMessage s@(Collect sc dm,i) = -        let newSeparator = Distribute sc dm (newMessage sc i)-        in -        (newSeparator,distribute (Just newSeparator) i)-      distributeMessage _ = error "Distribute message can only update a collect phase message"-      subTrees = map distributeMessage l+      newTree = setFactors g cTree   in -  Node nc subTrees---- | Depth first search in  tree-findInTree :: (Tree edge a -> Bool) -> Maybe edge -> Tree edge a -> Maybe (Maybe edge,Tree edge a)-findInTree cmp e n@(Node a []) = if (cmp n) then Just (e,n) else Nothing -findInTree cmp e n@(Node a l) = -  let findSome [] = Nothing-      findSome ((e',h):t) = -        case findInTree cmp (Just e') h of -          Nothing -> findSome t -          Just r -> Just r-  in-  case cmp n of -    True -> Just (e,n) -    False -> findSome l+  distribute . collect $ newTree   -- | Compute the marginal posterior (if some evidence is set on the junction tree)   -- otherwise compute just the marginal prior. posterior :: Factor f => JunctionTree f -> DV -> Maybe f-posterior t v = do -  (maybeEdge,Node n l) <- findInTree (clusterIsContainingVariable v . nodeCluster) Nothing t-  let receivedDownMessage = maybe (factorFromScalar 1.0) id $ -                               do-                                 e <- maybeEdge -                                 downMessage e-      upMessages = fromJust . mapM (upMessage . fst) $ l-      p = factorProduct (receivedDownMessage:nodeValueEvidence n:nodeValueFactor n:upMessages)-  return $ normedFactor $ factorProjectTo [v] p +posterior t v = +  case snd $ traverseTree (findClusterFor v) Nothing t of +    Nothing -> Nothing+    Just c -> let NodeValue f e = nodeValue t c +                  d = maybe (factorFromScalar 1.0) id $ downMessage t =<< (nodeParent t c)+                  u = map (upMessage t) (nodeChildren t c)+                  unNormalized = factorProjectTo [v] (factorProduct (f:e:d:u))+              in +              Just $ factorDivide unNormalized (factorNorm unNormalized) --- | Apply some evidence modifications in the tree-applyEvidenceWith :: (JunctionTree f -> JunctionTree f) -- ^ Node modification function. Only change node value. Not the children-                  -> JunctionTree f -- ^ Input tree-                  -> JunctionTree f-applyEvidenceWith nodeChange n@(Node _ []) = nodeChange n -applyEvidenceWith nodeChange n@(Node _ l) =-  let Node n' l' = nodeChange n -      changeChildren (e,c) = (e,applyEvidenceWith nodeChange c)-  in -  Node n' (map changeChildren l')+-- | Find a cluster containing the variable+findClusterFor :: DV +               -> Maybe Cluster+               -> Cluster -- ^ Current cluster+               -> NodeValue f -- ^ Current value+               -> Action (Maybe Cluster) (NodeValue f)+findClusterFor dv s c@(Cluster sc) v = +  case Set.member dv sc of +    False -> Skip s +    True -> Stop (Just c) --- | Change the evidence for a node-evidenceWith :: Factor f => DVISet Int -> JunctionTree f -> JunctionTree f-evidenceWith assignments t@(Node n l) = -  let n' = case evidenceForCluster assignments (nodeCluster t) of -             Nothing -> n -             Just e' -> nodeValueWithNewEvidence n e'-  in -  Node n' l --- | Remove the evidence for a node-clearNodeEvidence (Node n l) = Node (clearNodeValueEvidence n) l ---- | Remove evidence in the junction tree-clearEvidence :: Factor f => JunctionTree f -> JunctionTree f-clearEvidence = distribute Nothing . collect . applyEvidenceWith (clearNodeEvidence)---- | Update evidence in the tree-updateEvidence :: Factor f => DVISet Int -> JunctionTree f -> JunctionTree f-updateEvidence assignments = distribute Nothing . collect . applyEvidenceWith (evidenceWith assignments)---- | Used to implement quickcheck.--- The junction tree property is the property that CA intersection CB is included in all clusters in the path--- from CA to CB.-junctionTreeProperty :: [VertexCluster] -> Tree () VertexCluster -> Bool-junctionTreeProperty path (Node _ []) = True -junctionTreeProperty path (Node c l) = -  let children = map snd l -  in-  checkPath c (reverse path) && all (junctionTreeProperty (c:path)) children --junctionTreeProperty_prop :: DirectedSG () String -> Property +junctionTreeProperty_prop :: DirectedSG () CPT -> Property  junctionTreeProperty_prop g = (not . isEmpty) g && (not . hasNoEdges) g && connectedGraph g ==>    let cmp ug = (compare `on` (numberOfAddedEdges ug))+      t = createUninitializedJunctionTree cmp g   in-  junctionTreeProperty [] (createVerticesJunctionTree cmp g)+  junctionTreeProperty t [] (root t) --- | Check that the intersection of C with any parent in included in any cluster between the parent and C.-checkPath :: VertexCluster -> [VertexCluster] -> Bool -checkPath c l = -  let parentSets = map fromVertexCluster l-      allIntersections = map (Set.intersection (fromVertexCluster c)) parentSets-      pathsToEachParent = tail . inits $ parentSets-      isSubsetOfAllParents i parents = all (Set.isSubsetOf i) parents+junctionTreeProperty :: JTree Cluster CPT -> [Cluster] -> Cluster -> Bool+junctionTreeProperty t path c = +  let cl = map (separatorChild t) . nodeChildren t $ c+  in+  checkPath c path && all (junctionTreeProperty t (c:path)) cl +++-- | Check that the intersection of C with any parent in included in all cluster between the parent and C.+checkPath :: Cluster -> [Cluster] -> Bool +checkPath _ [] = True+checkPath (Cluster c) l = +  let clusterSet (Cluster s) = s -- x+      parentSets = map clusterSet l -- Example a b c d where a is the root+      allIntersectionsWithParents = map (Set.intersection c) parentSets -- a ^ x, b ^ x , c ^ x , d ^ x+      pathsToEachParent = tail . inits $ parentSets -- a, ab, abc, abcd+      isSubsetOfAllParents i path = all (Set.isSubsetOf i) path   in    -  and $ zipWith isSubsetOfAllParents allIntersections pathsToEachParent+  and $ zipWith isSubsetOfAllParents allIntersectionsWithParents pathsToEachParent {-  Moral graph@@ -625,8 +389,8 @@     (originGraph',dstGraph')  -- | Add the missing parent links-addMissingLinks :: DirectedGraph g => Vertex -> b -> g () b -> g () b-addMissingLinks v _ g = +addMissingLinks :: DirectedGraph g => g () b -> Vertex -> b -> g () b+addMissingLinks g v _ =      let (_,g') = connectAllNodesWith g g (\e m -> m) (\e m -> addEdge e () m) (parents g v)     in      g'@@ -637,11 +401,11 @@                     => g  () b                      -> g' () b  convertToUndirected m = -    let addVertexWithLabel v dat g = +    let addVertexWithLabel g v dat  =             let theName = fromJust $ vertexLabel m v            in             addLabeledVertex theName v dat g-        newDiscreteGraph = foldrWithVertex addVertexWithLabel emptyGraph m+        newDiscreteGraph = foldlWithVertex' addVertexWithLabel emptyGraph m         addEmptyEdge edge g = addEdge edge () g     in      foldr addEmptyEdge newDiscreteGraph . allEdges $ m@@ -651,4 +415,4 @@ moralGraph :: (NamedGraph g, FoldableWithVertex g, DirectedGraph g)             => g () b -> UndirectedSG () b  moralGraph g = -    convertToUndirected  . foldrWithVertex addMissingLinks g $ g+    convertToUndirected  . foldlWithVertex' addMissingLinks g $ g

@@ -0,0 +1,486 @@+{- | Junction Trees ++The Tree data structures are not working very well with message passing algorithms. So, junction trees are using+a different representation++-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+module Bayes.FactorElimination.JTree(+      IsCluster(..)+    , Cluster(..)+    , JTree(..)+    , JunctionTree(..)+    , setFactors+    , distribute +    , collect+    , fromCluster+    , changeEvidence+    , nodeIsMemberOfTree+    , singletonTree+    , addNode +    , addSeparator+    , leaves+    , nodeValue+    , NodeValue(..)+    , SeparatorValue(..)+    , downMessage+    , upMessage +    , nodeParent +    , nodeChildren+    , traverseTree+    , separatorChild+    , treeNodes+    , Action(..)+    ) where ++import qualified Data.Map as Map+import qualified Data.Tree as Tree+import Data.Maybe(fromJust,mapMaybe)+import qualified Data.Set as Set+import Data.Monoid+import Data.List((\\), intersect,partition, foldl')+import Bayes.PrivateTypes +import Bayes.Factor+import Bayes++import Debug.Trace +debug s a = trace (s ++ " " ++ show a ++ "\n") a++type UpMessage a = a +type DownMessage a = Maybe a++-- | Separator value+data SeparatorValue a = SeparatorValue !(UpMessage a) !(DownMessage a)+                      | EmptySeparator -- ^ Use to track the progress in the collect phase+                      deriving(Eq)++instance Show a => Show (SeparatorValue a) where +    show EmptySeparator = ""+    show (SeparatorValue u Nothing) = "u(" ++ show u ++ ")"+    show (SeparatorValue u (Just d)) = "u(" ++ show u ++ ") d(" ++ show d ++ ")"++type FactorValue a = a +type EvidenceValue a = a++-- | Node value+data NodeValue a = NodeValue !(FactorValue a) !(EvidenceValue a) deriving(Eq)++instance Show a => Show (NodeValue a) where +    show (NodeValue f e) = "f(" ++ show f ++ ") e(" ++ show e ++ ")"++-- | Junction tree.+-- 'c' is the node / separator identifier (for instance a set of 'DV')+-- a are the values for a node or separator+data JTree  c f = JTree {  root :: !c+                        -- | Leaves of the tree+                        ,  leavesSet :: !(Set.Set c)+                        -- | The children of a node are separators+                        ,  childrenMap :: !(Map.Map c [c])+                        -- | Parent of a node+                        ,  parentMap :: !(Map.Map c c)+                        -- | Parent of a separator+                        ,  separatorParentMap :: !(Map.Map c c)+                        -- | The child of a seperator is a node+                        ,  separatorChildMap :: !(Map.Map c c)+                        -- | Values for nodes and seperators+                        ,  nodeValueMap :: !(Map.Map c (NodeValue f))+                        ,  separatorValueMap :: !(Map.Map c (SeparatorValue f))+                        } deriving(Eq)++-- | Create a singleton tree with just one root node+singletonTree r factorValue evidenceValue = +    let t = JTree r Set.empty Map.empty Map.empty Map.empty Map.empty Map.empty Map.empty+    in +    addNode r factorValue evidenceValue t++-- | Leaves of the tree+leaves :: JTree c a -> [c]+leaves = Set.toList . leavesSet++-- | All nodes of the tree+treeNodes :: JTree c a -> [c]+treeNodes = Map.keys . nodeValueMap++-- | Value of a node+nodeValue :: Ord c => JTree c a -> c -> NodeValue a +nodeValue t e = fromJust $ Map.lookup e (nodeValueMap t)++-- | Change the value of a node+setNodeValue :: Ord c => c -> NodeValue a -> JTree c a -> JTree c a+setNodeValue c v t = t {nodeValueMap = Map.insert c v (nodeValueMap t)} ++-- | Parent of a node+nodeParent :: Ord c => JTree c a -> c -> Maybe c +nodeParent t e = Map.lookup e (parentMap t)++-- | Value of a node+separatorValue :: Ord c => JTree c a -> c -> SeparatorValue a +separatorValue t e = fromJust $ Map.lookup e (separatorValueMap t)++-- | Parent of a separator+separatorParent :: Ord c => JTree c a -> c -> c +separatorParent t e = fromJust $ Map.lookup e (separatorParentMap t)++-- | UpMessage for a separator node+upMessage :: Ord c => JTree c a -> c -> a+upMessage t c = case separatorValue t c of +                  SeparatorValue up _ -> up +                  _ -> error "Trying to get an up message on an empty seperator ! Should never occur !"++-- | DownMessage for a separator node+downMessage :: Ord c => JTree c a -> c -> Maybe a +downMessage t c = case separatorValue t c of +     SeparatorValue _ (Just down) -> Just down +     SeparatorValue _ Nothing -> Nothing+     _ -> error "Trying to get a down message on an empty separator ! Should never occur !"++-- | Return the separator childrens of a node+nodeChildren :: Ord c => JTree c a -> c -> [c]+nodeChildren t e = maybe [] id $ Map.lookup e (childrenMap t)++-- | Return the child of a separator+separatorChild :: Ord c => JTree c a -> c -> c +separatorChild t e = fromJust $ Map.lookup e (separatorChildMap t)++-- | Check if a node is member of the tree+nodeIsMemberOfTree :: Ord c => c -> JTree c a -> Bool +nodeIsMemberOfTree c t = Map.member c (nodeValueMap t)++-- | Add a separator between two nodes.+-- The nodes MUST already be in the tree+addSeparator :: (Ord c) +             => c -- ^ Origin node +             -> c -- ^ Separator+             -> c -- ^ Destination node +             -> JTree c a -- ^ Current tree +             -> JTree c a -- ^ Modified tree +addSeparator node sep dest t = +    t { childrenMap = Map.insertWith' (++) node [sep] (childrenMap t)+      , separatorChildMap = Map.insert sep dest (separatorChildMap t)+      , separatorValueMap = Map.insert sep EmptySeparator (separatorValueMap t)+      , leavesSet = Set.delete node (leavesSet t) +      , parentMap = Map.insert dest sep (parentMap t)+      , separatorParentMap = Map.insert sep node (separatorParentMap t)+      }++-- | Add a new node+addNode :: (Ord c) +        => c -- ^ Node+        -> a -- ^ Factor value +        -> a -- ^ Evidence value+        -> JTree c a +        -> JTree c a +addNode node factorValue evidenceValue t = +   t { nodeValueMap = Map.insert node (NodeValue factorValue evidenceValue) (nodeValueMap t)+     , leavesSet = Set.insert node (leavesSet t) +     }++-- | Update the up message of a separator+updateUpMessage :: Ord c +                => Maybe c -- ^ Separator node to update (if any : none for root node)+                -> a -- ^ New value+                -> JTree c a -- ^ Old tree+                -> JTree c a+updateUpMessage Nothing _ t = t+updateUpMessage (Just sep) newval t = +    let newSepValue =  case separatorValue t sep of +                         EmptySeparator -> SeparatorValue newval Nothing+                         SeparatorValue up down -> SeparatorValue newval down +    in +    t {separatorValueMap = Map.insert sep newSepValue (separatorValueMap t)}++-- | Update the down message of a separator+updateDownMessage :: Ord c +                  => c -- ^ Separator node to update+                  -> a -- ^ New value+                  -> JTree c a -- ^ Old tree+                  -> JTree c a+updateDownMessage sep newval t = +    let newSepValue = case separatorValue t sep of +                        EmptySeparator -> error "Can't set a down message on an empty separator"+                        SeparatorValue up _ -> SeparatorValue up (Just newval)+    in +    t {separatorValueMap = Map.insert sep newSepValue (separatorValueMap t)}++{-++Message passing algorithms+    +-}++-- | Functions used to generate new messages+class Message f c | f -> c where+    -- | Generate a new message from the received ones+    newMessage :: [f] -> NodeValue f -> c -> f +++-- | Check that a separator is initialized+separatorInitialized :: SeparatorValue a -> Bool+separatorInitialized EmptySeparator = False +separatorInitialized _ = True++allSeparatorsHaveReceivedAMessage :: Ord c+                                  => JTree c a -- ^ Tree+                                  -> [c] -- ^ Separators+                                  -> Bool +allSeparatorsHaveReceivedAMessage t seps = +  all separatorInitialized . map (separatorValue t) $ seps++-- | Update the up separator by sending a message+-- But only if all the down separators have received a message+updateUpSeparator :: (Message a c, Ord c) +                  => JTree c a +                  -> c -- ^ Node generating the new upMessage+                  -> JTree c a +updateUpSeparator t h  = +    let seps = nodeChildren t h+    in+    case allSeparatorsHaveReceivedAMessage t seps of +      False -> t +      True -> let incomingMessages = map (upMessage t) seps+                  currentValue = nodeValue t h+                  destinationNode = nodeParent t h+              in +              case destinationNode of +                Nothing -> t -- When root+                Just p -> let generatedMessage = newMessage incomingMessages currentValue p+                          in +                          updateUpMessage destinationNode generatedMessage t++-- | Update the down separator by sending a message+updateDownSeparator :: (Message a c, Ord c) +                    => c -- ^ Node generating the message +                    -> JTree c a +                    -> c -- ^ Child receiving the message+                    -> JTree c a +updateDownSeparator node t child  = +    let incomingMessagesFromBelow = map (upMessage t) (nodeChildren t node \\ [child])+        messageFromAbove = downMessage t =<< (nodeParent t node)+        incomingMessages = maybe incomingMessagesFromBelow (\x -> x:incomingMessagesFromBelow) messageFromAbove+        currentValue = nodeValue t node+        generatedMessage = newMessage incomingMessages currentValue child+    in +    updateDownMessage child generatedMessage t++unique :: Ord c => [c] -> [c]+unique = Set.toList . Set.fromList++data TraversalState = ACluster | ASeparator++-- | Collect message taking into account that the tree depth may be different for different leaves.+collect :: (Ord c, Message a c) +        => JTree c a +        -> JTree c a+collect t = _collect ACluster (leaves t) t++_collect :: (Ord c, Message a c) +         => TraversalState -- ^ Node processing phase or separator processing phase+         -> [c]+         -> JTree c a -- ^ Tree+         -> JTree c a -- ^ Modified tree +_collect _ [] t = t+_collect ACluster l t = +    let newTree = foldl' updateUpSeparator t l+    in +    _collect ASeparator (mapMaybe (nodeParent t) l) newTree+_collect ASeparator l t = _collect ACluster (unique . map (separatorParent t) $ l) t+    +distribute :: (Ord c, Message a c)+           => JTree c a +           -> JTree c a+distribute t = _distribute ACluster t (root t) ++_distribute :: (Ord c, Message a c)+            => TraversalState -- ^ True if node+            -> JTree c a +            -> c -- ^ Destination of the distribute+            -> JTree c a +_distribute ACluster t node = +    let children = nodeChildren t node+        newTree = foldl' (updateDownSeparator node) t $ children+    in+    foldl' (_distribute ASeparator) newTree children+_distribute ASeparator t node = _distribute ACluster t (separatorChild t node)++{-++Factors and evidence modifications++-}++-- | This class is used to check if evidence or a factor is relevant+-- for a cluster+class IsCluster c where +  overlappingEvidence :: c -> [DVI Int] -> [DVI Int]+  clusterVariables :: c -> [DV]+  mkSeparator :: c -> c -> c++instance IsCluster [DV] where +  overlappingEvidence c e = filter (\x -> instantiationVariable x `elem` c) e+  clusterVariables = id+  mkSeparator = intersect++data Action s a = Skip !s +                | ModifyAndStop !s !a+                | Modify !s !a+                | Stop !s++-- | Traverse a tree and modify it+traverseTree :: Ord c +             => (s -> c -> NodeValue f -> Action s (NodeValue f)) -- ^ Modification function+             -> s -- ^ Current state+             -> JTree c f -- ^ Input tree+             -> (JTree c f,s)+traverseTree action state t = _traverseTree True action (t,state) (root t)++_traverseTree False action (t,state) current  = _traverseTree True action (t,state)  (separatorChild t current) +_traverseTree True action (t,state) current = +  case action state current (nodeValue t current) of +     Stop newState -> (t,newState)+     ModifyAndStop _ newValue -> (setNodeValue current newValue t, state) +     Skip newState -> foldl' (_traverseTree False action) (t,newState) (nodeChildren t current)+     Modify newState newValue -> +         let newTree = setNodeValue current newValue t +         in +         foldl' (_traverseTree False action) (newTree,newState) (nodeChildren newTree current)++mapWithCluster :: Ord c +               => (c -> NodeValue f -> NodeValue f)+               -> JTree c f +               -> JTree c f        +mapWithCluster f t = t {nodeValueMap = Map.mapWithKey f (nodeValueMap t)}++-- | Set the factors in the tree +setFactors :: (Graph g, Factor f, Show f, IsCluster c, Ord c)+           => BayesianNetwork g f +           -> JTree c f +           -> JTree c f+setFactors g t = +  let factors = allVertexValues g +  in+  fst . traverseTree updateFactor factors $ t +++-- | Update factors in a cluster+updateFactor :: (Factor f, IsCluster c) +             => [f] -- ^ Remaining list of factors to attribute+             -> c -- ^ Current cluster+             -> NodeValue f -- ^ Current value+             -> Action [f] (NodeValue f)+updateFactor lf c (NodeValue _ evidence) | null lf = Stop lf+                                         | otherwise =+  let isFactorIncluded l = all (`elem` clusterVariables c) (factorVariables l)+      (attributedFactors,remainingFactors) = partition isFactorIncluded lf +  in +  Modify remainingFactors  (NodeValue (factorProduct attributedFactors) evidence)++-- | Change evidence in the network+changeEvidence :: (IsCluster c, Ord c, Factor f, Message f c)+               => [DVI Int] -- ^ Evidence+               -> JTree c f +               -> JTree c f +changeEvidence e t =  distribute . +                      collect . fst .+                      traverseTree changeNodeEvidence e $ +                      t { separatorValueMap = Map.map (const EmptySeparator) (separatorValueMap t)}++changeNodeEvidence :: (IsCluster c, Factor f) +                   => [DVI Int] -- ^ Evidence+                   -> c  -- ^ Current cluster+                   -> NodeValue f -- ^ Current value+                   -> Action [DVI Int] (NodeValue f)+changeNodeEvidence [] c v = Stop []+changeNodeEvidence e c (NodeValue f olde) = +  let oe = overlappingEvidence c e+      ns = e \\ oe+      newEvidence = factorProduct $ map factorFromInstantiation oe+  in +  Modify ns (NodeValue f newEvidence)+++-- | Cluster of discrete variables.+-- Discrete variables instead of vertices are needed because the+-- factor are using 'DV' and we need to find+-- which factors must be contained in a given cluster.+newtype Cluster = Cluster (Set.Set DV) deriving(Eq,Ord)++instance Show Cluster where +  show (Cluster s) = show . Set.toList $ s++fromCluster (Cluster s) = Set.toList s +++instance Factor f => Message f Cluster where +  newMessage input (NodeValue f e) dv = factorProjectTo (fromCluster dv) (factorProduct (f:e:input))+++type JunctionTree f = JTree Cluster f++{-++Implement the show function to see the structure of the tree+(without the values)+    +-}++data NodeKind c = N !c | S !c++label True c a = c ++ "=" ++ show a +label False c _ = c++-- | Convert the JTree into a tree of string+-- using the cluster.+toTree :: (Ord c, Show c, Show a) +       => Bool -- ^ True if the data must be displayed+       -> JTree c a +       -> Tree.Tree String+toTree d t = +    let r = root t+        v = nodeValue t r+        nodec = map S (nodeChildren t r)+    in +    Tree.Node (label d (show r) v) (_toTree d t nodec)+++_toTree :: (Ord c, Show c, Show a) +        => Bool+        -> JTree c a +        -> [NodeKind c] +        -> [Tree.Tree String]+_toTree _ _ [] = []+_toTree d t ((N h):l) = +    let nodec = map S (nodeChildren t h) -- Node children are separators+        v = nodeValue t h+    in+    Tree.Node (label d (show h) v) (_toTree d t nodec):_toTree d t l+_toTree d t ((S h):l) = +    let separatorc = [N $ separatorChild t h] -- separator child is a node+        v = separatorValue t h+    in+    Tree.Node (label d ("<" ++ show h ++ ">") v ) (_toTree d t separatorc):_toTree d t l++instance (Ord c, Show c, Show a) => Show (JTree c a) where +    show = Tree.drawTree . toTree False++displayTree b = Tree.drawTree . toTree b++{-++Debug functions for tests++-}++--instance Message (Sum Int) String where +--    newMessage l (NodeValue a b) _ = mconcat (a:b:l)+--+--testTree :: JTree String (Sum Int)+--testTree = let s a= Sum a+--           in+--           addSeparator "ROOT" "RB" "B" .+--           addNode "B" (s 3) (s 3) . +--           addSeparator "ROOT" "RA" "A"  . +--           addNode "A" (s 2) (s 2) $ +--           singletonTree "ROOT" (s 4) (s 5)+--

@@ -12,6 +12,7 @@ import qualified Data.Map as Map import Bayes.Factor import Bayes+import Bayes.PrivateTypes  --import Debug.Trace  @@ -139,7 +140,7 @@         conds = tail dvs         oldOrder = conds ++ [dst]         dvalues = map read values :: [Double]-        newvalues = changeVariableOrder oldOrder dvs dvalues+        newvalues = changeVariableOrder (DVSet oldOrder) (DVSet dvs) dvalues     cpt dst conds ~~ newvalues     return () 

@@ -0,0 +1,184 @@+{- | Private types for Bayes and Factors.++Those type are not exported++-}++module Bayes.PrivateTypes( + -- * Classes+   BayesianDiscreteVariable(..)+ , Set(..)+ -- * Variables+ , DV(..)+ , DVSet(..)+ -- * Instantiations+ , DVI(..)+ , setDVValue+ , (=:)+ , instantiationValue+ , instantiationVariable+ -- * Vertices + , Vertex(..)+ -- * Misc+ , getMinBound+ -- * Indices + , Index(..)+ , forAllInstantiations + , indicesForDomain+ , fromIndex+ , instantiationDetails+ , allInstantiationsForOneVariable+ ) where+++import qualified Data.List as L ++{-+	Set+-}++-- | A Set of variables used in a factor. s is the set and a the variable+class Set s where+    -- | Empty set+    emptySet :: s a+    -- | Union of two sets+    union :: Eq a => s a -> s a -> s a+    -- | Intersection of two sets+    intersection :: Eq a => s a -> s a -> s a+    -- | Difference of two sets+    difference :: Eq a => s a -> s a -> s a+    -- | Check if the set is empty+    isEmpty :: s a -> Bool+    -- | Check if an element is member of the set+    isElem :: Eq a => a -> s a -> Bool+    -- | Add an element to the set+    addElem :: Eq a => a -> s a -> s a+    -- | Number of elements in the set+    nbElements :: s a -> Int++    -- | Check if a set is subset of another one+    subset :: Eq a => s a -> s a -> Bool++    -- | Check set equality+    equal :: Eq a => s a -> s a -> Bool+    equal sa sb = (sa `subset` sb) && (sb `subset` sa)++instance Set [] where+    emptySet = []+    union = L.union+    intersection = L.intersect+    difference a b = a L.\\ b+    isEmpty [] = True +    isEmpty _ = False+    isElem = L.elem +    addElem a l = if a `elem` l then l else a:l+    nbElements = length+    subset sa sb = all (`elem` sb) sa++{-++Misc++-}+-- | Vertex type used to identify a vertex in a graph+newtype Vertex = Vertex {vertexId :: Int} deriving(Eq,Ord)++instance Show Vertex where +    show (Vertex v) = "v" ++ show v++-- | A discrete variable has a number of levels which is required to size the factors+class BayesianDiscreteVariable v where+    dimension :: v -> Int ++-- | Get the minimum bound for a type+getMinBound :: Bounded a => a -> a +getMinBound _ = minBound++++{-++Variables+	+-}++-- | A discrete variable+data DV = DV !Vertex !Int deriving(Eq,Ord)++-- | A set of discrete variables+-- The tag is used to check that an index is used with the right set of DV+newtype DVSet s = DVSet [DV] deriving(Eq)++-- | Remove the type tag when not needed+fromDVSet :: DVSet s -> [DV]+fromDVSet (DVSet l) = l++instance Show DV where+    show (DV v d) = show v ++ "(" ++ show d ++ ")"++instance BayesianDiscreteVariable DV where+    dimension (DV _ d) = d++{-++Index++-}++-- | An index with meaning only for a given DVSet+newtype Index s = Index Int deriving(Eq)++-- | Used to forget the type tag+fromIndex :: Index s -> Int +fromIndex (Index i) = i ++-- | Generate all the indices for a set of variables+indicesForDomain :: DVSet s -> [[Index s]]+indicesForDomain (DVSet l) = mapM indicesForOneDomain l+ where + 	indicesForOneDomain (DV _ d) = map Index [0..d-1]++allInstantiationsForOneVariable :: DV -> [DVI Int]+allInstantiationsForOneVariable v@(DV _ d) = map (setDVValue v) [0..d-1]++-- | Generate all instantiations of variables+-- The DVInt can be in any order so the tag s is not used+forAllInstantiations :: DVSet s -> [[DVI Int]]+forAllInstantiations (DVSet l) = mapM allInstantiationsForOneVariable l+ ++{- ++Instantiations++-}+-- | Discrete Variable instantiation. A variable and its value+data DVI a = DVI DV !a deriving(Eq)++instance Show a => Show (DVI a) where +   show (DVI (DV v _) i) = show v ++ "=" ++ show i++   -- | A set of variable instantiations+type DVISet a = [DVI a]++-- | Create a discrete variable instantiation for a given discrete variable+setDVValue :: DV -> a -> DVI a+setDVValue v a = DVI v a++-- | Create a variable instantiation using values from+-- an enumeration+(=:) :: (Bounded b, Enum b) => DV -> b -> DVI Int +(=:) a b = setDVValue a (fromEnum b - fromEnum (getMinBound b))++instance BayesianDiscreteVariable (DVI a) where+    dimension (DVI v _) = dimension v++-- | Get the variables and their values with a type constraint+instantiationDetails :: [DVI Int] -> (DVSet s, [Index s])+instantiationDetails l = (DVSet $ map instantiationVariable l, map (Index . instantiationValue) l)++-- | Extract value of the instantiation+instantiationValue (DVI _ v) = v++-- | Discrete variable from the instantiation+instantiationVariable (DVI dv _) = dv

@@ -13,6 +13,10 @@ import Bayes.Factor(testProductProject_prop,testScale_prop,testProjectCommut_prop,testScalarProduct_prop,testProjectionToScalar_prop) import Bayes.FactorElimination(junctionTreeProperty_prop) +#ifdef LOCAL+import Bayes.Test.ReferencePatterns(compareAsiaReference,compareCancerReference,comparePokerReference,compareFarmReference)+#endif + -- | Run all the tests runTests = defaultMain tests @@ -32,6 +36,14 @@                 testProperty "Test the junction tree property" junctionTreeProperty_prop,                 testCase "Test variable elimination == factor elimination" compareVariableFactor             ]+#ifdef LOCAL+        , testGroup "Reference patterns" [ +                testCase "Asia reference pattern" compareAsiaReference,+                testCase "Cancer reference pattern" compareCancerReference,+                testCase "Poker reference pattern" comparePokerReference,+                testCase "Farm reference pattern" compareFarmReference+        ]+#endif      ] 

@@ -32,8 +32,8 @@         jt = createJunctionTree nodeComparisonForTriangulation exampleG     compareFactors "PRIOR FOR RAIN" (posterior jt rain) (priorMarginal exampleG [winter,sprinkler,wet,road] [rain]) -    let jt1 = updateEvidence [wet =: True] jt -        jt2 = updateEvidence [wet =: True, sprinkler =: True] jt1 +    let jt1 = changeEvidence [wet =: True] jt +        jt2 = changeEvidence [wet =: True, sprinkler =: True] jt1       compareFactors "POSTERIOR RAIN FOR WET" (posterior jt1 rain)           (posteriorMarginal exampleG [winter,sprinkler,wet,road] [rain]  [wet =: True])

@@ -0,0 +1,122 @@+{- | A comparison of factor elimination with reference values generated with another bayesian network software++It is a non regression test. The test patterns are not provided with this package.+So, those tests are disabled by default in the hackage version.++-}+module Bayes.Test.ReferencePatterns(+#ifdef LOCAL+   compareAsiaReference+ , compareCancerReference+ , comparePokerReference+ , compareFarmReference+#endif+ ) where++import Test.HUnit.Base(assertBool)+import Data.Maybe(fromJust)+import qualified Data.Map as Map+import Bayes.Factor+import Bayes+import Bayes.FactorElimination+import Bayes.Examples(anyExample)+import Bayes.FactorElimination.JTree(root)++value varmap jt s = +  let v =  fromJust $ Map.lookup s varmap+    in +    factorToList (fromJust $ posterior jt v) ++testWithRef varmap jt s l = assertBool s $ value varmap jt s ~=~ l+testWithRefAndPrint varmap jt s l = do+  let r = value varmap jt s +  putStrLn $ "Computed:" ++ show r+  putStrLn $ "Reference:" ++ show l+  putStrLn ""+  assertBool s $ r ~=~ l++-- Check that the float values are equal with an accuracy < 0.01%+comparePercent :: Double -> Double -> Bool+comparePercent a b = abs (a-b) < 1e-4++(~=~) a b = and (zipWith comparePercent a b)++#ifdef LOCAL+compareFarmReference = do +  (varmap,g) <- anyExample "studfarm.net"+  let jt = createJunctionTree nodeComparisonForTriangulation g+  +  assertBool "Junction Tree property" $ junctionTreeProperty jt [] (root jt)+  testWithRef varmap jt "L"  [0.01,0.99]+  testWithRef varmap jt "Ann"  [0.01,0.99]+  testWithRef varmap jt "Brian"  [0.01,0.99]+  testWithRef varmap jt "Cecily"  [0.01,0.99]+  testWithRef varmap jt "K"  [0.01,0.99]+  testWithRef varmap jt "Fred"  [0.01,0.99]+  testWithRef varmap jt "Dorothy"  [0.01,0.99]+  testWithRef varmap jt "Eric"  [0.01,0.99]+  testWithRef varmap jt "Gwenn"  [0.01,0.99]+  testWithRef varmap jt "Henry"  [0.0091,0.9909]+  testWithRef varmap jt "Irene"  [0.0099,0.9901]+  testWithRef varmap jt "John"  [0.0004,0.0087,0.9909]+++comparePokerReference = do +  (varmap,g) <- anyExample "poker.net"+  let jt = createJunctionTree nodeComparisonForTriangulation g++  assertBool "Junction Tree property" $ junctionTreeProperty jt [] (root jt)+  testWithRef varmap jt "OH0"  [0.1672, 0.0445,0.0635,0.4659,0.1694,0.0494,0.0353,0.0024,0.0024]+  testWithRef varmap jt "OH1"  [0.0265,0.0170,0.0357,0.4125,0.2633,0.1599,0.0676,0.0098,0.0077]+  testWithRef varmap jt "OH2"  [0.2472,0.0628,0.2903,0.0258,0.2526,0.0881,0.0212,0.0121]+  testWithRef varmap jt "SC"  [0.2450,0.7116,0.0435]+  testWithRef varmap jt "FC"  [0.0895,0.6988,0.0445,0.1672]+  testWithRef varmap jt "Besthand"  [0.6396,0.3604]+  testWithRef varmap jt "MH"  [0.1250,0.1250,0.1250,0.1250,0.1250,0.1250,0.1250,0.1250]+++compareAsiaReference = do +  (varmap,g) <- anyExample "asia.net"+  let jt = createJunctionTree nodeComparisonForTriangulation g++  assertBool "Junction Tree property" $ junctionTreeProperty jt [] (root jt)+  testWithRef varmap jt "A"  [0.0100, 0.9900]+  testWithRef varmap jt "S"  [0.5000, 0.5000]+  testWithRef varmap jt "T"  [0.0104, 0.9896]+  testWithRef varmap jt "L"  [0.0550, 0.9450]+  testWithRef varmap jt "B"  [0.4500, 0.5500]+  testWithRef varmap jt "E"  [0.0648, 0.9352]+  testWithRef varmap jt "X"  [0.1103, 0.8897]+  testWithRef varmap jt "D"  [0.4360, 0.5640]++-- | Type defined to set the evidence on the Coma variable+-- from the cancer network.+data Coma = Present | Absent deriving(Eq,Enum,Bounded)++compareCancerReference = do +  (varmap,g) <- anyExample "cancer.net"+  let jt = createJunctionTree nodeComparisonForTriangulation g++  assertBool "Junction Tree property" $ junctionTreeProperty jt [] (root jt)+  testWithRef varmap jt "A"  [0.2000, 0.8000]+  testWithRef varmap jt "B"  [0.3200, 0.6800]+  testWithRef varmap jt "C"  [0.0800, 0.9200]+  testWithRef varmap jt "D"  [0.3200, 0.6800]+  testWithRef varmap jt "E"  [0.6160, 0.3840]++  let varD = fromJust $ Map.lookup "D" varmap+  let jt' = changeEvidence [varD =: Present] jt +  testWithRef varmap jt' "A"  [0.4250, 0.5750]+  testWithRef varmap jt' "B"  [0.8000, 0.2000]+  testWithRef varmap jt' "C"  [0.2000, 0.8000]+  testWithRef varmap jt' "D"  [1.0000, 0.0000]+  testWithRef varmap jt' "E"  [0.6400, 0.3600]++  let jt'' = changeEvidence [varD =: Absent] jt'+  testWithRef varmap jt'' "A"  [0.0941, 0.9059]+  testWithRef varmap jt'' "B"  [0.0941, 0.9059]+  testWithRef varmap jt'' "C"  [0.0235, 0.9765]+  testWithRef varmap jt'' "D"  [0.0000, 1.0000]+  testWithRef varmap jt'' "E"  [0.6047, 0.3953]++#endif

@@ -26,12 +26,12 @@ --debug s a = trace (s  ++ "\n" ++ show a ++ "\n") a  -- | Elimination order-type EliminationOrder = DVSet+type EliminationOrder = [DV]  -- | Get all variables from a Bayesian Network allVariables :: (Graph g, Factor f)               => BayesianNetwork g f -             -> DVSet+             -> [DV] allVariables g =    let s = allVertexValues g        createDV = factorMainVariable 

@@ -7,7 +7,7 @@ -- The package version. See the Haskell package versioning policy -- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for -- standards guiding when and how versions should be incremented.-Version:             0.1+Version:             0.2  -- A short (one-line) description of the package. Synopsis:            Inference with Discrete Bayesian Networks@@ -15,9 +15,7 @@ -- A longer description of the package. Description:  Algorithms for inference with Discrete Bayesian Networks.    It is a very preliminary version. It has only been tested on very simple- examples where it worked. On bigger networks, imported from Hugin files, it was very very very slow.- So, you can use this software as a toy. Much more work is needed to validate- and optimize it.     + examples where it worked. This 0.2 version is using new faster and cleaner algorithms.  -- URL for the project homepage or repository. Homepage:            http://www.alpheccar.org@@ -48,7 +46,16 @@  data-files: cancer.net +Flag local {+  Description: Enable local tests by the author. They can only be run with some test patterns not distributed with this package.+  Default:     False+} +Flag profile {+  Description: Build profiling version of the library too.+  Default:     False+}+ Library   -- Modules exported by the library.   Exposed-modules:@@ -61,11 +68,19 @@     Bayes.Test.CompareEliminations     Bayes.Examples     Bayes.Examples.Tutorial+    Bayes.Test.ReferencePatterns   other-modules:     Paths_hbayes     Bayes.ImportExport.HuginNet.Splitting+    Bayes.PrivateTypes+    Bayes.FactorElimination.JTree    GHC-Options: -O2 -funbox-strict-fields+  Extensions: CPP+  if flag(local)+    cpp-options: -DLOCAL+  if flag(profile)+    ghc-options: -auto-all       -- Packages needed in order to build this package.