hbayes (empty) → 0.1
raw patch · 14 files changed
+3331/−0 lines, 14 filesdep +HUnitdep +QuickCheckdep +arraysetup-changed
Dependencies added: HUnit, QuickCheck, array, base, boxes, containers, directory, filepath, mtl, parsec, pretty, random, split, test-framework, test-framework-hunit, test-framework-quickcheck2, vector
Files
- Bayes.hs +881/−0
- Bayes/Examples.hs +186/−0
- Bayes/Examples/Tutorial.hs +244/−0
- Bayes/Factor.hs +610/−0
- Bayes/FactorElimination.hs +654/−0
- Bayes/ImportExport/HuginNet.hs +192/−0
- Bayes/ImportExport/HuginNet/Splitting.hs +9/−0
- Bayes/Test.hs +38/−0
- Bayes/Test/CompareEliminations.hs +47/−0
- Bayes/VariableElimination.hs +219/−0
- LICENSE +30/−0
- Setup.hs +4/−0
- cancer.net +128/−0
- hbayes.cabal +89/−0
+ Bayes.hs view
@@ -0,0 +1,881 @@+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{- | Discrete Bayesian Network Library.++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. ++Look at the "Bayes.Examples" and "Bayes.Examples.Tutorial" in this package +to see how to use the library.++-}+module Bayes(+ -- * Graph+ -- ** Graph classes+ Graph(..)+ , UndirectedGraph(..)+ , DirectedGraph(..)+ , FoldableWithVertex(..)+ , NamedGraph(..)+ -- ** Graph Monad+ , GraphMonad+ , GMState(..)+ , graphNode+ , runGraph+ , execGraph+ , evalGraph+ -- ** Support functions for Graph constructions+ , Vertex+ , Edge + , edge+ , newEdge+ , edgeEndPoints+ , connectedGraph+ -- * SimpleGraph implementation+ -- ** The SimpleGraph type+ , DirectedSG+ , UndirectedSG+ -- ** Bayesian network+ , SBN+ , BayesianNetwork(..)+ -- * Bayesian Monad used to ease creation of Bayesian Networks+ , BNMonad+ , runBN + , evalBN+ , execBN+ , variable+ , variableWithSize+ , cpt+ , proba+ , t+ , (~~)+ -- * Testing+ , testEdgeRemoval_prop+ , testVertexRemoval_prop+) where++import qualified Data.IntMap as IM+import qualified Data.Map as M+import Control.Monad.State.Strict+import Control.Monad.Writer.Strict+import Control.Applicative((<$>))+import Bayes.Factor+import Data.Maybe+import qualified Data.Map as Map+import qualified Data.Foldable as F+import qualified Data.Traversable as T +import Control.Applicative +import qualified Data.Set as Set++import Test.QuickCheck+import Test.QuickCheck.Arbitrary+import Data.List(sort,intercalate,nub)++--import Debug.Trace+--debug a = trace (show a) a++-- | Bayesian network. g must be a directed graph and f a factor+type BayesianNetwork g f = g () f++instance Arbitrary (DirectedSG String String) where+ arbitrary = do + let createVertex g i = do + name <- arbitrary :: Gen String+ return $ addVertex (Vertex i) name g+ createEdge g (va,vb) = do + name <- arbitrary :: Gen String+ return $ addEdge (edge va vb) name 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++instance Arbitrary (DirectedSG () String) where+ arbitrary = do + let createVertex g i = do + name <- arbitrary :: Gen String+ return $ addVertex (Vertex i) name 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+ Just (vs,ve) = edgeVertices g e+ Just bi = ingoing g ve+ Just bo = outgoing g vs+ g' = removeEdge e g + Just bi' = ingoing g' ve+ Just bo' = outgoing g' vs+ in+ (map (sort . (:) e ) [bi', bo'] == map sort [bi,bo]) &&+ (sort (allEdges g) == sort (e:allEdges g'))++testVertexRemoval_prop :: DirectedSG String String -> Property+testVertexRemoval_prop g = (not . hasNoVertices) g ==>+ let Just v = someVertex g+ Just bi = ingoing g v + Just bo = outgoing g v+ g' = removeVertex v g+ srcVertices = mapMaybe (startVertex g') bi+ dstVertices = mapMaybe (endVertex g') bo + isNotDstVertex = not . (v `elem`) . mapMaybe (endVertex g') . fromJust . outgoing g'+ isNotStartVertex = not . (v `elem`) . mapMaybe (startVertex g') . fromJust . ingoing g'+ in + (sort (allVertices g) == sort (v:allVertices g')) &&+ (all isNotDstVertex srcVertices) && (all isNotStartVertex dstVertices)+++-- | Graph class used for graph processing algorithms.+-- A graph processing algorithm does not have to know how the graph is implemented nor if it is+-- directed or undirected+class Graph g where+ -- | Add a new vertex+ addVertex :: Vertex -> b -> g a b -> g a b+ -- | Remove a vertex+ removeVertex :: Vertex -> g a b -> g a b+ -- | Get the vertex value if the vertex is found in the graph+ vertexValue :: g a b -> Vertex -> Maybe b+ -- | Change the vertex value if the vertex is found in the graph+ changeVertexValue :: Vertex -> b -> g a b -> Maybe (g a b)+ -- | Generate a \"random\" vertex+ someVertex :: g a b -> Maybe Vertex++ -- | Check is the graph has no vertrex+ hasNoVertices :: g a b -> Bool++ -- | Generate all vertices+ allVertices :: g a b -> [Vertex]++ -- | Get all the values+ allVertexValues :: g a b -> [b]++ -- | Get all nodes+ allNodes :: g a b -> [(Vertex,b)]++ -- | Check if two vertices are linked by a vertex+ isLinkedWithAnEdge :: g a b -> Vertex -> Vertex -> Bool++ -- | Add an edge+ addEdge :: Edge -> a -> g a b -> g a b++ -- | Remove an dedge+ removeEdge :: Edge -> g a b -> g a b++ -- | Vertices for an edge+ edgeVertices :: g a b -> Edge -> Maybe (Vertex,Vertex)++ -- | Edge value if the edge is found in the graph+ edgeValue :: g a b -> Edge -> Maybe a++ -- | Return a \"random\" edge+ someEdge :: g a b -> Maybe Edge++ -- | Check if the graph has no edges+ hasNoEdges :: g a b -> Bool++ -- | One extremity of the edge (which is the end only for directed edge)+ endVertex :: g a b -> Edge -> Maybe Vertex+ endVertex g e = do + (_,ve) <- edgeVertices g e+ return ve + + -- | One extremity of the edge (which is the start only for directed edge)+ startVertex :: g a b -> Edge -> Maybe Vertex+ startVertex g e = do + (vs,_) <- edgeVertices g e+ return vs++ -- | All edges of the graph+ allEdges :: g a b -> [Edge]++ -- | All values of the graph+ allEdgeValues :: g a b -> [a]+ + -- | Returns an empty graph+ emptyGraph :: g a b++ -- | Check if the graph is empty+ isEmpty :: g a b -> Bool+ isEmpty g = hasNoVertices g && hasNoEdges g++ -- | Check if the graph is oriented+ oriented :: g a b -> Bool++ -- | All the neighbors of a vertex+ neighbors :: g a b -> Vertex -> Maybe [Vertex]++-- | A named graph is a graph where the vertices have a name.+-- This name is not a vertex value. Putting this name in the vertex value+-- would make algorithm less readable.+-- A vertex name is only useful to display the graph.+-- Labeled graph has a different meaning in graph theory.+class Graph g => NamedGraph g where+ -- | Add a vertex with a vertex name in addition to the value+ addLabeledVertex :: String -> Vertex -> b -> g a b -> g a b+ -- | Returns the vertex label+ vertexLabel :: g a b -> Vertex -> Maybe String+++-- | Undirected graph+class Graph g => UndirectedGraph g where+ edges :: g a b -> Vertex -> Maybe [Edge]++-- | Directed graph+class Graph g => DirectedGraph g where+ ingoing :: g a b -> Vertex -> Maybe [Edge]+ outgoing :: g a b -> Vertex -> Maybe [Edge]+++-- | Check if the graph is connected+connectedGraph :: Graph g => g a b -> Bool +connectedGraph g = + let visited = visitVertex g (Set.empty) ([fromJust $ someVertex g])+ vertices = Set.fromList $ allVertices g+ equalSets a b = Set.isSubsetOf a b && Set.isSubsetOf b a+ in + equalSets visited vertices+ where + visitVertex _ visited [] = visited+ visitVertex theGraph visited (current:n) = + if Set.member current visited+ then + visitVertex theGraph visited n+ else+ let n' = fromJust $ neighbors theGraph current+ in+ visitVertex theGraph (Set.insert current visited) (n ++ n')++++ +++-- | Edge type used to identify and edge in a graph+data Edge = Edge !Vertex !Vertex deriving(Eq,Ord,Show)++-- | Create an edge description+edge :: Vertex -> Vertex -> Edge +edge a b = Edge a b++-- | Endpoints of an edge+edgeEndPoints :: Edge -> (Vertex,Vertex)+edgeEndPoints (Edge va vb) = (va,vb)+++-- | Synonym for undefined because it is clearer to use t to set the Enum bounds of a variable+t = undefined++-- | Neighborhood structure for directed or undirected edges+-- | Directed edges+data DE = DE ![Edge] ![Edge] deriving(Eq,Show)++-- | Undirected edges+data UE = UE ![Edge] deriving(Eq,Show)++-- | Class used to share as much code as possible between+-- directed and undirected graphs without+-- implementing an undirected graph as a graph where+-- we have a directed edge in both directions +class NeighborhoodStructure n where+ -- | Return an empty neighborhood+ emptyNeighborhood :: n + -- | Ingoing edges+ ingoingNeighbors :: n -> [Edge]+ -- | Outgoing edge+ outgoingNeighbors :: n -> [Edge]+ -- | Remove an edge+ removeNeighborsEdge :: Edge -> n -> n+ -- | Add an outgoing edge+ addOutgoingEdge :: Edge -> n -> n+ -- Add in ingoing edge+ addIngoingEdge :: Edge -> n -> n++-- | Directed neighborhood structure for a vertex+instance NeighborhoodStructure DE where+ emptyNeighborhood = DE [] []+ ingoingNeighbors (DE i _) = i+ outgoingNeighbors (DE _ o) = o + removeNeighborsEdge e (DE i o) = + let i' = filter (/= e) i+ o' = filter (/= e) o + in + DE i' o'+ addOutgoingEdge e (DE i o) = DE i (e:o)+ addIngoingEdge e (DE i o) = DE (e:i) o++-- | Undirected neighborhood structure for a vertex+instance NeighborhoodStructure UE where+ emptyNeighborhood = UE []+ ingoingNeighbors (UE e) = e+ outgoingNeighbors (UE e) = e+ removeNeighborsEdge e (UE l) = + let l' = filter (/= e) l+ in + UE l'+ addOutgoingEdge e (UE l) = UE (e:l)+ addIngoingEdge e (UE l) = UE (e:l)++-- | Implementtaion of a SimpleGraph+data SimpleGraph local edgedata vertexdata = SP {+ -- | Mapping of edge to edge data+ edgeMap :: !(M.Map Edge edgedata) + -- ^ Mapping of vertex number to vertex neighborhood and vertex data+ , vertexMap :: !(IM.IntMap (local, vertexdata))+ -- ^ Vertex names. Used only to generate the graphviz representtaion. Names are useless for the algorithms+ -- and I don't want them to appear in the vetex values which should only be factor. Otherwise, the algorithms+ -- are less elegant since I have to extract the factors from the values+ , nameMap :: !(IM.IntMap String)+ } ++-- | Directed simple graph+type DirectedSG = SimpleGraph DE++-- | Undirected simple graph+type UndirectedSG = SimpleGraph UE++instance (Eq a, Eq b) => Eq (SimpleGraph DE a b) where+ (==) (SP a b _) (SP a' b' _) = a == a' && b == b'++-- | An empty simple graph+emptySimpleGraph = SP M.empty IM.empty IM.empty++-- | Used to prevent adding duplicates to a graph+noRedundancy new old = old++instance Functor (SimpleGraph local edge) where + fmap f (SP em vm nm) = SP em (IM.map (\(l,d) -> (l, f d)) vm) nm++instance F.Foldable (SimpleGraph local edge) where+ foldr f c (SP _ vm _) = IM.foldr (\(_,d) s -> f d s) c vm++instance T.Traversable (SimpleGraph local edge) where+ traverse f (SP em vm nm) = + let l = IM.toList vm -- [(IM.Key, (DE, String))]+ onTriple f (k,(l,v)) = (\z -> (k,(l,z))) <$> f v+ l' = T.traverse (onTriple f) l -- f [(k,(l,z))]+ result y = (\x -> SP em (IM.fromList x) nm) <$> y+ in + result l'++-- | The foldable class is limited. For a graph g we may need the vertex in addition to the value+class FoldableWithVertex g where+ -- | Fold with vertex + foldrWithVertex :: (Vertex -> a -> b -> 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++_addLabeledVertex vertexName vert@(Vertex v) value (SP em vm name) =+ let vm' = IM.insertWith noRedundancy v (emptyNeighborhood,value) vm+ name' = IM.insert v vertexName name + in+ SP em vm' name'++_vertexLabel (SP _ _ name) (Vertex v) = IM.lookup v name++instance NamedGraph DirectedSG where+ addLabeledVertex = _addLabeledVertex+ vertexLabel = _vertexLabel++instance NamedGraph UndirectedSG where+ addLabeledVertex = _addLabeledVertex+ vertexLabel = _vertexLabel++-- | SimpleGraph is an instance of Graph.+instance Graph DirectedSG where+ addVertex = _addVertex+ removeVertex = _removeVertex+ vertexValue = _vertexValue+ changeVertexValue = _changeVertexValue+ someVertex = _someVertex+ hasNoVertices = _hasNoVertices+ allVertices = _allVertices+ allVertexValues = _allVertexValues+ allNodes = _allNodes+ isLinkedWithAnEdge = _isLinkedWithAnEdge+ addEdge = _addEdge+ removeEdge = _removeEdge+ edgeVertices = _edgeVertices+ edgeValue = _edgeValue+ someEdge = _someEdge+ hasNoEdges = _hasNoEdges+ allEdges = _allEdges+ allEdgeValues = _allEdgeValues+ emptyGraph = _emptyGraph+ oriented _ = True+ neighbors g v = nub <$> liftA2 (++) + (map (\(Edge _ e) -> e) <$> (outgoing g v)) + (map (\(Edge s _) -> s) <$> (ingoing g v))++-- | Reverse edge direction+reverseEdge :: Edge -> Edge +reverseEdge (Edge va vb) = edge vb va++-- | SimpleGraph is an instance of Graph.+instance Graph UndirectedSG where+ addVertex = _addVertex+ removeVertex = _removeVertex+ vertexValue = _vertexValue+ changeVertexValue = _changeVertexValue+ someVertex = _someVertex+ hasNoVertices = _hasNoVertices+ allVertices = _allVertices+ allVertexValues = _allVertexValues+ allNodes = _allNodes+ isLinkedWithAnEdge = _isLinkedWithAnEdge+ addEdge = _addEdge+ removeEdge e g = _removeEdge (reverseEdge e) (_removeEdge e g)+ edgeVertices = _edgeVertices+ edgeValue g e = case _edgeValue g e of + Nothing -> _edgeValue g (reverseEdge e) + r@(Just _) -> r+ someEdge = _someEdge+ hasNoEdges = _hasNoEdges+ allEdges = _allEdges+ allEdgeValues = _allEdgeValues+ emptyGraph = _emptyGraph+ oriented _ = False+ -- in undirected graphs the edge direction does not count so we need to get both+ -- ends to be sure we don not forget a vertex. In addition to that, an end may be the current vertex which+ -- is not part of the neighbors. So it has to be filtered out. Obviously, a better solution will+ -- have to be designed.+ neighbors g v = filter (/= v) <$> nub <$> liftA2 (++) + (map (\(Edge _ e) -> e) <$> (edges g v)) + (map (\(Edge s _) -> s) <$> (edges g v))++_emptyGraph = emptySimpleGraph++_hasNoVertices (SP _ vm _) = IM.null vm++_hasNoEdges (SP em _ _) = M.null em++_allVertices (SP _ vm _) = map Vertex . IM.keys $ vm++_allEdges (SP em _ _) = M.keys $ em++_allNodes (SP _ vm _) = map (\(k,(_,v)) -> (Vertex k,v)) . IM.assocs $ vm++_allVertexValues (SP _ vm _) = map snd (IM.elems vm)++_allEdgeValues (SP em _ _) = M.elems em++_isLinkedWithAnEdge (SP em _ _) va vb = M.member (edge va vb) em || M.member (edge vb va) em++_someVertex (SP _ vm _) = + if IM.null vm + then + Nothing + else + Just . Vertex . head . IM.keys $ vm++_someEdge (SP em _ _) = + if M.null em + then + Nothing + 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++_removeVertex v@(Vertex vertex) g@(SP _ vm _) = maybe g removeVertexWithValue (IM.lookup vertex vm)+ where+ removeVertexWithValue (n,_) = let g' = foldr _removeEdge g (ingoingNeighbors n)+ SP em vm' nm' = foldr _removeEdge g' (outgoingNeighbors n)+ in + SP em (IM.delete vertex vm') nm'+_vertexValue g@(SP _ vm _) (Vertex i) = maybe Nothing (Just . extractValue) (IM.lookup i vm)+ where+ extractValue (_,d) = d++_changeVertexValue v@(Vertex vi) newValue g@(SP e vm nm) = + let newVertexMap = do+ (n,_) <- IM.lookup vi vm+ return $ IM.insert vi (n,newValue) vm+ in + case newVertexMap of + Nothing -> Just g+ Just nvm -> Just $ SP e nvm nm++_removeEdge e@(Edge (Vertex vs) (Vertex ve)) g@(SP em vm nm) = + let r = do + _ <- M.lookup e em -- Check e is member of the graph+ (ns,vsdata) <- IM.lookup vs vm+ (ne,vedata) <- IM.lookup ve vm+ return ((vs,(removeNeighborsEdge e ns,vsdata)),(ve,(removeNeighborsEdge e ne,vedata)))+ updateGraph ((vs,vsdata),(ve,vedata)) =+ let vm' = IM.insert ve vedata . IM.insert vs vsdata $ vm+ em' = M.delete e em + in + SP em' vm' nm+ in + maybe g updateGraph r++_edgeVertices (SP em _ _) e@(Edge vs ve) =+ if M.member e em + then + Just (vs,ve)+ else+ Nothing++_edgeValue (SP em _ _) e = do+ v <- M.lookup e em+ return v++_addEdge newEdge@(Edge vs ve) value g@(SP em vm nm) = + if testEdgeExistence g em vs ve + then + g+ else+ SP (M.insert newEdge value em) (addEdgeReference vm vs ve) nm+ where+ testEdgeExistence g em va vb = + if (oriented g)+ then + M.member (Edge va vb) em+ else + M.member (Edge va vb) em || M.member (Edge vb va) em + addEdgeReference vm (Vertex vsi) (Vertex vei) = IM.adjust addi vei (IM.adjust addo vsi vm)+ addi (n,v) = (addIngoingEdge newEdge n,v) + addo (n,v) = (addOutgoingEdge newEdge n,v) ++instance UndirectedGraph UndirectedSG where+ edges g@(SP _ vm _) v@(Vertex vi) =+ do + (n,_) <- IM.lookup vi vm+ return (ingoingNeighbors n)++instance DirectedGraph DirectedSG where+ ingoing g@(SP _ vm _) v@(Vertex vi) =+ do + (n,_) <- IM.lookup vi vm+ return (ingoingNeighbors n)++ outgoing g@(SP _ vm _) v@(Vertex vi) =+ do + (n,_) <- IM.lookup vi vm+ return (outgoingNeighbors n) ++{-+ +Following code is used to display a graph in a form adapted to humans.++-}+printNode nm (Vertex k,v) = do + tell "\n"+ let r = IM.lookup k nm+ when (isJust r) $ do+ tell $ fromJust r+ tell "\n"+ tell $ show v+ tell "\n"+addVertexToGraphviz nm (k,(_,v)) = do+ tell $ show k+ let r = IM.lookup k $ nm + when (isJust r) $ do+ tell " [label=\""+ tell $ fromJust r+ tell "\"] ;" + tell "\n"++instance (Show b, Show e) => Show (DirectedSG e b)where+ show g@(SP em vm nm) = execWriter $ do+ tell "digraph dot {\n"+ mapM_ (addVertexToGraphviz nm) $ IM.toList vm+ tell "\n"+ mapM_ addEdgeToGraphviz $ M.toList em+ tell "}\n"+ mapM_ (printNode nm) (allNodes g)+ where+ addEdgeToGraphviz (Edge (Vertex vs) (Vertex ve),l) = do+ tell $ show vs + tell " -> "+ tell $ show ve+ tell " [label=\""+ tell $ show l+ tell "\"]"+ tell ";\n"++instance (Show b, Show e) => Show (UndirectedSG e b)where+ show g@(SP em vm nm) = execWriter $ do+ tell "graph dot {\n"+ mapM_ (addVertexToGraphviz nm) $ IM.toList vm+ tell "\n"+ mapM_ addEdgeToGraphviz $ M.toList em+ tell "}\n"+ mapM_ (printNode nm) (allNodes g)+ where+ addEdgeToGraphviz (Edge (Vertex vs) (Vertex ve),l) = do+ tell $ show vs + tell " -- "+ tell $ show ve+ tell " [label=\""+ tell $ show l+ tell "\"]"+ tell ";\n"+++-- | Bayesian variable : name,dimension, factor+-- When initialized it is using a factor with bayesian variables.+-- But the factor value are not yet set+data MaybeBNode f = UninitializedBNode String Int+ | InitializedBNode String Int f+++displayFactors :: (NeighborhoodStructure n, Show f, Factor f, Graph (SimpleGraph n)) => SimpleGraph n a f -> String+displayFactors g@(SP _ _ nm) = + let nodes = allNodes g+ displayFactor (Vertex i,f) = + let s = fromJust . IM.lookup i $ nm+ in+ s ++ "\n" ++ show f+ in+ intercalate "\n" $ map displayFactor nodes++-- | An implementation of the BayesianNetwork using the simple graph and no value of edges+type SBN f = DirectedSG () f++-- | State used for the construction of the graph in the monad and containing+-- auxiliary informations like vertex name to vertex id and vertex count+type AuxiliaryState = (M.Map String Int, Int)++emptyAuxiliaryState = (M.empty,0)++-- | The State for the monad with a mapping from variable name to variable ID.+type BNState g f = GMState g () (MaybeBNode f)++-- | The Bayesian monad+type BNMonad g f a = GraphMonad g () (MaybeBNode f) a++-- | The state of the graph monad : the graph and auxiliary data+-- useful during the construction+type GMState g e f = (AuxiliaryState,g e f)++-- | Graph monad.+-- The monad used to simplify the description of a new graph+-- g is the graph type. e the edge type. f the node type (generally a 'Factor')+newtype GraphMonad g e f a = GM {runGraphMonad :: State (GMState g e f) a} deriving(Monad, MonadState (GMState g e f))++-- | Get the Bayesian Discrete Variable for a vertex.+-- It works because we keep the variable dimension+factorVariable :: Graph g => Vertex -> BNMonad g f (Maybe DV) +factorVariable v = do + g <- gets snd + let value = vertexValue g v+ case value of+ Nothing -> return Nothing+ Just (UninitializedBNode _ d) -> return $ Just $ DV v d+ Just (InitializedBNode _ d _) -> return $ Just $ DV v d+ ++-- | Get a named vertex from the graph monad+getVertex :: Graph g => String -> GraphMonad g e f (Maybe Vertex)+getVertex a = do+ (namemap,_) <- gets fst+ return $ do+ i <- M.lookup a namemap+ return (Vertex i)++-- | Create an edge between two vertex of the Bayesian network+(<--) :: Graph g => DV -> DV -> BNMonad g f ()+DV va _ <-- DV vb _ = newEdge vb va ()++-- | Add a new labeled edge to the graph+newEdge :: Graph g => Vertex -> Vertex -> e -> GraphMonad g e f ()+newEdge va vb e = do+ (aux,g) <- get + let g1 = addEdge (edge va vb) e g+ put $! (aux,g1)+ return ()++whenJust Nothing _ = return ()+whenJust (Just i) f = f i >> return ()++-- | Get the node of a bayesian network under creation+getBayesianNode :: Graph g => Vertex -> BNMonad g f (Maybe (MaybeBNode f))+getBayesianNode v = do+ g <- gets snd+ return $ vertexValue g v++-- | Set the node of a bayesian network under creation+setBayesianNode :: Graph g => Vertex -> MaybeBNode f -> BNMonad g f ()+setBayesianNode v newValue = do+ (aux,oldGraph) <- get+ let newGraph = changeVertexValue v newValue oldGraph+ + whenJust newGraph $ \nvm -> do+ put $! (aux, nvm)++-- | Initialize the values of a factor+(~~) :: (DirectedGraph g, Factor f) + => BNMonad g f DV -- ^ Discrete variable in the graph+ -> [Double] -- ^ List of values+ -> BNMonad g f ()+(~~) mv l = do + (DV v _) <- mv -- This is updating the state and so the graph+ g <- gets snd+ current <- factorVariable v+ mvalue <- getBayesianNode v+ maybe (return ()) (setCpt g v current) mvalue+ where+ setCpt g _ _ (InitializedBNode _ _ _) = return ()+ setCpt g v current (UninitializedBNode s dim) = do + let vertices = map (fromJust . startVertex g) . fromJust . ingoing g $ v+ fv <- mapM factorVariable vertices+ let cpt = factorWithVariables (map fromJust (current:fv)) l+ newValue r = InitializedBNode s dim r+ maybe (return ()) (setBayesianNode v . newValue) cpt++ +minBoundForEnum :: Bounded a => a -> a+minBoundForEnum _ = minBound++maxBoundForEnum :: Bounded a => a -> a+maxBoundForEnum _ = maxBound++intValue :: Enum a => a -> Int+intValue = fromEnum+++-- | Set the bound of a bayesian variable (number of levels)+setVariableBoundWithSize :: Graph g+ => Vertex -- ^ Vertex+ -> Int -- ^ Inf limit (0 for instance)+ -> Int -- ^ Sup limit (1 for instance for 2 elements)+ -> BNMonad g f ()+setVariableBoundWithSize a bmin bmax = do+ v <- getBayesianNode a+ whenJust v $ \(UninitializedBNode s _) -> do+ setBayesianNode a (UninitializedBNode s (bmax - bmin + 1))++setVariableBound :: (Enum a, Bounded a, Graph g) + => Vertex -- ^ Vertex+ -> a -- ^ Bounded variable (t :: type where t is undefined)+ -> BNMonad g f ()+setVariableBound a e = + let bmin = intValue $ minBoundForEnum e+ bmax = intValue $ maxBoundForEnum e+ in + setVariableBoundWithSize a bmin bmax++-- | Create a new named Bayesian variable if not found.+-- Otherwise, return the found one.+addVariableIfNotFound :: NamedGraph g => String -> BNMonad g f Vertex+addVariableIfNotFound vertexName = graphNode vertexName (UninitializedBNode vertexName 0)++-- | Add a node in the graph using the graph monad+graphNode :: NamedGraph g => String -> f -> GraphMonad g e f Vertex +graphNode vertexName initValue = do+ (aux@(namemap,_),g) <- get+ maybe (createAndReturnVertex aux g) returnVertex (M.lookup vertexName namemap)+ where+ returnVertex i = return (Vertex i)+ createAndReturnVertex (namemap,count) g = do+ let g1 = addLabeledVertex vertexName (Vertex count) initValue g+ namemap1 = M.insert vertexName count namemap+ put $! ((namemap1,count+1),g1)+ return (Vertex count)++-- | Define a Bayesian variable (name and bounds)+variable :: (Enum a, Bounded a, NamedGraph g) + => String -- ^ Variable name+ -> a -- ^ Variable bounds+ -> BNMonad g f DV+variable name e = do+ va <- addVariableIfNotFound name+ setVariableBound va e+ maybeValue <- getBayesianNode va + setBayesianNode va (fromJust maybeValue)+ case fromJust maybeValue of + UninitializedBNode s d -> return (DV va d)+ InitializedBNode _ d _ -> return (DV va d)++-- | Define a Bayesian variable (name and bounds)+variableWithSize :: NamedGraph g+ => String -- ^ Variable name+ -> Int -- ^ Variable size+ -> BNMonad g f DV+variableWithSize name e = do+ va <- addVariableIfNotFound name+ setVariableBoundWithSize va 0 (e-1)+ maybeValue <- getBayesianNode va + setBayesianNode va (fromJust maybeValue)+ case fromJust maybeValue of + UninitializedBNode s d -> return (DV va d)+ InitializedBNode _ d _ -> return (DV va d)++-- | Define a conditional probability between different variables+-- Variables are ordered like+-- FFF FFT FTF FTT TFF TFT TTF TTT+-- and same for other enumeration keeping enumeration order+cpt :: DirectedGraph g => DV -> [DV] -> BNMonad g f DV+cpt node conditions = do+ mapM_ (node <--) (reverse conditions)+ return node++-- | Define proba for a variable+-- Values are ordered like+-- FFF FFT FTF FTT TFF TFT TTF TTT+-- and same for other enumeration keeping enumeration order+proba :: DirectedGraph g => DV -> BNMonad g f DV+proba node = cpt node []+++runGraph :: Graph g => GraphMonad g e f a -> (a,g e f)+runGraph = removeAuxiliaryState . flip runState (emptyAuxiliaryState,emptyGraph) . runGraphMonad + where + removeAuxiliaryState (r,(_,g)) = (r,g)++evalGraph :: Graph g => GraphMonad g e f a -> a+evalGraph = flip evalState (emptyAuxiliaryState,emptyGraph) . runGraphMonad ++execGraph :: Graph g => GraphMonad g e f a -> g e f+execGraph = snd . flip execState (emptyAuxiliaryState,emptyGraph) . runGraphMonad ++-- | Create a bayesian network using the simple graph implementation+-- The initialized nodes are replaced by the factor.+-- Returns the monad values and the built graph.+runBN :: BNMonad DirectedSG f a -> (a,DirectedSG () f)+runBN x = + let (r,g) = runGraph x+ convertBNodes (InitializedBNode s d f) = f + convertBNodes (UninitializedBNode s d) = error $ "All variables must be initialized with a factor: " ++ s ++ "(" ++ show d ++ ")"+ in + (r,fmap convertBNodes g)++-- | Create a bayesian network but only returns the monad value.+-- Mainly used for testing.+evalBN :: BNMonad DirectedSG f a -> a+evalBN = evalGraph++-- | Create a bayesian network but only returns the monad value.+-- Mainly used for testing.+execBN :: BNMonad DirectedSG f a -> DirectedSG () f+execBN x = + let g = execGraph x+ convertBNodes (InitializedBNode s d f) = f + convertBNodes (UninitializedBNode s d) = error $ "All variables must be initialized with a factor: " ++ s ++ "(" ++ show d ++ ")"+ in + fmap convertBNodes g
+ Bayes/Examples.hs view
@@ -0,0 +1,186 @@+{- | Examples of networks++/Creating a simple network/++The 'example' function is the typical example.+It is using the monad 'BNMonad'. The goal of this monad is to offer+a way of describing the network which is natural.++There are only three functions to understand inside the monad:++ * 'variable' to create a discrete variable of type 'DV'. Creating a discrete+ variable is using a 'Bounded' and 'Enum' type like for instance 'Bool'.++ * 'proba' to define the probability P(A) of a variable A++ * 'cpt' to define the conditional probability table P(A | BC)++It is important to understand how the values are organized. If you define+P( wet | sprinkler road) then you have to give the values in the order:++@+wet=False, sprinkler=False, road=False+wet=False, sprinkler=False, road=True+wet=False, sprinkler=True, road=False+wet=False, sprinkler=True, road=True+@++Finally, don't forget to return the discrete variables at the end of your network+construction because those variables are used for making inferences.++@+example :: ('DVSet','SBN' 'CPT')+example = 'runBN' $ do + winter <- 'variable' \"winter\" (t :: Bool)+ sprinkler <- 'variable' \"sprinkler\" (t :: Bool) + wet <- 'variable' \"wet grass\" (t :: Bool) + rain <- 'variable' \"rain\" (t :: Bool) + road <- 'variable' \"slippery road\" (t :: Bool) +--+ 'proba' winter ~~ [0.4,0.6]+ 'cpt' sprinkler [winter] ~~ [0.25,0.8,0.75,0.2]+ 'cpt' rain [winter] ~~ [0.9,0.2,0.1,0.8]+ 'cpt' wet [sprinkler,rain] ~~ [1,0.2,0.1,0.05,0,0.8,0.9,0.95]+ 'cpt' road [rain] ~~ [1,0.3,0,0.7]+ return [winter,sprinkler,rain,wet,road]+@++/Importing a network from a Hugin file/++The 'exampleImport' function can be used to import a file in Hugin format.+Only a subset of the format is supported.+The function will return a mapping from node names to Discrete Variables 'DV'.+The node name is used and not the node's label.+The function is also returning a simple bayesian network 'SBN' using 'CPT'+as factors.++The implementation is using 'getDataFileName' to find the path of the+test pattern installed by cabal.++@+exampleImport :: IO (Map.Map String 'DV','SBN' 'CPT')+exampleImport = do + path <- 'getDataFileName' \"cancer.net\"+ r <- 'importBayesianGraph' path+ return ('runBN' $ fromJust r)+@++-}+module Bayes.Examples(+ example+ , exampleJunction+ , exampleImport+ , exampleDiabete+ , exampleAsia+ , examplePoker+ , exampleFarm+ , examplePerso+ , testJunction+ , anyExample+ ) where ++import Bayes+import Bayes.Factor+import Bayes.ImportExport.HuginNet+import Data.Maybe(fromJust)+import qualified Data.Map as Map+import System.Directory(getHomeDirectory)+import System.FilePath((</>))+import Paths_hbayes++-- | Example showing how to import a graph described into+-- a Hugin file.+exampleImport :: IO (Map.Map String DV,SBN CPT)+exampleImport = do + path <- getDataFileName "cancer.net"+ r <- importBayesianGraph path+ return (runBN $ fromJust r)++-- | Genereic loading functions to load some other+-- examples from the author's dropbox.+-- Those additional examples are not distributed with this package.+-- They are used only for testing and debugging purposes+genericExample :: String -> IO (Map.Map String DV,SBN CPT)+genericExample s = do + r <- importBayesianGraph s+ return (runBN $ fromJust r)++anyExample s = do+ h <- getHomeDirectory+ genericExample $ h </> "Dropbox/bayes_examples" </> s+ +-- | Diabete example (not provided with this package)+exampleDiabete = do + h <- getHomeDirectory+ genericExample $ h </> "Dropbox/bayes_examples/Diabetes.hugin"++-- | Asia example (not provided with this package)+exampleAsia = do + h <- getHomeDirectory+ genericExample $ h </> "Dropbox/bayes_examples/asia.net"++-- | Poker example (not provided with this package)+examplePoker = do + h <- getHomeDirectory+ genericExample $ h </> "Dropbox/bayes_examples/poker.net"++-- | Farm example (not provided with this package)+exampleFarm = do + h <- getHomeDirectory+ genericExample $ h </> "Dropbox/bayes_examples/studfarm.net"++-- | Perso example (not provided with this package)+examplePerso = do + h <- getHomeDirectory+ genericExample $ h </> "Dropbox/bayes_examples/mytest.net"+++-- | Standard example found in many books about Bayesian Networks.+example :: (DVSet,SBN CPT)+example = runBN $ do + winter <- variable "winter" (t :: Bool)+ sprinkler <- variable "sprinkler" (t :: Bool) + wet <- variable "wet grass" (t :: Bool) + rain <- variable "rain" (t :: Bool) + road <- variable "slippery road" (t :: Bool) ++ proba winter ~~ [0.4,0.6]+ cpt sprinkler [winter] ~~ [0.25,0.8,0.75,0.2]+ cpt rain [winter] ~~ [0.9,0.2,0.1,0.8]+ cpt wet [sprinkler,rain] ~~ [1,0.2,0.1,0.05,0,0.8,0.9,0.95]+ cpt road [rain] ~~ [1,0.3,0,0.7]+ return [winter,sprinkler,rain,wet,road]++testJunction :: DirectedSG () Vertex+testJunction = execGraph $ do+ a <- graphNode "A" (Vertex 0) + b <- graphNode "B" (Vertex 1) + c <- graphNode "C" (Vertex 2) + newEdge a b () + newEdge a c ()++exampleJunction :: UndirectedSG () Vertex+exampleJunction = execGraph $ do + a <- graphNode "A" (Vertex 0) + b <- graphNode "B" (Vertex 1) + c <- graphNode "C" (Vertex 2) + d <- graphNode "D" (Vertex 3) + e <- graphNode "E" (Vertex 4) + f <- graphNode "F" (Vertex 5) + g <- graphNode "G" (Vertex 6) + h <- graphNode "H" (Vertex 7) ++ newEdge a b () + newEdge a c ()+ newEdge b d ()+ newEdge c e () + newEdge d e ()+ newEdge d f ()+ newEdge e f ()+ newEdge c g ()+ newEdge e h ()+ newEdge g h ()+ newEdge g e ()+ + return ()+
+ Bayes/Examples/Tutorial.hs view
@@ -0,0 +1,244 @@+{- | Tutorial explaining how to make infereces with the library.++Thus tutorial is using examples from the module "Bayes.Examples". Please,+refer to this module for documentation about how the example bayesian networks are+created or loaded.++/Inferences/++The function 'inferencesOnStandardNetwork' is showing how to use variable elimination+and factor elimination to make inferences.++First, the 'example' is loaded to make its variables and its bayesian network available:++@+ let ([winter,sprinkler,rain,wet,road],exampleG) = 'example'+@++Then, we compute a prior marginal. Prior means that no evidence is used. A bayesian+network is a factorisation of a distribution P(A B C ...). If you want to know the+probability of only A, you need to sum out the other variables to eliminate them and get+P(A). To compute this prior marginal using variable elimnation, you need to give an elimination+order. The complexity of the computation is depending on the elimination order chosen.++For instance, if you want to compute the prior probability of rain, you can write:++@+ 'priorMarginal' exampleG [winter,sprinkler,wet,road] [rain] +@++Now, if you have observed that the grass is wet and want to take into account thios observation+to compute the posterior probability of rain (after observation):++@+ 'posteriorMarginal' exampleG [winter,sprinkler,wet,road] [rain] [wet '=:' True]+@ ++If you want to combine several observations:++@+ 'posteriorMarginal' exampleG [winter,sprinkler,wet,road] [rain] [wet '=:' True, sprinkler '=:' True]+@++There are several problems with variable elimination:++ * You have to specify an elimination order ++ * If you want to compute another marginal (for instance probability of winter), you have+ to recompute everything.++But, there exists another category of elimination algorithms based upon factor elimination. +They require the creation of an auxiliary data structure : the junction tree.++This tree is then used for computing all marginals (without having to recompute everything).+The junction tree is equivalent to giving an elimination order.++So, the previous examples can also be computed with factor elimination. First, the +junction tree must created:++@+ let jt = 'createJunctionTree' 'nodeComparisonForTriangulation' exampleG+@++The junction tree being equivalent to an elimination order, the order chosen will+depend on a cost function. In the previous example, the cost function 'nodeComparisonForTriangulation'+is used. Other cost functions may be introduced in a futute version of this library.++Once the junction tree has been computd, it can be used to compute several marginals:++@+ 'posterior' jt rain+@++The function is called posterior and will compute posterior only when solme evidence has+been introduced into the tree. Otherwise it is computing a prior.++To set evidence, you need to update the junction tree with new evidence:++@+ let jt' = 'updateEvidence' [wet '=:'' True] jt + 'posterior' jt' rain+@++/Inferences with an imported network/++There is a slight additional difficulty with imported networks : you need+to create new data type to be able to set evidence.++For instance, in the cancer network there is a Coma variable with levels Present or Absent.+When imported, those levels are imported as number. But, the evidence API in this library is+requiring enumerations.++So, you need to create a 'Coma' type:++@+ data Coma = Present | Absent deriving(Eq,Enum,Bounded)+@++and check that 'Present' is corresponding to the level 0 in the importd network.++Once this datatype is created, you can easily use the cancer network. First we load+the network and import the discrete variables of type 'DV' from the names of the nodes in the+network (not the label of the nodes)++@+ print \"CANCER NETWORK\"+ (varmap,cancer) <- 'exampleImport'+ print cancer+ let [varA,varB,varC,varD,varE] = fromJust $ mapM (flip Map.lookup varmap) ["A","B","C","D","E"]+@++Once the variables are available, you can create the junction tree and start making inferences:++@+ let jtcancer = 'createJunctionTree' 'nodeComparisonForTriangulation' cancer+--+ mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer $ x)) [varA,varB,varC,varD,varE]+--+ print \"UPDATED EVIDENCE\"+ let jtcancer' = 'updateEvidence' [varD '=:' Present] jtcancer + mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer' $ x)) [varA,varB,varC,varD,varE]+@++-}+module Bayes.Examples.Tutorial(+ -- * Tests with the standard network + inferencesOnStandardNetwork+ -- * Tests with the cancer network+ , inferencesOnCancerNetwork+ , Coma(..)+ , miscTest+ ) where ++import Bayes.Factor+import Bayes+import Bayes.VariableElimination+import Bayes.Examples(example, exampleJunction,exampleImport,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)+import Bayes.FactorElimination+import Data.Function(on)+import qualified Data.Map as Map+import Data.Maybe(fromJust,mapMaybe)+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++miscTest s = do + (varmap,perso) <- anyExample s+ let names = Map.keys varmap+ l = mapMaybe (flip Map.lookup varmap) names+ jtperso = createJunctionTree nodeComparisonForTriangulation perso+ print perso+ print jtperso+ print "FACTOR ELIMINATION"+ let post (v,name) = do + putStrLn name + print $ posterior jtperso v+ mapM_ post (zip l names)++ print "VARIABLE ELIMINATION"+ let prior (v,name) = do + putStrLn name + print $ priorMarginal perso (l L.\\ [v]) [v]+ mapM_ prior (zip l names)+++-- | Type defined to set the evidence on the Coma variable+-- from the cancer network.+data Coma = Present | Absent deriving(Eq,Enum,Bounded)++-- | Inferences with the cancer network+inferencesOnCancerNetwork = do + print "CANCER NETWORK"+ (varmap,cancer) <- exampleImport+ print cancer+ let [varA,varB,varC,varD,varE] = fromJust $ mapM (flip Map.lookup varmap) ["A","B","C","D","E"]+ let jtcancer = createJunctionTree nodeComparisonForTriangulation cancer++ 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 + 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 + mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varD,varE]++-- | Inferences with the standard network+inferencesOnStandardNetwork = do+ let ([winter,sprinkler,rain,wet,road],exampleG) = example++ putStrLn ""+ print "VARIABLE ELIMINATION"+ putStrLn ""+ print "Prior Marginal : probability of rain"+ let m = priorMarginal exampleG [winter,sprinkler,wet,road] [rain] + print m+ putStrLn ""++ print "Posterior Marginal : probability of rain if grass wet"+ let m = posteriorMarginal exampleG [winter,sprinkler,wet,road] [rain] [wet =: True]+ print m+ putStrLn ""++ print "Posterior Marginal : probability of rain if grass wet and sprinkler used"+ let m = posteriorMarginal exampleG [winter,sprinkler,wet,road] [rain] [wet =: True, sprinkler =: True]+ print m+ putStrLn ""++ let jt = createJunctionTree nodeComparisonForTriangulation exampleG++ putStrLn ""+ print "FACTOR ELIMINATION"+ putStrLn ""+ print "Prior Marginal : probability of rain"+ let m = posterior jt rain+ print m+ putStrLn ""++ let jt' = updateEvidence [wet =: True] jt ++ print "Posterior Marginal : probability of rain if grass wet"+ let m = posterior jt' rain+ print m+ putStrLn ""++ let jt'' = clearEvidence jt'+ print "Prior Marginal : probability of rain"+ let m = posterior jt rain+ print m+ putStrLn ""++ let jt3 = updateEvidence [wet =: True, sprinkler =: True] jt'++ print "Posterior Marginal : probability of rain if grass wet and sprinkler used"+ let m = posterior jt3 rain+ print m+ putStrLn ""++ return ()
+ Bayes/Factor.hs view
@@ -0,0 +1,610 @@+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{- | Conditional probability table++Conditional Probability Tables and Probability tables++-}+module Bayes.Factor(+ -- * Factor+ Factor(..)+ , isomorphicFactor+ , normedFactor+ -- * Set of variables + , Set(..)+ , BayesianDiscreteVariable(..)+ -- * Implementation+ , Vertex(..)+ -- ** Discrete variables and instantiations+ , DV(..)+ , DVSet(..)+ , DVI+ , DVISet(..)+ , setDVValue+ , instantiationValue+ , instantiationVariable+ , variableVertex+ , (=:)+ , forAllInstantiations+ , factorFromInstantiation+ , changeVariableOrder+ -- ** Factor+ , CPT+ -- * Tests+ , testProductProject_prop+ , testScale_prop+ , testProjectCommut_prop+ , testScalarProduct_prop+ , testProjectionToScalar_prop+ ) where++import qualified Data.Vector.Unboxed as V+import Data.Vector.Unboxed((!))+import Data.Maybe(fromJust,mapMaybe)+import qualified Data.List as L+import Text.PrettyPrint.Boxes hiding((//))+import Test.QuickCheck+import Test.QuickCheck.Arbitrary+import qualified Data.IntMap as IM+import Control.Monad+import System.Random(Random)++--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+factorFromInstantiation (DVI dv a) = + let setValue i = if i == a then 1.0 else 0.0 + 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)++-- | 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+ -- | 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+ emptyFactor :: f+ -- | 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 + -- | 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+ -- | 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+ -- | Value of factor for a given set of variable instantitation.+ -- The variable instantion is like a multi-dimensional index.+ factorValue :: f -> DVISet 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+ -- | Dimension of the factor (number of floating point values)+ factorDimension :: f -> Int+ + -- | Norm of the factor = sum of its values+ factorNorm :: f -> Double + ++ -- | Scale the factor values by a given scaling factor+ factorScale :: Double -> f -> f++ -- | Create a scalar factor with no variables+ factorFromScalar :: Double -> f++ -- | Create an evidence factor from an instantiation.+ -- If the instantiation is empty then we get nothing+ evidenceFrom :: DVISet Int -> Maybe f+ ++ -- | Divide all the factor values+ factorDivide :: f -> Double -> f+ factorDivide f d = (1.0 / d) `factorScale` f ++ -- | 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+ -- | Project to. The variable are kept and other variables are removed+ factorProjectTo :: DVSet -> 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 + -> [Double] -- ^ Old values+ -> [Double] -- ^ New values+changeVariableOrder 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.+-- The variables on the left or right of the condition bar are not tracked. What's matter is that+-- it is encoding a function of several variables.+-- Marginalization of variables will be computed from the bayesian graph where+-- 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+ }+ | Scalar Double++debugCPT (Scalar d) = do + putStrLn "SCALAR CPT"+ print d+ putStrLn ""++debugCPT (CPT d m v) = do + putStrLn "CPT"+ print d + putStrLn ""+ print m + putStrLn ""+ print v+ putStrLn ""+{-++CPT can't have same same vertex values but with different sizes.+But, arbitrary CPT generation will general several vertex with same vertex id+and different vertex size.++So, we introduce a function mapping a vertex ID to a vertex size. So, vertex size are hard coded++-}++quickCheckVertexSize :: Int -> Int+quickCheckVertexSize 0 = 2+quickCheckVertexSize 1 = 2+quickCheckVertexSize 2 = 2+quickCheckVertexSize _ = 2++-- | Generate a random value until this value is not already present in the list+whileIn :: (Arbitrary a, Eq a) => [a] -> Gen a -> Gen a+whileIn l m = do + newVal <- m + if newVal `elem` l + then+ whileIn l m + else + return newVal++-- | Generate a random vector of n elements without replacement (no duplicate)+-- May loop if the range is too small !+generateWithoutReplacement :: (Random a, Arbitrary a, Eq a) + => Int -- ^ Vector size+ -> (a,a) -- ^ Bounds+ -> Gen [a]+generateWithoutReplacement n b | n == 1 = generateSingle b + | n > 1 = generateMultiple n b + | otherwise = return []+ where+ generateSingle b = do + r <- choose b+ return [r]+ generateMultiple n b = do + l <- generateWithoutReplacement (n-1) b+ newelem <- whileIn l $ choose b+ return (newelem:l)++++instance Arbitrary CPT where+ arbitrary = do + nbVertex <- choose (1,4) :: Gen Int+ vertexNumbers <- generateWithoutReplacement nbVertex (0,50)+ let dimensions = map (\i -> (DV (Vertex i) (quickCheckVertexSize i))) vertexNumbers+ let valuelen = product (map dimension dimensions)+ rndValues <- vectorOf valuelen (choose (0.0,1.0) :: Gen Double)+ return . fromJust . factorWithVariables dimensions $ rndValues++-- | Test product followed by a projection when the factors have no+-- common variables++-- | Floating point number comparisons which should take into account+-- all the subtleties of that kind of comparison+nearlyEqual :: Double -> Double -> Bool+nearlyEqual a b = + let absA = abs a + absB = abs b + diff = abs (a-b)+ epsilon = 2e-5+ in+ case (a,b) of + (x,y) | x == y -> True -- handle infinities+ | x*y == 0 -> diff < (epsilon * epsilon)+ | otherwise -> diff / (absA + absB) < epsilon++testScale_prop :: Double -> CPT -> Bool+testScale_prop s f = (factorNorm (s `factorScale` f)) `nearlyEqual` (s * (factorNorm f))++testProductProject_prop :: CPT -> CPT -> Property+testProductProject_prop fa fb = isEmpty ((factorVariables fa) `intersection` (factorVariables fb)) ==> + let r = factorProjectOut (factorVariables fb) (factorProduct [fa,fb])+ fa' = r `factorDivide` (factorNorm fb)+ in+ fa' `isomorphicFactor` fa++testScalarProduct_prop :: Double -> CPT -> Bool +testScalarProduct_prop v f = (factorProduct [(Scalar v),f]) `isomorphicFactor` (v `factorScale` f)++testProjectionToScalar_prop :: CPT -> Bool +testProjectionToScalar_prop f = + let allVars = factorVariables f + in+ (factorProjectOut allVars f) `isomorphicFactor` (factorFromScalar (factorNorm f))++testProjectCommut_prop:: CPT -> Property +testProjectCommut_prop f = length (factorVariables f) >= 3 ==>+ let a = take 1 (factorVariables f)+ b = take 1 . drop 1 $ factorVariables f + commuta = factorProjectOut a (factorProjectOut b f)+ commutb = factorProjectOut b (factorProjectOut a f)+ in+ commuta `isomorphicFactor` commutb++-- | Test equality of two factors taking into account the fact+-- that the variables may be in a different order.+-- In case there is a distinction between conditionned variable and+-- conditionning variables (imposed from the exterior) then this+-- comparison may not make sense. It is a comparison of+-- function of several variables which no special interpretation of the+-- 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)+ guard (factorDimension fa == factorDimension fb)+ guard $ and [factorValue fa ia `nearlyEqual` factorValue fb ia | ia <- forAllInstantiations va]+ return ()++{-++Following functions are used to typeset the factor when displaying it++-}+vname :: Int -> Int -> Box+vname vc i = text $ "v" ++ show vc ++ "=" ++ show i++dispFactor :: FactorPrivate f => f -> DV -> [Int] -> DVSet -> Box+dispFactor cpt h c [] = + let dstIndexes = allValues h+ dependentIndexes = reverse c+ factorValueAtPosition p = + let v = factorValuePrivate 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)++instance Show CPT where+ show (Scalar v) = "\nScalar Factor:\n" ++ show v+ show c@(CPT [] _ v) = "\nEmpty CPT:\n"++ 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)+ in+ "\n" ++ show d ++ "\n" ++ render (hsep 1 top [dstColumn,table])++instance Factor CPT where+ emptyFactor = emptyCPT+ isScalarFactor (Scalar _) = True+ isScalarFactor _ = False+ factorFromScalar v = Scalar v+ factorDimension f@(CPT _ _ _) = product . map dimension . factorVariables$ f+ factorDimension _ = 1+ containsVariable (CPT _ m _) (DV (Vertex i) _) = IM.member i m+ containsVariable (Scalar _) _ = False+ factorWithVariables = createCPTWithDims+ factorVariables (CPT v _ _) = v+ factorVariables (Scalar _) = []+ factorNorm f@(CPT _ _ _) = sum [ factorValuePrivate f x | x <- forAllIndices (factorVariables f)]+ 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)]+ in + fromJust $ factorWithVariables (factorVariables f) newValues+ factorValue (Scalar v) _ = v + factorValue f i = + let multiIndex = reorder i f+ in + factorValuePrivate f multiIndex+ evidenceFrom [] = Nothing + evidenceFrom l = + let index = map instantiationValue l + variables = map instantiationVariable l+ setValueForIndex i = if i == index then 1.0 else 0.0 + in+ factorWithVariables variables . map setValueForIndex $ forAllIndices variables++instance FactorPrivate CPT where+ factorValuePrivate = getCPTValue+++-- | 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+ pos' = scanr (*) (1::Int) (tail dim)+ c = sum . map (\(x,y) -> x * y) $ (zip pos' pos)+ in + c++-- | 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 dims values = + let createDVIndex i (DV (Vertex v) _) = (v,i)+ m = IM.fromList . zipWith createDVIndex ([0,1..]::[Int]) $ dims+ p = product (map dimension dims) + in+ if length values == p+ then+ Just $ CPT dims m (V.fromList values)+ else + Nothing+
+ Bayes/FactorElimination.hs view
@@ -0,0 +1,654 @@+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{- | Algorithms for factor elimination++-}+module Bayes.FactorElimination(+ -- * Moral graph+ moralGraph+ -- * Triangulation+ , nodeComparisonForTriangulation+ , numberOfAddedEdges+ , triangulate+ -- * Junction tree+ , minimumSpanningTree+ , createClusterGraph+ , Cluster+ , createJunctionTree+ , JunctionTree+ -- * Shenoy-Shafer message passing+ , collect + , distribute+ , posterior + -- * Evidence+ , clearEvidence+ , updateEvidence+ -- * Test + , junctionTreeProperty_prop+ , createVerticesJunctionTree+ , VertexCluster+ ) where++import Bayes+import qualified Data.Foldable as F+import Data.Maybe(fromJust,mapMaybe,isJust)+import Control.Monad(mapM)+import Bayes.Factor hiding (isEmpty)+import Data.Function(on)+import Data.List(minimumBy,maximumBy,inits)+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 Test.QuickCheck hiding ((.||.), collect)+import Test.QuickCheck.Arbitrary++--import Debug.Trace+--debug s a = trace (s ++ " " ++ show a ++ "\n") a++{-+ +Comparison functions for graph triangulation++-}++-- | Number of edges added when connecting all neighbors+numberOfAddedEdges :: UndirectedGraph g + => g a b + -> Vertex + -> Int +numberOfAddedEdges g v = + let nodes = fromJust $ neighbors g v+ in + length [edge x y | x <- nodes, y <- nodes, x /= y, not (isLinkedWithAnEdge g x y)]++-- | Weight of a node+weight :: (UndirectedGraph g, Factor f)+ => g a f + -> Vertex + -> Int +weight g v = + factorDimension . fromJust . vertexValue g $ v++(.||.) :: (a -> a -> Ordering)+ -> (a -> a -> Ordering) + -> (a -> a -> Ordering)+f .||. g = + \a b -> case f a b of+ EQ -> g a b + r -> r++-- | Node selection comparison function used for triangulating the graph+nodeComparisonForTriangulation :: (UndirectedGraph g, Factor f)+ => g a f+ -> Vertex + -> Vertex + -> Ordering +nodeComparisonForTriangulation g = (compare `on` (numberOfAddedEdges g)) .||. (compare `on` (weight g))++{-++Graph triangulation++-}++-- | A cluster containing only the vertices and not yet the factors+newtype VertexCluster = VertexCluster (Set.Set Vertex) deriving(Eq)++fromVertexCluster (VertexCluster s) = s++instance Show VertexCluster where + show (VertexCluster s) = show . Set.toList $ s++-- | Triangulate a graph using a cost function+-- The result is the triangulated graph and the list of clusters+-- which may not be maximal.+triangulate :: Graph g+ => (Vertex -> Vertex -> Ordering) -- ^ Criterion function for triangulation+ -> g () b+ -> ([VertexCluster],g () b) -- ^ Returns the clusters and the triangulated graph+triangulate cmp g = + -- At start, gsrc and gdst are the same+ -- gsrc is modified. It is where vertex elimination is taking place.+ -- The edges are added to gdst+ let processAllNodes gsrc gdst l | hasNoVertices gsrc = (keepMaximalClusters (reverse l),gdst)+ | otherwise = + let selectedNode = minimumBy cmp (allVertices gsrc)+ theNeighbors = selectedNode : (fromJust $ neighbors gsrc selectedNode)+ addEmptyEdge e g = addEdge e () g+ (gsrc',gdst') = connectAllNodesWith gsrc gdst addEmptyEdge addEmptyEdge theNeighbors+ gsrc'' = removeVertex selectedNode gsrc' + in + processAllNodes gsrc'' gdst' ((VertexCluster . Set.fromList $ theNeighbors) : l)++ in + processAllNodes g g []+++-- | Find for a containing cluster. +findContainingCluster :: VertexCluster -- ^ Cluster processed+ -> [VertexCluster] -- ^ Cluster list where to look for a containing cluster+ -> (Maybe VertexCluster,[VertexCluster]) -- ^ Return the containing cluster and a new list without the containing cluster+findContainingCluster cluster l = + let clusterIsNotASubsetOf s = (Set.isSubsetOf (fromVertexCluster cluster) (fromVertexCluster s))+ (prefix,suffix) = break clusterIsNotASubsetOf l+ in + case suffix of + [] -> (Nothing,l)+ _ -> (Just (head suffix),prefix ++ tail suffix)+++-- | Remove clusters already contained in a previous clusters+keepMaximalClusters :: [VertexCluster] -> [VertexCluster]+keepMaximalClusters [] = []+keepMaximalClusters l = checkIfMaximal [] (head l) (tail l)+ where + checkIfMaximal reversedPrefix current [] = + case findContainingCluster current (reverse reversedPrefix) of + (Nothing,_) -> reverse (current:reversedPrefix) + (Just r,l) -> reverse (r:reverse l)+ checkIfMaximal reversedPrefix current suffix = + case findContainingCluster current (reverse reversedPrefix) of + (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 + -> VertexCluster + -> Cluster +vertexClusterToCluster g c = + let vertices = Set.toList . fromVertexCluster $ c+ variables = map factorMainVariable . mapMaybe (vertexValue g) $ vertices+ 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+ in + evidenceFrom selectedVariables+++-- | 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 ++downMessage (Distribute _ _ d) = Just d +downMessage _ = Nothing ++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+ in+ vertexClusterToCluster g separatorVertices++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)++setFactorEdgeUpdate :: (Graph g, Factor f) + => BayesianNetwork g f + -> VertexCluster + -> VertexCluster+ -> () + -> Separator f+setFactorEdgeUpdate g parentvalue childvalue _ = NoMessage $ computeSeparatorCluster g parentvalue childvalue ++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++-- | 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) +++-- | 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 = + let theMoralGraph = moralGraph g+ (clusters,_) = triangulate (cmp theMoralGraph) theMoralGraph+ g'' = createClusterGraph clusters :: UndirectedSG Int VertexCluster+ in + minimumSpanningTree g''++-- | Create a function tree+createJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f, Show f)+ => (UndirectedSG () f -> Vertex -> Vertex -> Ordering) -- ^ Weight function on the moral graph+ -> 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+ -- 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+ 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+++-- | 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 ++-- | 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')++-- | 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 g = (not . isEmpty) g && (not . hasNoEdges) g && connectedGraph g ==> + let cmp ug = (compare `on` (numberOfAddedEdges ug))+ in+ junctionTreeProperty [] (createVerticesJunctionTree cmp g)++-- | 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+ in + and $ zipWith isSubsetOfAllParents allIntersections pathsToEachParent+{-++Moral graph++-}+-- | Get the parents of a vertex+parents :: DirectedGraph g => g a b -> Vertex -> [Vertex]+parents g v = fromJust $ ingoing g v >>= mapM (startVertex g) ++-- | Get the children of a vertex+children :: DirectedGraph g => g a b -> Vertex -> [Vertex]+children g v = fromJust $ outgoing g v >>= mapM (endVertex g) ++-- | Connect all the nodes which are not connected and apply the function f for each new connection+-- The origin and dest graph must share the same vertex.+connectAllNodesWith :: (Graph g, Graph g') + => g a b -- ^ Graph containing the nodes+ -> g' a b -- ^ Graph to be modified+ -> (Edge -> g a b -> g a b) -- ^ Function used to modify the source graph+ -> (Edge -> g' a b -> g' a b) -- ^ Function used to modify a new graph+ -> [Vertex] -- ^ List of nodes to connect+ -> (g a b,g' a b) -- ^ Result graph+connectAllNodesWith originGraph dstGraph g f nodes = + let h e (x,y) = (g e x, f e y)+ (originGraph',dstGraph') = + foldr h (originGraph,dstGraph) [edge x y | x <- nodes, y <- nodes, x /= y, not (isLinkedWithAnEdge originGraph x y)]+ in + (originGraph',dstGraph')++-- | Add the missing parent links+addMissingLinks :: DirectedGraph g => Vertex -> b -> g () b -> g () b+addMissingLinks v _ g = + let (_,g') = connectAllNodesWith g g (\e m -> m) (\e m -> addEdge e () m) (parents g v)+ in + g'+++-- | Convert the graph to an undirected form+convertToUndirected :: (FoldableWithVertex g, Graph g, NamedGraph g, NamedGraph g',UndirectedGraph g')+ => g () b + -> g' () b +convertToUndirected m = + let addVertexWithLabel v dat g = + let theName = fromJust $ vertexLabel m v+ in + addLabeledVertex theName v dat g+ newDiscreteGraph = foldrWithVertex addVertexWithLabel emptyGraph m+ addEmptyEdge edge g = addEdge edge () g+ in + foldr addEmptyEdge newDiscreteGraph . allEdges $ m++-- | For the junction tree construction, only the vertices are needed during the intermediate steps.+-- So, the moral graph is returned without any vertex data.+moralGraph :: (NamedGraph g, FoldableWithVertex g, DirectedGraph g) + => g () b -> UndirectedSG () b +moralGraph g = + convertToUndirected . foldrWithVertex addMissingLinks g $ g
+ Bayes/ImportExport/HuginNet.hs view
@@ -0,0 +1,192 @@+-- | Parser for a subset of the Hugin Net language+module Bayes.ImportExport.HuginNet( + importBayesianGraph+ ) where ++import Text.ParserCombinators.Parsec.Prim+import Text.ParserCombinators.Parsec.Char +import Text.ParserCombinators.Parsec.Combinator ++import Data.Maybe(mapMaybe,fromJust)+import Bayes.ImportExport.HuginNet.Splitting +import qualified Data.Map as Map+import Bayes.Factor+import Bayes++--import Debug.Trace ++--debug a = trace (show a) a++data Section = Net + | Node String [String] Int+ | Potential [String] [String]+ deriving(Eq,Show)++name :: Parser String +name = many1 (alphaNum <|> oneOf "_-")++sectionContent :: Parser ()+sectionContent = do + string "{"+ newline+ many1 (noneOf "}")+ string "}" + optional newline+ return ()++net :: Parser Section+net = do + string "net"+ newline+ sectionContent+ return Net++levelName = do + char '"'+ s <- many1 (noneOf "\"")+ char '"'+ return s ++-- | Node states+state :: Parser [String]+state = do + spaces+ string "states"+ spaces+ char '='+ spaces + char '('+ spaces+ levels <- sepEndBy1 levelName (many1 space)+ char ')'+ spaces+ char ';'+ spaces+ optional newline+ return levels++factorValues :: Parser String+factorValues = do + spaces + string "data"+ spaces+ char '='+ spaces+ r <- many1 (noneOf ";")+ spaces + optional newline + return r++unknownCommand = do + manyTill (noneOf "}") newline + return Nothing++recognizedCommand :: Parser a -> Parser (Maybe a)+recognizedCommand c = choice [try c >>= return . Just, unknownCommand]++node :: Parser Section+node = do + string "node"+ spaces+ n <- name+ newline+ string "{"+ newline+ l <- many (recognizedCommand state)+ string "}" + optional newline+ let r = concat . mapMaybe id $ l+ return $ Node n r (length r)++potential :: Parser Section+potential = do + string "potential"+ spaces + conditions <- manyTill anyChar newline+ string "{"+ newline+ l <- many (recognizedCommand factorValues)+ string "}" + optional newline+ let r = concat . mapMaybe id $ l+ return $ Potential (splitCPT conditions) (splitValues r)++section :: Parser Section+section = choice [try net,try node,try potential]++comment = do + string "%%"+ manyTill anyChar newline+ return () ++manyEmpty = skipMany (space <|> newline)++netParser :: Parser [Section]+netParser = do+ many comment+ manyEmpty+ sepEndBy1 section manyEmpty++addVertexName (Node s _ d) (c,m) = (c+1,Map.insert s (DV (Vertex c) d) m)+addVertexName _ (c,m) = (c,m)++addSection m (Node _ _ _) = return ()++addSection m (Net) = return ()+addSection m (Potential conditions values) = do + let dvs = fromJust . mapM (flip Map.lookup m) $ conditions+ dst = head dvs + conds = tail dvs+ oldOrder = conds ++ [dst]+ dvalues = map read values :: [Double]+ newvalues = changeVariableOrder oldOrder dvs dvalues+ cpt dst conds ~~ newvalues+ return ()++addVariables (Node s _ d) = do + v <- variableWithSize s d+ return $ Just (s,v)++addVariables _ = return Nothing++-- | Import a bayesian network form a Hugin file.+-- Only a subset of the file format is supported.+-- You may have to convert the line endings to be able to parse a file+-- When it is succeeding, it is returing a bayesian network monad and+-- a mapping from node names to discrete variables.+importBayesianGraph :: Factor f + => String + -> IO (Maybe (BNMonad DirectedSG f (Map.Map String DV)))+importBayesianGraph s = do + r <- readBayesianNetwork s + case r of + Nothing -> return Nothing + Just s -> return . Just $ createBayesianGraph s++mapMaybeM :: Monad m => (a -> m (Maybe b)) -> [a] -> m [b]+mapMaybeM f l = mapM f l >>= return . mapMaybe id++createBayesianGraph :: Factor f => [Section] -> BNMonad DirectedSG f (Map.Map String DV)+createBayesianGraph s = do + vars <- mapMaybeM addVariables s+ let m = Map.fromList vars+ mapM_ (addSection m) s+ return m++-- | Horrible way to remove the comments+filterComment :: Bool -> String -> String+filterComment False ('%':l) = filterComment True l+filterComment False (a:l) = a:filterComment False l +filterComment False [] = []+filterComment True ('\n':l) = '\n':filterComment False l +filterComment True (a:l) = filterComment True l +filterComment True [] = []++readBayesianNetwork s = do + f <- readFile s+ let result = runParser netParser () s (filterComment False f)+ case result of + Left err -> do + print err + return Nothing+ Right a -> return (Just a)
+ Bayes/ImportExport/HuginNet/Splitting.hs view
@@ -0,0 +1,9 @@+module Bayes.ImportExport.HuginNet.Splitting ( + splitCPT+ , splitValues+ ) where ++import Data.List.Split++splitCPT = split (dropBlanks . dropDelims $ oneOf "() |") +splitValues = split (dropBlanks . dropDelims $ oneOf "() \n\t")
+ Bayes/Test.hs view
@@ -0,0 +1,38 @@+{- | Testing of the implementation.++-}+module Bayes.Test (+ runTests+ ) where+import Test.Framework (defaultMain, testGroup)+import Test.Framework.Providers.QuickCheck2 (testProperty)+import Test.Framework.Providers.HUnit(testCase)+import Bayes.Test.CompareEliminations(compareVariableFactor)++import Bayes(testEdgeRemoval_prop,testVertexRemoval_prop)+import Bayes.Factor(testProductProject_prop,testScale_prop,testProjectCommut_prop,testScalarProduct_prop,testProjectionToScalar_prop)+import Bayes.FactorElimination(junctionTreeProperty_prop)++-- | Run all the tests+runTests = defaultMain tests++tests = [+ testGroup "Graph" [+ testProperty "Edge Removal" testEdgeRemoval_prop,+ testProperty "Vertex Removal" testVertexRemoval_prop+ ]+ , testGroup "Factor" [+ testProperty "Factor scaling and norm" testScale_prop,+ testProperty "Product / Project" testProductProject_prop,+ testProperty "Commutativity of project" testProjectCommut_prop,+ testProperty "Product with scalar factor" testScalarProduct_prop,+ testProperty "Test projection to scalar" testProjectionToScalar_prop+ ]+ , testGroup "Junction Tree" [+ testProperty "Test the junction tree property" junctionTreeProperty_prop,+ testCase "Test variable elimination == factor elimination" compareVariableFactor+ ]++ ]++
+ Bayes/Test/CompareEliminations.hs view
@@ -0,0 +1,47 @@+{- | A comparison of variable elimination and factor elimination on a simple graph.++It is a non regression test.++-}+module Bayes.Test.CompareEliminations(+ compareVariableFactor+ ) where++import Test.HUnit.Lang(assertFailure)++import Bayes.Examples(example)+import Bayes.Factor+import Bayes+import Bayes.VariableElimination+import Bayes.FactorElimination++compareFactors :: String -> Maybe CPT -> CPT -> IO ()+compareFactors s Nothing _ = assertFailure s+compareFactors s (Just a) b = + if a `isomorphicFactor` b + then + return () + else + assertFailure s++-- | Compare that variable elemination and factor elimination are giving+-- similar results on a simple example+compareVariableFactor :: IO ()+compareVariableFactor = do + let ([winter,sprinkler,rain,wet,road],exampleG) = example+ 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 ++ compareFactors "POSTERIOR RAIN FOR WET" (posterior jt1 rain) + (posteriorMarginal exampleG [winter,sprinkler,wet,road] [rain] [wet =: True])+ compareFactors "POSTERIOR RAIN FOR WET" (posterior jt2 rain) + (posteriorMarginal exampleG [winter,sprinkler,wet,road] [rain] [wet =: True, sprinkler =: True])++ compareFactors "PRIOR FOR WINTER" (posterior jt winter) (priorMarginal exampleG [sprinkler,wet,road,rain] [winter])+ compareFactors "PRIOR FOR SPRINKLER" (posterior jt sprinkler) (priorMarginal exampleG [winter,wet,road,rain] [sprinkler])+ compareFactors "PRIOR FOR WET" (posterior jt wet) (priorMarginal exampleG [winter,sprinkler,road,rain] [wet])+ compareFactors "PRIOR FOR ROAD" (posterior jt road) (priorMarginal exampleG [winter,sprinkler,wet,rain] [road])+
+ Bayes/VariableElimination.hs view
@@ -0,0 +1,219 @@+{- | Algorithms for variable elimination++-}+module Bayes.VariableElimination(+ -- * Inferences+ priorMarginal+ , posteriorMarginal+ -- * Interaction graph and elimination order+ , interactionGraph+ , degreeOrder+ , minDegreeOrder+ , minFillOrder+ , allVariables+ , EliminationOrder+ ) where++import Bayes+import Bayes.Factor+import Data.List(partition,minimumBy,(\\),find)+import Data.Maybe(fromJust)+import Data.Function(on)+import qualified Data.Map as M++--import Debug.Trace ++--debug s a = trace (s ++ "\n" ++ show a ++ "\n") a++-- | Elimination order+type EliminationOrder = DVSet++-- | Get all variables from a Bayesian Network+allVariables :: (Graph g, Factor f) + => BayesianNetwork g f + -> DVSet+allVariables g = + let s = allVertexValues g + createDV = factorMainVariable + in + map createDV s++-- | Used for bucket elimination. Factor are organized by their first DV+type Buckets f = (EliminationOrder,M.Map DV [f])++createBuckets :: (Graph g, Factor f, Show f) + => BayesianNetwork g f -- ^ Bayesian Network+ -> EliminationOrder -- ^ Variables to eliminate+ -> EliminationOrder -- ^ Remaining variables+ -> Buckets f +createBuckets g e r = + let s = allVertexValues g+ -- We put the selected variables for elimination in the right order at the beginning+ -- Which means the function can work with a partial order which is completed with other+ -- variables by default.+ theOrder = e ++ r+ addDVToBucket dv (rf, m) =+ let (fk,remaining) = partition (flip containsVariable dv) rf+ in + (remaining, M.insert dv fk m)+ (_,b) = foldr addDVToBucket (s,M.empty) (reverse theOrder)+ in+ (tail theOrder,b)++-- | Get the factors for a bucket+getBucket :: DV + -> Buckets f + -> [f]+getBucket dv (_,m) = fromJust $ M.lookup dv m++-- | Update bucket+updateBucket :: Factor f => DV -> f -> Buckets f -> Buckets f +updateBucket dv f b@(e,m) = + if isScalarFactor f + then + (tail e,M.insert dv [f] m)+ else+ let b' = removeFromBucket dv b+ (e',m') = addBucket f b'+ in + (tail e',m')++-- | Add a factor to the right bucket+addBucket :: Factor f => f -> Buckets f -> Buckets f+addBucket f (e,b) = + let inBucket = find (f `containsVariable`) e+ in + case inBucket of + Nothing -> (e,b)+ Just bucket -> (e, M.insertWith' (++) bucket [f] b)++-- | Remove a variable from the bucket+removeFromBucket :: DV -> Buckets f -> Buckets f +removeFromBucket dv (e,m) = (e,M.delete dv m) ++-- | Compute the prior marginal. All the variables in the+-- elimination order are conditionning variables ( p( . | conditionning variables) )+posteriorMarginal :: (Graph g, Factor f, Show f) + => BayesianNetwork g f -- ^ Bayesian Network+ -> EliminationOrder -- ^ Ordering of variables to marginzalie+ -> EliminationOrder -- ^ Ordering of remaining variables+ -> [DVI Int] -- ^ Assignment for some factors in vaiables to marginalize+ -> f+posteriorMarginal n p r assignment = + -- The elimintation order are the variables to eliminate.+ -- But the algorithm also needs the remaining variables+ let bucket = createBuckets n p r+ assignmentFactors = map factorFromInstantiation assignment+ bucket' = foldr addBucket bucket assignmentFactors+ (_,resultBucket) = foldr marginalizeOneVariable bucket' (reverse p)+ resultFactor = factorProduct . concat . M.elems $ resultBucket+ -- The norm is P(e) and result factor is P(Q,e)+ norm = factorNorm resultFactor+ in+ -- We get P(Q | e)+ resultFactor `factorDivide` norm + where + marginalizeOneVariable dv currentBucket = + let fk = getBucket dv currentBucket+ p = factorProduct fk+ f' = factorProjectOut [dv] p+ in+ updateBucket dv f' currentBucket++-- | Compute the prior marginal. All the variables in the+-- elimination order are conditionning variables ( p( . | conditionning variables) )+priorMarginal :: (Graph g, Factor f, Show f) + => BayesianNetwork g f -- ^ Bayesian Network+ -> EliminationOrder -- ^ Ordering of variables to marginalize+ -> EliminationOrder -- ^ Ordering of remaining to keep in result+ -> f+priorMarginal g ea eb = posteriorMarginal g ea eb []++-- | Compute the interaction graph of the BayesianNetwork+interactionGraph :: (FoldableWithVertex g,Factor f, UndirectedGraph g')+ => BayesianNetwork g f+ -> g' () DV+interactionGraph g = + foldrWithVertex addFactor emptyGraph g + where+ addFactor vertex factor graph = + let allvars = factorVariables factor+ edges = [(x,y) | x <- allvars, y <- allvars , x /= y]+ addNewEdge (va,vb) g = + let g' = addVertex (variableVertex vb) vb . addVertex (variableVertex va) va $ g + in+ addEdge (edge (variableVertex va) (variableVertex vb)) () $ g'+ in + foldr addNewEdge graph edges++-- | Number of neighbors for a variable in the bayesian network+nbNeighbors :: UndirectedSG () DV + -> DV + -> Int +nbNeighbors g dv = + let r = fromJust $ neighbors g (variableVertex dv)+ in + length r++-- | Number of missing links between the neighbors of the graph+nbMissingLinks :: UndirectedSG () DV + -> DV + -> Int +nbMissingLinks g dv = + let r = fromJust $ neighbors g (variableVertex dv)+ edges = [(x,y) | x <- r, y <- r , x /= y, not (isLinkedWithAnEdge g x y)]+ in + length edges++-- | Compute the degree order of an elimination order+degreeOrder :: (FoldableWithVertex g, Factor f, Graph g)+ => BayesianNetwork g f+ -> EliminationOrder + -> Int +degreeOrder g p =+ let ig = interactionGraph g :: UndirectedSG () DV+ (_,w) = foldr processVariable (ig,0) p + in + w + where + addAnEdge (va,vb) g = addEdge (edge va vb) () g+ processVariable bdv (g,w) = + let r = fromJust $ neighbors g (variableVertex bdv)+ nbNeighbors = length r+ edges = [(x,y) | x <- r, y <- r , x /= y, not (isLinkedWithAnEdge g x y)]+ g' = removeVertex (variableVertex bdv) (foldr addAnEdge g edges)+ in+ if nbNeighbors > w + then + (g',nbNeighbors) + else + (g',w)+ +-- | Find an elimination order minimizing a metric+eliminationOrderForMetric :: (Graph g, Factor f, FoldableWithVertex g, UndirectedGraph g')+ => (g' () DV -> DV -> Int)+ -> BayesianNetwork g f + -> EliminationOrder +eliminationOrderForMetric metric g = + let ig = interactionGraph g+ s = allVertexValues ig+ getOptimalNode _ [] = []+ getOptimalNode g l = + let (optimalNode,_) = minimumBy (compare `on` snd) . map (\v -> (v,metric g v)) $ l+ g' = removeVertex (variableVertex optimalNode) g+ in + optimalNode : getOptimalNode g' (l \\ [optimalNode])+ in + getOptimalNode ig s++-- | Elimination order minimizing the degree+minDegreeOrder :: (Graph g, Factor f, FoldableWithVertex g)+ => BayesianNetwork g f + -> EliminationOrder +minDegreeOrder = eliminationOrderForMetric nbNeighbors++-- | Elimination order minimizing the filling+minFillOrder :: (Graph g, Factor f, FoldableWithVertex g)+ => BayesianNetwork g f + -> EliminationOrder +minFillOrder = eliminationOrderForMetric nbMissingLinks
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c)2012, alpheccar++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 alpheccar nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,4 @@+import Distribution.Simple+++main = defaultMain
+ cancer.net view
@@ -0,0 +1,128 @@+net+{+ HR_Compile_TriangMethod = "0";+ HR_Monitor_GraphPrecision = "100";+ HRUNTIME_Grid_GridSnap = "1";+ HRUNTIME_Propagate_AutoSum = "1";+ HR_Propagate_AutoNormal = "1";+ HRUNTIME_Monitor_OpenGraph = "0";+ HR_Font_Size = "-12";+ HR_Monitor_AutoUpdGraph = "0";+ node_size = (100.0 40.0);+ HRUNTIME_Propagate_AutoNormal = "1";+ HR_Grid_GridSnap = "1";+ HR_Compile_Compress = "0";+ HRUNTIME_Compile_Compress = "0";+ HR_Compile_ApproxEpsilon = "0.00001";+ jenginegenerator6060830225489488864L = "edu.ucla.belief.inference.JoinTreeSettings@7d399ae5";+ HR_Color_DiscreteChance = "16";+ HR_Propagate_AutoSum = "1";+ HR_Propagate_Auto = "0";+ HR_Compile_Approximate = "0";+ HRUNTIME_Monitor_InitStates = "5";+ HR_Grid_GridShow = "0";+ HRUNTIME_Font_Name = "Arial";+ HR_Groups_GroupColors = "";+ HR_Groups_GroupNames = "";+ HR_Color_ContinuosChance = "48";+ HR_Groups_UserGroupsNo = "0";+ HRUNTIME_Grid_GridShow = "0";+ HRUNTIME_Propagate_Auto = "0";+ HR_Color_Decision = "17";+ HR_Monitor_InitSD = "2";+ HRUNTIME_Grid_X = "10";+ HRUNTIME_Grid_Y = "10";+ HR_Grid_X = "10";+ HRUNTIME_Compile_TriangMethod = "0";+ HRUNTIME_Font_Size = "-12";+ HR_Grid_Y = "10";+ HR_Font_Name = "Arial";+ HR_Font_Weight = "400";+ HR_Monitor_InitStates = "5";+ HRUNTIME_Monitor_AutoUpdGraph = "0";+ HR_Font_Italic = "0";+ HR_Monitor_OpenGraph = "0";+ HRUNTIME_Font_Weight = "400";+ HRUNTIME_Font_Italic = "0";+ HRUNTIME_Compile_Approximate = "0";+ HR_Color_Utility = "36";+ HRUNTIME_Monitor_GraphPrecision = "100";+ HRUNTIME_Compile_ApproxEpsilon = "0.00001";+}++node D+{+ states = ("Present" "Absent" );+ position = (147 -256);+ excludepolicy = "include whole CPT";+ ismapvariable = "false";+ ID = "D";+ label = "D: Coma";+ diagnosistype = "AUXILIARY";+}+node E+{+ states = ("Present" "Absent" );+ position = (414 -266);+ excludepolicy = "include whole CPT";+ ismapvariable = "false";+ ID = "E";+ label = "E: Severe Headaches";+ diagnosistype = "AUXILIARY";+}+node A+{+ states = ("Present" "Absent" );+ position = (131 0);+ excludepolicy = "include whole CPT";+ ismapvariable = "false";+ ID = "A";+ label = "A:Metastatic Cancer";+ diagnosistype = "AUXILIARY";+}+node C+{+ states = ("Present" "Absent" );+ position = (255 -128);+ excludepolicy = "include whole CPT";+ ismapvariable = "false";+ ID = "C";+ label = "C: Brain Tumor";+ diagnosistype = "AUXILIARY";+}+node B+{+ states = ("Increased" "Not increased" );+ position = (0 -128);+ excludepolicy = "include whole CPT";+ ismapvariable = "false";+ ID = "B";+ label = "B: Serum Calcium";+ diagnosistype = "AUXILIARY";+}+potential ( D | C B )+{+ data = ((( 0.8 0.2 )+ ( 0.8 0.2 ))+ (( 0.8 0.2 )+ ( 0.05 0.95 )));+}+potential ( E | C )+{+ data = (( 0.8 0.2 )+ ( 0.6 0.4 ));+}+potential ( A | )+{+ data = ( 0.2 0.8 );+}+potential ( C | A )+{+ data = (( 0.2 0.8 )+ ( 0.05 0.95 ));+}+potential ( B | A )+{+ data = (( 0.8 0.2 )+ ( 0.2 0.8 ));+}
+ hbayes.cabal view
@@ -0,0 +1,89 @@+-- hbayes.cabal auto-generated by cabal init. For additional options,+-- see+-- http://www.haskell.org/cabal/release/cabal-latest/doc/users-guide/authors.html#pkg-descr.+-- The name of the package.+Name: hbayes++-- 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++-- A short (one-line) description of the package.+Synopsis: Inference with Discrete Bayesian Networks++-- 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. ++-- URL for the project homepage or repository.+Homepage: http://www.alpheccar.org++-- The license under which the package is released.+License: BSD3++-- The file containing the license text.+License-file: LICENSE++-- The package author(s).+Author: alpheccar++-- An email address to which users can send suggestions, bug reports,+-- and patches.+Maintainer: misc@alpheccar.org++-- A copyright notice.+Copyright: Copyright (c) 2012, alpheccar ++Category: Math++Build-type: Simple+tested-with: GHC==7.4.1 ++-- Constraint on the version of Cabal needed to build this package.+Cabal-version: >=1.8++data-files: cancer.net+++Library+ -- Modules exported by the library.+ Exposed-modules:+ Bayes+ Bayes.Factor+ Bayes.ImportExport.HuginNet+ Bayes.VariableElimination+ Bayes.FactorElimination+ Bayes.Test+ Bayes.Test.CompareEliminations+ Bayes.Examples+ Bayes.Examples.Tutorial+ other-modules:+ Paths_hbayes+ Bayes.ImportExport.HuginNet.Splitting++ GHC-Options: -O2 -funbox-strict-fields++ + -- Packages needed in order to build this package.+ Build-depends: + base < 5,+ mtl == 2.0.1.0,+ containers == 0.4.2.1,+ array == 0.4.0.0,+ QuickCheck == 2.4.2,+ pretty == 1.1.1.0,+ boxes,+ vector,+ random,+ split,+ parsec,+ filepath,+ directory,+ test-framework-quickcheck2,+ test-framework,+ test-framework-hunit,+ HUnit