diff --git a/Bayes.hs b/Bayes.hs
new file mode 100644
--- /dev/null
+++ b/Bayes.hs
@@ -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
diff --git a/Bayes/Examples.hs b/Bayes/Examples.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Examples.hs
@@ -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 ()
+
diff --git a/Bayes/Examples/Tutorial.hs b/Bayes/Examples/Tutorial.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Examples/Tutorial.hs
@@ -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 ()
diff --git a/Bayes/Factor.hs b/Bayes/Factor.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Factor.hs
@@ -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
+
diff --git a/Bayes/FactorElimination.hs b/Bayes/FactorElimination.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/FactorElimination.hs
@@ -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
diff --git a/Bayes/ImportExport/HuginNet.hs b/Bayes/ImportExport/HuginNet.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/ImportExport/HuginNet.hs
@@ -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)
diff --git a/Bayes/ImportExport/HuginNet/Splitting.hs b/Bayes/ImportExport/HuginNet/Splitting.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/ImportExport/HuginNet/Splitting.hs
@@ -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") 
diff --git a/Bayes/Test.hs b/Bayes/Test.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Test.hs
@@ -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
+            ]
+
+    ]
+
+
diff --git a/Bayes/Test/CompareEliminations.hs b/Bayes/Test/CompareEliminations.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Test/CompareEliminations.hs
@@ -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])
+
diff --git a/Bayes/VariableElimination.hs b/Bayes/VariableElimination.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/VariableElimination.hs
@@ -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
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,4 @@
+import Distribution.Simple
+
+
+main = defaultMain
diff --git a/cancer.net b/cancer.net
new file mode 100644
--- /dev/null
+++ b/cancer.net
@@ -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	));
+}
diff --git a/hbayes.cabal b/hbayes.cabal
new file mode 100644
--- /dev/null
+++ b/hbayes.cabal
@@ -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