diff --git a/Bayes.hs b/Bayes.hs
--- a/Bayes.hs
+++ b/Bayes.hs
@@ -36,6 +36,7 @@
   , edgeEndPoints
   , connectedGraph
   , dag
+  , printGraphValues
   -- * SimpleGraph implementation
   -- ** The SimpleGraph type
   , DirectedSG
@@ -49,6 +50,7 @@
   , evalBN
   , execBN
   , variable
+  , unamedVariable
   , variableWithSize
   , cpt
   , proba
@@ -646,6 +648,10 @@
     tell "\"] ;" 
   tell "\n"
 
+-- | Print the values of the graph vertices
+printGraphValues :: (Graph (SimpleGraph n), Show b) => SimpleGraph n e b -> IO () 
+printGraphValues g@(SP _ _ nm) = putStrLn . execWriter $ mapM_ (printNode nm) (allNodes g)
+
 instance (Show b, Show e) => Show (DirectedSG e b)where
   show g@(SP em vm nm) = execWriter $ do
   tell "digraph dot {\n"
@@ -653,7 +659,6 @@
   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 
@@ -671,7 +676,6 @@
   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 
@@ -834,16 +838,42 @@
 -- | 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)
+  ((namemap,_),_) <- get
+  maybe (getNewEmptyVariable (Just vertexName) initValue) 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)
 
+-- | Generate a new unique unamed empty variable
+getNewEmptyVariable :: NamedGraph g => Maybe String -> f -> GraphMonad g e f Vertex  
+getNewEmptyVariable name initValue = do 
+  ((namemap,count),g) <- get 
+  let vertexName = maybe ("unamed" ++ show count) id name
+      g1 = addLabeledVertex vertexName (Vertex count) initValue g
+      namemap1 = M.insert vertexName count namemap
+  put $! ((namemap1,count+1),g1)
+  return (Vertex count)
+
+-- | Initialize a new variable
+_initializeNewVariable :: (Enum a, Bounded a, NamedGraph g)
+                       => Vertex 
+                       -> a 
+                       -> BNMonad g f DV
+_initializeNewVariable va e = do 
+  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) 
+
+-- | Create a new unamed variable
+unamedVariable :: (Enum a, Bounded a, NamedGraph g)
+               => a -- ^ Variable bounds 
+               -> BNMonad g f DV 
+unamedVariable e = do 
+  va <- getNewEmptyVariable Nothing (UninitializedBNode "unamed" 0)
+  _initializeNewVariable va e
+
 -- | Define a Bayesian variable (name and bounds)
 variable :: (Enum a, Bounded a, NamedGraph g) 
         => String -- ^ Variable name
@@ -851,12 +881,7 @@
         -> 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)
+  _initializeNewVariable va e
 
 -- | Define a Bayesian variable (name and bounds)
 variableWithSize :: NamedGraph g
diff --git a/Bayes/Examples/Tutorial.hs b/Bayes/Examples/Tutorial.hs
--- a/Bayes/Examples/Tutorial.hs
+++ b/Bayes/Examples/Tutorial.hs
@@ -128,9 +128,11 @@
     -- * Tests with the cancer network
     , inferencesOnCancerNetwork
 #endif
+#ifdef LOCAL
+    , miscDiabete
+#endif
     , Coma(..)
     , miscTest
---    , miscDiabete
 	) where 
 
 import Bayes.Factor
@@ -148,11 +150,13 @@
 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
+#ifdef LOCAL
+miscDiabete = do 
+  (varmap,perso) <- exampleDiabete
+  let jtperso = createJunctionTree nodeComparisonForTriangulation perso
+      cho0 = fromJust . Map.lookup "cho_0" $ varmap
+  print $ posterior jtperso cho0
+#endif
 
 miscTest s = do 
   (varmap,perso) <- anyExample s
@@ -202,6 +206,7 @@
 inferencesOnStandardNetwork = do
     let ([winter,sprinkler,rain,wet,road],exampleG) = example
 
+    print exampleG
     putStrLn ""
     print "VARIABLE ELIMINATION"
     putStrLn ""
@@ -221,7 +226,8 @@
     putStrLn ""
 
     let jt = createJunctionTree nodeComparisonForTriangulation exampleG
-
+    print jt
+    displayTreeValues jt
     putStrLn ""
     print "FACTOR ELIMINATION"
     putStrLn ""
diff --git a/Bayes/Factor.hs b/Bayes/Factor.hs
--- a/Bayes/Factor.hs
+++ b/Bayes/Factor.hs
@@ -31,6 +31,7 @@
  , CPT
  -- * Tests
  , testProductProject_prop
+ , testAssocProduct_prop
  , testScale_prop
  , testProjectCommut_prop
  , testScalarProduct_prop
@@ -277,6 +278,10 @@
 
 testScalarProduct_prop :: Double -> CPT -> Bool 
 testScalarProduct_prop v f = (factorProduct [(Scalar v),f]) `isomorphicFactor` (v `factorScale` f)
+
+testAssocProduct_prop :: CPT -> CPT -> CPT -> Bool
+testAssocProduct_prop a b c = (factorProduct [factorProduct [a,b],c] `isomorphicFactor` factorProduct [a,factorProduct [b,c]]) &&
+  (factorProduct [a,b,c] `isomorphicFactor` (factorProduct [factorProduct [a,b],c]) )
 
 testProjectionToScalar_prop :: CPT -> Bool 
 testProjectionToScalar_prop f = 
diff --git a/Bayes/FactorElimination.hs b/Bayes/FactorElimination.hs
--- a/Bayes/FactorElimination.hs
+++ b/Bayes/FactorElimination.hs
@@ -18,6 +18,7 @@
     , createJunctionTree
     , createUninitializedJunctionTree
     , JunctionTree
+    , displayTreeValues
     -- * Shenoy-Shafer message passing
     , collect 
     , distribute
@@ -26,11 +27,13 @@
     , changeEvidence
     -- * Test 
     , junctionTreeProperty_prop
+    , junctionTreeAllClusters_prop
     , VertexCluster
     -- * For debug 
     , junctionTreeProperty
     , maximumSpanningTree
     , fromVertexCluster
+    , triangulatedebug
     ) where
 
 import Bayes
@@ -39,19 +42,20 @@
 import Control.Monad(mapM,guard)
 import Bayes.Factor hiding (isEmpty)
 import Data.Function(on)
-import Data.List(minimumBy,maximumBy,inits,foldl')
+import Data.List(minimumBy,maximumBy,inits,foldl',nub,(\\))
 import qualified Data.Set as Set
 import qualified Data.Map as Map
 import qualified Data.Functor as Functor
 import qualified Data.Tree as T 
 import Bayes.FactorElimination.JTree
 import Control.Applicative((<$>))
+import Bayes.VariableElimination(marginal)
 
 import Test.QuickCheck hiding ((.||.), collect)
 import Test.QuickCheck.Arbitrary
 
 --import Debug.Trace
---debug s a = trace (s ++ " " ++ show a ++ "\n") a
+--debug s a = trace (s ++ "\n" ++ show a ++ "\n") a
 
 {-
  
@@ -63,16 +67,16 @@
 numberOfAddedEdges :: UndirectedGraph g 
                    => g a b 
                    -> Vertex 
-                   -> Int 
+                   -> Integer 
 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)]
+    fromIntegral $ length [edge x y | x <- nodes, y <- nodes, x /= y, not (isLinkedWithAnEdge g x y)]
 
 weightedEdges :: (UndirectedGraph g, Factor f) 
               => g a f 
               -> Vertex 
-              -> Int 
+              -> Integer 
 weightedEdges g v = 
     let nodes = fromJust $ neighbors g v
     in 
@@ -82,9 +86,9 @@
 weight :: (UndirectedGraph g, Factor f)
        => g a f 
        -> Vertex 
-       -> Int 
+       -> Integer 
 weight g v = 
-    factorDimension . fromJust . vertexValue g $ v
+    fromIntegral $ factorDimension . fromJust . vertexValue g $ v
 
 (.||.) :: (a -> a -> Ordering)
        -> (a -> a -> Ordering) 
@@ -109,43 +113,45 @@
 -}
 
 -- | A cluster containing only the vertices and not yet the factors
-newtype VertexCluster = VertexCluster (Set.Set Vertex) deriving(Eq)
+newtype VertexCluster = VertexCluster (Set.Set Vertex) deriving(Eq,Ord)
 
 fromVertexCluster (VertexCluster s) = s
 
 instance Show VertexCluster where 
     show (VertexCluster s) = show . Set.toList $ s
 
-
-
-instance IsCluster Cluster where 
-  overlappingEvidence (Cluster c) e = filter (\x -> Set.member (instantiationVariable x) c) e
-  clusterVariables (Cluster s) = Set.toList s
-  mkSeparator (Cluster sa) (Cluster sb) = Cluster $ Set.intersection sa sb
-
 -- | Triangulate a graph using a cost function
 -- The result is the triangulated graph and the list of clusters
 -- which may not be maximal.
 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)
+            -> [VertexCluster] -- ^ Returns the clusters and the triangulated graph
+triangulate cmp gr = removeNodes cmp gr []
+ where 
+  removeNodes cmp g l | hasNoVertices g = keepMaximalClusters (reverse l)
+                      | otherwise = 
+                          let selectedNode = minimumBy cmp (allVertices g)
+                              theNeighbors = fromJust $ neighbors g selectedNode
+                              g' = removeVertex selectedNode . connectAllNonAdjacentNodes theNeighbors $ g 
+                              newCluster = VertexCluster . Set.fromList $ (selectedNode:theNeighbors)
+                          in 
+                          removeNodes cmp g' (newCluster:l)
 
-    in 
-    processAllNodes g g []
+triangulatedebug :: Graph g
+            => (Vertex -> Vertex -> Ordering) -- ^ Criterion function for triangulation
+            -> g () b
+            -> ([VertexCluster],[g () b]) -- ^ Returns the clusters and the triangulated graph
+triangulatedebug cmp gr = removeNodes cmp gr [] []
+ where 
+  removeNodes cmp g l gl | hasNoVertices g = (reverse l,reverse gl)
+                         | otherwise = 
+                             let selectedNode = minimumBy cmp (allVertices g)
+                                 theNeighbors = fromJust $ neighbors g selectedNode
+                                 g' = removeVertex selectedNode . connectAllNonAdjacentNodes theNeighbors $ g 
+                                 newCluster = VertexCluster . Set.fromList $ (selectedNode:theNeighbors)
+                             in 
+                             removeNodes cmp g' (newCluster:l) (g:gl)
 
 
 -- | Find for a containing cluster. 
@@ -212,23 +218,23 @@
 -- | Get all possible edges between the leaves and the remaining nodes
 possibilities :: (Ord c , UndirectedGraph g) 
               => g Int c -- ^ Original graph to get the edge value 
-              -> JTree c (Vertex,f) -- ^ Tree to get the vertex for a leaf
+              -> JTree c f -- ^ Tree to get the vertex for a leaf
               -> [Vertex] -- ^ Vertices to add to the tree
               -> [c] -- ^ List of leaves
               -> [(Vertex,c,Int)] -- ^ Found edge to add
 possibilities g currentT remaining leavesClusters = do 
   rv <- remaining
   lv <- leavesClusters
-  let NodeValue (lvVertex,lvCluster) _ = nodeValue currentT lv
+  let NodeValue lvVertex lvCluster _ = nodeValue currentT lv
   guard (isLinkedWithAnEdge g rv lvVertex)
   let ev = fromJust $ edgeValue g (edge rv lvVertex)
   return $ (rv,lv,ev)
 
 -- | Find the max edge to add to the tree
-findMax :: (UndirectedGraph g, Ord c, Factor f)
+findMax :: (UndirectedGraph g, Ord c, Factor f,Show c)
         => g Int c -- ^ Graph
         -> [Vertex] -- ^ Nodes to add 
-        -> JTree c (Vertex,f)
+        -> JTree c f
         -> ([Vertex],(Vertex,c),c) 
 findMax g remaining currentT = 
   let leavesClusters = treeNodes currentT
@@ -239,41 +245,29 @@
   in
   (remaining', (rf, foundCluster), lf)
 
-removeVertices :: JTree c (Vertex,f) -> JTree c f
-removeVertices t = t { nodeValueMap = Map.map removeVertexFromNode (nodeValueMap t)
-                     , separatorValueMap = Map.map removeVertexFromSeparator (separatorValueMap t)
-                     }
- where 
-  removeVertexFromNode (NodeValue (_,f) (_,e)) = NodeValue f e 
-  removeVertexFromSeparator (SeparatorValue (_,u) (Just (_,d))) = SeparatorValue u (Just d)
-  removeVertexFromSeparator (SeparatorValue (_,u) Nothing) = SeparatorValue u Nothing 
-  removeVertexFromSeparator EmptySeparator = EmptySeparator
-
 -- | Implementing the Prim's algorithm for minimum spanning tree
-maximumSpanningTree :: (UndirectedGraph g, IsCluster c, Factor f, Ord c) 
+maximumSpanningTree :: (UndirectedGraph g, IsCluster c, Factor f, Ord c, Show c, Show f) 
                     => g Int c 
                     -> JTree c f
 maximumSpanningTree g = 
     let rootNodeVertex = fromJust $ someVertex g 
         rootNodeValue = fromJust $ vertexValue g rootNodeVertex
-        unitFactor = factorFromScalar 1.0 
-        startTree = singletonTree rootNodeValue (rootNodeVertex,unitFactor) (rootNodeVertex,unitFactor) 
+        startTree = singletonTree rootNodeValue rootNodeVertex [] [] 
         remainingVertices = filter (/= rootNodeVertex) (allVertices g) 
     in 
-    removeVertices $ buildTree g remainingVertices startTree 
+    buildTree g remainingVertices startTree 
 
-buildTree :: (UndirectedGraph g , IsCluster c, Factor f, Ord c)
+buildTree :: (UndirectedGraph g , IsCluster c, Factor f, Ord c, Show c, Show f)
           => g Int c 
           -> [Vertex]
-          -> JTree c (Vertex,f) 
-          -> JTree c (Vertex,f) 
+          -> JTree c f 
+          -> JTree c f
 buildTree g [] currentT = currentT 
-buildTree g l@(h:t) currentT = 
-    let unitFactor = factorFromScalar 1.0
-        (l',(foundElemVertex,foundElemValue),leaf) = findMax g l currentT
+buildTree g l currentT = 
+    let (l',(foundElemVertex,foundElemValue),leaf) = findMax g l currentT
         sep = mkSeparator foundElemValue leaf
         newTree = addSeparator leaf sep foundElemValue . 
-                  addNode foundElemValue (foundElemVertex,unitFactor) (foundElemVertex,unitFactor) $ currentT
+                  addNode foundElemValue foundElemVertex [] [] $ currentT
     in 
     buildTree g l' newTree
    
@@ -285,13 +279,13 @@
 
 
 -- | Create a junction tree with only the clusters and no factors
-createUninitializedJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f)
+createUninitializedJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f, Show f)
                                 => (UndirectedSG () f -> Vertex -> Vertex -> Ordering) -- ^ Weight function on the moral graph
                                 -> g () f -- ^ Input directed graph
                                 -> JunctionTree f -- ^ Junction tree
 createUninitializedJunctionTree cmp g =  
   let theMoralGraph = moralGraph g
-      (clusters,_) = triangulate (cmp theMoralGraph) theMoralGraph
+      clusters = triangulate (cmp theMoralGraph) theMoralGraph
       g'' = createClusterGraph g clusters :: UndirectedSG Int Cluster
   in 
   maximumSpanningTree g''
@@ -315,10 +309,12 @@
 posterior t v = 
   case snd $ traverseTree (findClusterFor v) Nothing t of 
     Nothing -> Nothing
-    Just c -> let NodeValue f e = nodeValue t c 
+    Just c -> let NodeValue ver f e = nodeValue t c 
                   d = maybe (factorFromScalar 1.0) id $ downMessage t =<< (nodeParent t c)
                   u = map (upMessage t) (nodeChildren t c)
-                  unNormalized = factorProjectTo [v] (factorProduct (f:e:d:u))
+                  allFactors = d:u ++ f ++ e
+                  variablesToRemove = (nub (concatMap factorVariables allFactors)) \\ [v]
+                  unNormalized = marginal allFactors variablesToRemove [v] []
               in 
               Just $ factorDivide unNormalized (factorNorm unNormalized)
 
@@ -341,6 +337,19 @@
   in
   junctionTreeProperty t [] (root t)
 
+junctionTreeAllClusters_prop :: DirectedSG () CPT -> Property 
+junctionTreeAllClusters_prop g = (not . isEmpty) g && (not . hasNoEdges) g && connectedGraph g ==> 
+      let theMoralGraph = moralGraph g
+          cmp ug = (compare `on` (numberOfAddedEdges ug))
+          clusters = triangulate (cmp theMoralGraph) theMoralGraph
+          g'' = createClusterGraph g clusters :: UndirectedSG Int Cluster
+          jt = maximumSpanningTree g'' :: JunctionTree CPT
+          treeClusters = treeNodes jt 
+          sa = Set.fromList (map (vertexClusterToCluster g) clusters) 
+          sb = Set.fromList treeClusters 
+      in 
+      Set.isSubsetOf sa sb && Set.isSubsetOf sb sa
+
 junctionTreeProperty :: JTree Cluster CPT -> [Cluster] -> Cluster -> Bool
 junctionTreeProperty t path c = 
   let cl = map (separatorChild t) . nodeChildren t $ c
@@ -374,26 +383,18 @@
 
 -- | 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)]
+connectAllNonAdjacentNodes :: (Graph g) 
+                           => [Vertex]  -- ^ List of nodes to connect
+                           -> g () b -- ^ Graph containing the nodes
+                           -> g () b
+connectAllNonAdjacentNodes nodes originGraph   =  
+    let addEmptyEdge g e = addEdge e () g
     in 
-    (originGraph',dstGraph')
-
+    foldl' addEmptyEdge originGraph [edge x y | x <- nodes, y <- nodes, x /= y, not (isLinkedWithAnEdge originGraph x y)]
+   
 -- | Add the missing parent links
 addMissingLinks :: DirectedGraph g => g () b -> Vertex -> b -> g () b
-addMissingLinks g v _ = 
-    let (_,g') = connectAllNodesWith g g (\e m -> m) (\e m -> addEdge e () m) (parents g v)
-    in 
-    g'
+addMissingLinks g v _ = connectAllNonAdjacentNodes (parents g v) g 
 
 
 -- | Convert the graph to an undirected form
diff --git a/Bayes/FactorElimination/JTree.hs b/Bayes/FactorElimination/JTree.hs
--- a/Bayes/FactorElimination/JTree.hs
+++ b/Bayes/FactorElimination/JTree.hs
@@ -7,11 +7,13 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
 module Bayes.FactorElimination.JTree(
       IsCluster(..)
     , Cluster(..)
     , JTree(..)
     , JunctionTree(..)
+    , Sep
     , setFactors
     , distribute 
     , collect
@@ -32,6 +34,8 @@
     , traverseTree
     , separatorChild
     , treeNodes
+    , treeValues
+    , displayTreeValues
     , Action(..)
     ) where 
 
@@ -40,10 +44,12 @@
 import Data.Maybe(fromJust,mapMaybe)
 import qualified Data.Set as Set
 import Data.Monoid
-import Data.List((\\), intersect,partition, foldl')
+import Data.List((\\), intersect,partition, foldl',minimumBy,nub)
 import Bayes.PrivateTypes 
 import Bayes.Factor
 import Bayes
+import Data.Function(on)
+import Bayes.VariableElimination(marginal)
 
 import Debug.Trace 
 debug s a = trace (s ++ " " ++ show a ++ "\n") a
@@ -61,40 +67,57 @@
     show (SeparatorValue u Nothing) = "u(" ++ show u ++ ")"
     show (SeparatorValue u (Just d)) = "u(" ++ show u ++ ") d(" ++ show d ++ ")"
 
-type FactorValue a = a 
-type EvidenceValue a = a
+type FactorValues a = [a]
+type EvidenceValues a = [a]
 
 -- | Node value
-data NodeValue a = NodeValue !(FactorValue a) !(EvidenceValue a) deriving(Eq)
+data NodeValue a = NodeValue !Vertex !(FactorValues a) !(EvidenceValues a) deriving(Eq)
 
 instance Show a => Show (NodeValue a) where 
-    show (NodeValue f e) = "f(" ++ show f ++ ") e(" ++ show e ++ ")"
+    show (NodeValue v f e) = "f(" ++ show f ++ ") e(" ++ show e ++ ")"
 
+newtype Sep = Sep Int deriving(Eq,Ord,Show,Num)
+
 -- | Junction tree.
 -- 'c' is the node / separator identifier (for instance a set of 'DV')
 -- a are the values for a node or separator
+-- Cluster are unique sor the cluster value is also the cluster key
+-- Separator values are not unique. Two different seperators can be the same
+-- cluster. So, separator unicity is enforced with a number
 data JTree  c f = JTree {  root :: !c
                         -- | Leaves of the tree
                         ,  leavesSet :: !(Set.Set c)
                         -- | The children of a node are separators
-                        ,  childrenMap :: !(Map.Map c [c])
+                        ,  childrenMap :: !(Map.Map c [Sep])
                         -- | Parent of a node
-                        ,  parentMap :: !(Map.Map c c)
+                        ,  parentMap :: !(Map.Map c Sep)
                         -- | Parent of a separator
-                        ,  separatorParentMap :: !(Map.Map c c)
+                        ,  separatorParentMap :: !(Map.Map Sep c)
                         -- | The child of a seperator is a node
-                        ,  separatorChildMap :: !(Map.Map c c)
+                        ,  separatorChildMap :: !(Map.Map Sep c)
                         -- | Values for nodes and seperators
                         ,  nodeValueMap :: !(Map.Map c (NodeValue f))
-                        ,  separatorValueMap :: !(Map.Map c (SeparatorValue f))
+                        ,  separatorValueMap :: !(Map.Map Sep (SeparatorValue f))
+                        ,  separatorCurrentKey :: !Sep
+                        ,  separatorClusterMap :: !(Map.Map Sep c)
                         } deriving(Eq)
 
 -- | Create a singleton tree with just one root node
-singletonTree r factorValue evidenceValue = 
-    let t = JTree r Set.empty Map.empty Map.empty Map.empty Map.empty Map.empty Map.empty
+singletonTree r rootVertex factorValue evidenceValue = 
+    let t = JTree r Set.empty Map.empty Map.empty Map.empty Map.empty Map.empty Map.empty (Sep 0) Map.empty
     in 
-    addNode r factorValue evidenceValue t
+    addNode r rootVertex factorValue evidenceValue t
 
+-- | Reset all evidences to 1 in the network
+resetEvidences :: Factor f => JTree c f -> JTree c f 
+resetEvidences t = t {nodeValueMap = Map.map resetNodeEvidence (nodeValueMap t)}
+ where 
+  resetNodeEvidence (NodeValue v f _) = NodeValue v f []
+
+-- | Get the cluster for a separator
+separatorCluster :: JTree c a -> Sep -> c 
+separatorCluster t s = fromJust $ Map.lookup s (separatorClusterMap t)
+
 -- | Leaves of the tree
 leaves :: JTree c a -> [c]
 leaves = Set.toList . leavesSet
@@ -103,6 +126,9 @@
 treeNodes :: JTree c a -> [c]
 treeNodes = Map.keys . nodeValueMap
 
+treeValues :: JTree c f -> [(c,NodeValue f)]
+treeValues = Map.toList . nodeValueMap
+
 -- | Value of a node
 nodeValue :: Ord c => JTree c a -> c -> NodeValue a 
 nodeValue t e = fromJust $ Map.lookup e (nodeValueMap t)
@@ -112,36 +138,36 @@
 setNodeValue c v t = t {nodeValueMap = Map.insert c v (nodeValueMap t)} 
 
 -- | Parent of a node
-nodeParent :: Ord c => JTree c a -> c -> Maybe c 
+nodeParent :: Ord c => JTree c a -> c -> Maybe Sep 
 nodeParent t e = Map.lookup e (parentMap t)
 
 -- | Value of a node
-separatorValue :: Ord c => JTree c a -> c -> SeparatorValue a 
+separatorValue :: Ord c => JTree c a -> Sep -> SeparatorValue a 
 separatorValue t e = fromJust $ Map.lookup e (separatorValueMap t)
 
 -- | Parent of a separator
-separatorParent :: Ord c => JTree c a -> c -> c 
+separatorParent :: Ord c => JTree c a -> Sep -> c 
 separatorParent t e = fromJust $ Map.lookup e (separatorParentMap t)
 
 -- | UpMessage for a separator node
-upMessage :: Ord c => JTree c a -> c -> a
+upMessage :: Ord c => JTree c a -> Sep -> a
 upMessage t c = case separatorValue t c of 
                   SeparatorValue up _ -> up 
                   _ -> error "Trying to get an up message on an empty seperator ! Should never occur !"
 
 -- | DownMessage for a separator node
-downMessage :: Ord c => JTree c a -> c -> Maybe a 
+downMessage :: Ord c => JTree c a -> Sep -> Maybe a 
 downMessage t c = case separatorValue t c of 
      SeparatorValue _ (Just down) -> Just down 
      SeparatorValue _ Nothing -> Nothing
      _ -> error "Trying to get a down message on an empty separator ! Should never occur !"
 
 -- | Return the separator childrens of a node
-nodeChildren :: Ord c => JTree c a -> c -> [c]
+nodeChildren :: Ord c => JTree c a -> c -> [Sep]
 nodeChildren t e = maybe [] id $ Map.lookup e (childrenMap t)
 
 -- | Return the child of a separator
-separatorChild :: Ord c => JTree c a -> c -> c 
+separatorChild :: Ord c => JTree c a -> Sep -> c 
 separatorChild t e = fromJust $ Map.lookup e (separatorChildMap t)
 
 -- | Check if a node is member of the tree
@@ -152,34 +178,39 @@
 -- The nodes MUST already be in the tree
 addSeparator :: (Ord c) 
              => c -- ^ Origin node 
-             -> c -- ^ Separator
+             -> c -- ^ Separator value
              -> c -- ^ Destination node 
              -> JTree c a -- ^ Current tree 
              -> JTree c a -- ^ Modified tree 
-addSeparator node sep dest t = 
-    t { childrenMap = Map.insertWith' (++) node [sep] (childrenMap t)
-      , separatorChildMap = Map.insert sep dest (separatorChildMap t)
-      , separatorValueMap = Map.insert sep EmptySeparator (separatorValueMap t)
+addSeparator node sepCluster dest t = 
+    let newSep = (separatorCurrentKey t) + 1 
+    in
+    t { childrenMap = Map.insertWith' (++) node [newSep] (childrenMap t)
+      , separatorChildMap = Map.insert newSep dest (separatorChildMap t)
+      , separatorValueMap = Map.insert newSep EmptySeparator (separatorValueMap t)
+      , separatorClusterMap = Map.insert newSep sepCluster (separatorClusterMap t)
       , leavesSet = Set.delete node (leavesSet t) 
-      , parentMap = Map.insert dest sep (parentMap t)
-      , separatorParentMap = Map.insert sep node (separatorParentMap t)
+      , parentMap = Map.insert dest newSep (parentMap t)
+      , separatorParentMap = Map.insert newSep node (separatorParentMap t)
+      , separatorCurrentKey = newSep
       }
 
 -- | Add a new node
 addNode :: (Ord c) 
         => c -- ^ Node
-        -> a -- ^ Factor value 
-        -> a -- ^ Evidence value
+        -> Vertex
+        -> [a] -- ^ Factor value 
+        -> [a] -- ^ Evidence value
         -> JTree c a 
         -> JTree c a 
-addNode node factorValue evidenceValue t = 
-   t { nodeValueMap = Map.insert node (NodeValue factorValue evidenceValue) (nodeValueMap t)
+addNode node vertex factorValue evidenceValue t = 
+   t { nodeValueMap = Map.insert node (NodeValue vertex factorValue evidenceValue) (nodeValueMap t)
      , leavesSet = Set.insert node (leavesSet t) 
      }
 
 -- | Update the up message of a separator
 updateUpMessage :: Ord c 
-                => Maybe c -- ^ Separator node to update (if any : none for root node)
+                => Maybe Sep -- ^ Separator node to update (if any : none for root node)
                 -> a -- ^ New value
                 -> JTree c a -- ^ Old tree
                 -> JTree c a
@@ -193,7 +224,7 @@
 
 -- | Update the down message of a separator
 updateDownMessage :: Ord c 
-                  => c -- ^ Separator node to update
+                  => Sep -- ^ Separator node to update
                   -> a -- ^ New value
                   -> JTree c a -- ^ Old tree
                   -> JTree c a
@@ -223,7 +254,7 @@
 
 allSeparatorsHaveReceivedAMessage :: Ord c
                                   => JTree c a -- ^ Tree
-                                  -> [c] -- ^ Separators
+                                  -> [Sep] -- ^ Separators
                                   -> Bool 
 allSeparatorsHaveReceivedAMessage t seps = 
   all separatorInitialized . map (separatorValue t) $ seps
@@ -245,7 +276,8 @@
               in 
               case destinationNode of 
                 Nothing -> t -- When root
-                Just p -> let generatedMessage = newMessage incomingMessages currentValue p
+                Just p -> let  sepC = separatorCluster t p
+                               generatedMessage = newMessage incomingMessages currentValue sepC
                           in 
                           updateUpMessage destinationNode generatedMessage t
 
@@ -253,56 +285,63 @@
 updateDownSeparator :: (Message a c, Ord c) 
                     => c -- ^ Node generating the message 
                     -> JTree c a 
-                    -> c -- ^ Child receiving the message
+                    -> Sep -- ^ Child receiving the message
                     -> JTree c a 
 updateDownSeparator node t child  = 
     let incomingMessagesFromBelow = map (upMessage t) (nodeChildren t node \\ [child])
         messageFromAbove = downMessage t =<< (nodeParent t node)
         incomingMessages = maybe incomingMessagesFromBelow (\x -> x:incomingMessagesFromBelow) messageFromAbove
         currentValue = nodeValue t node
-        generatedMessage = newMessage incomingMessages currentValue child
+        childC = separatorCluster t child
+        generatedMessage = newMessage incomingMessages currentValue childC
     in 
     updateDownMessage child generatedMessage t
 
 unique :: Ord c => [c] -> [c]
 unique = Set.toList . Set.fromList
 
-data TraversalState = ACluster | ASeparator
-
 -- | Collect message taking into account that the tree depth may be different for different leaves.
 collect :: (Ord c, Message a c) 
         => JTree c a 
         -> JTree c a
-collect t = _collect ACluster (leaves t) t
+collect t = _collectNodes (leaves t) t
 
-_collect :: (Ord c, Message a c) 
-         => TraversalState -- ^ Node processing phase or separator processing phase
-         -> [c]
-         -> JTree c a -- ^ Tree
-         -> JTree c a -- ^ Modified tree 
-_collect _ [] t = t
-_collect ACluster l t = 
+_collectSeparators :: (Ord c, Message a c) 
+                   => [Sep]
+                   -> JTree c a -- ^ Tree
+                   -> JTree c a -- ^ Modified tree
+_collectSeparators l t = _collectNodes (unique . map (separatorParent t) $ l) t
+
+_collectNodes :: (Ord c, Message a c) 
+              => [c]
+              -> JTree c a -- ^ Tree
+              -> JTree c a -- ^ Modified tree 
+_collectNodes  [] t = t
+_collectNodes  l t = 
     let newTree = foldl' updateUpSeparator t l
     in 
-    _collect ASeparator (mapMaybe (nodeParent t) l) newTree
-_collect ASeparator l t = _collect ACluster (unique . map (separatorParent t) $ l) t
+    _collectSeparators (mapMaybe (nodeParent t) l) newTree
     
 distribute :: (Ord c, Message a c)
            => JTree c a 
            -> JTree c a
-distribute t = _distribute ACluster t (root t) 
+distribute t = _distributeNodes t (root t) 
 
-_distribute :: (Ord c, Message a c)
-            => TraversalState -- ^ True if node
-            -> JTree c a 
-            -> c -- ^ Destination of the distribute
-            -> JTree c a 
-_distribute ACluster t node = 
+_distributeSeparators :: (Ord c, Message a c)
+                      => JTree c a 
+                      -> Sep -- ^ Destination of the distribute
+                      -> JTree c a 
+_distributeSeparators t node = _distributeNodes t (separatorChild t node)
+
+_distributeNodes :: (Ord c, Message a c)
+                 => JTree c a 
+                 -> c -- ^ Destination of the distribute
+                 -> JTree c a 
+_distributeNodes t node = 
     let children = nodeChildren t node
         newTree = foldl' (updateDownSeparator node) t $ children
     in
-    foldl' (_distribute ASeparator) newTree children
-_distribute ASeparator t node = _distribute ACluster t (separatorChild t node)
+    foldl' _distributeSeparators newTree children
 
 {-
 
@@ -313,8 +352,11 @@
 -- | This class is used to check if evidence or a factor is relevant
 -- for a cluster
 class IsCluster c where 
+  -- | Evidence contained in the cluster
   overlappingEvidence :: c -> [DVI Int] -> [DVI Int]
+  -- | Cluser variables
   clusterVariables :: c -> [DV]
+  -- | Intersection of two clusters
   mkSeparator :: c -> c -> c
 
 instance IsCluster [DV] where 
@@ -333,18 +375,19 @@
              -> s -- ^ Current state
              -> JTree c f -- ^ Input tree
              -> (JTree c f,s)
-traverseTree action state t = _traverseTree True action (t,state) (root t)
+traverseTree action state t = _traverseTreeNodes action (t,state) (root t)
 
-_traverseTree False action (t,state) current  = _traverseTree True action (t,state)  (separatorChild t current) 
-_traverseTree True action (t,state) current = 
+_traverseTreeSeparators action (t,state) current = _traverseTreeNodes  action (t,state)  (separatorChild t current) 
+
+_traverseTreeNodes action (t,state) current = 
   case action state current (nodeValue t current) of 
      Stop newState -> (t,newState)
      ModifyAndStop _ newValue -> (setNodeValue current newValue t, state) 
-     Skip newState -> foldl' (_traverseTree False action) (t,newState) (nodeChildren t current)
+     Skip newState -> foldl' (_traverseTreeSeparators action) (t,newState) (nodeChildren t current)
      Modify newState newValue -> 
          let newTree = setNodeValue current newValue t 
          in 
-         foldl' (_traverseTree False action) (newTree,newState) (nodeChildren newTree current)
+         foldl' (_traverseTreeSeparators action) (newTree,newState) (nodeChildren newTree current)
 
 mapWithCluster :: Ord c 
                => (c -> NodeValue f -> NodeValue f)
@@ -353,52 +396,49 @@
 mapWithCluster f t = t {nodeValueMap = Map.mapWithKey f (nodeValueMap t)}
 
 -- | Set the factors in the tree 
-setFactors :: (Graph g, Factor f, Show f, IsCluster c, Ord c)
+updateTreeValues :: (Factor f, IsCluster c, Ord c, Show c, Show f)
+                 => (f -> NodeValue f -> NodeValue f) 
+                 -> [f]
+                 -> JTree c f 
+                 -> JTree c f
+updateTreeValues change factors t = 
+  let allNodes = treeNodes t
+      factorIncludedInCluster f c = all (`elem` clusterVariables c) (factorVariables f)
+      coveringClusters f = filter (f `factorIncludedInCluster`) allNodes
+      clusterSize a = product . map (fromIntegral . dimension) . clusterVariables $ a :: Integer
+      addFactor t newFactor = 
+        let minimumCluster = minimumBy (compare `on` clusterSize) (coveringClusters newFactor)
+            clusterValue = nodeValue t minimumCluster
+        in 
+        setNodeValue minimumCluster (change newFactor clusterValue) t
+  in 
+  foldl' addFactor t factors
+
+-- | Set the factors in the tree 
+setFactors :: (Graph g, Factor f, IsCluster c, Ord c, Show c, Show f)
            => BayesianNetwork g f 
            -> JTree c f 
            -> JTree c f
 setFactors g t = 
   let factors = allVertexValues g 
-  in
-  fst . traverseTree updateFactor factors $ t 
-
-
--- | Update factors in a cluster
-updateFactor :: (Factor f, IsCluster c) 
-             => [f] -- ^ Remaining list of factors to attribute
-             -> c -- ^ Current cluster
-             -> NodeValue f -- ^ Current value
-             -> Action [f] (NodeValue f)
-updateFactor lf c (NodeValue _ evidence) | null lf = Stop lf
-                                         | otherwise =
-  let isFactorIncluded l = all (`elem` clusterVariables c) (factorVariables l)
-      (attributedFactors,remainingFactors) = partition isFactorIncluded lf 
+      changeFactor f (NodeValue v oldf e) = NodeValue v (f:oldf) e
   in 
-  Modify remainingFactors  (NodeValue (factorProduct attributedFactors) evidence)
+  updateTreeValues  changeFactor factors t  
 
 -- | Change evidence in the network
-changeEvidence :: (IsCluster c, Ord c, Factor f, Message f c)
+changeEvidence :: (IsCluster c, Ord c, Factor f, Message f c, Show c, Show f)
                => [DVI Int] -- ^ Evidence
                -> JTree c f 
                -> JTree c f 
-changeEvidence e t =  distribute . 
-                      collect . fst .
-                      traverseTree changeNodeEvidence e $ 
-                      t { separatorValueMap = Map.map (const EmptySeparator) (separatorValueMap t)}
-
-changeNodeEvidence :: (IsCluster c, Factor f) 
-                   => [DVI Int] -- ^ Evidence
-                   -> c  -- ^ Current cluster
-                   -> NodeValue f -- ^ Current value
-                   -> Action [DVI Int] (NodeValue f)
-changeNodeEvidence [] c v = Stop []
-changeNodeEvidence e c (NodeValue f olde) = 
-  let oe = overlappingEvidence c e
-      ns = e \\ oe
-      newEvidence = factorProduct $ map factorFromInstantiation oe
+changeEvidence e t =  
+  let evidences = map factorFromInstantiation e
+      changeEvidence newe (NodeValue v f olde) = NodeValue v f (newe:olde)
   in 
-  Modify ns (NodeValue f newEvidence)
-
+  distribute . 
+  collect .
+  updateTreeValues changeEvidence evidences .
+  resetEvidences $
+  t { separatorValueMap = Map.map (const EmptySeparator) (separatorValueMap t)}                
 
 -- | Cluster of discrete variables.
 -- Discrete variables instead of vertices are needed because the
@@ -406,6 +446,11 @@
 -- which factors must be contained in a given cluster.
 newtype Cluster = Cluster (Set.Set DV) deriving(Eq,Ord)
 
+instance IsCluster Cluster where 
+  overlappingEvidence c = overlappingEvidence (fromCluster c)
+  clusterVariables c = clusterVariables (fromCluster c)
+  mkSeparator (Cluster a) (Cluster b) = Cluster (Set.intersection a b)
+
 instance Show Cluster where 
   show (Cluster s) = show . Set.toList $ s
 
@@ -413,7 +458,12 @@
 
 
 instance Factor f => Message f Cluster where 
-  newMessage input (NodeValue f e) dv = factorProjectTo (fromCluster dv) (factorProduct (f:e:input))
+  newMessage input (NodeValue _ f e) dv = 
+    let allFactors = f ++ e ++ input 
+        variablesToKeep = fromCluster dv 
+        variablesToRemove = (nub (concatMap factorVariables allFactors)) \\ variablesToKeep
+    in 
+    marginal allFactors variablesToRemove variablesToKeep []
 
 
 type JunctionTree f = JTree Cluster f
@@ -439,48 +489,53 @@
 toTree d t = 
     let r = root t
         v = nodeValue t r
-        nodec = map S (nodeChildren t r)
+        nodec = nodeChildren t r
     in 
-    Tree.Node (label d (show r) v) (_toTree d t nodec)
+    Tree.Node (label d (show r) v) (_toTreeSeparators d t nodec)
 
 
-_toTree :: (Ord c, Show c, Show a) 
-        => Bool
-        -> JTree c a 
-        -> [NodeKind c] 
-        -> [Tree.Tree String]
-_toTree _ _ [] = []
-_toTree d t ((N h):l) = 
-    let nodec = map S (nodeChildren t h) -- Node children are separators
+_toTreeNodes :: (Ord c, Show c, Show a) 
+             => Bool
+             -> JTree c a 
+             -> [c] 
+             -> [Tree.Tree String]
+_toTreeNodes _ _ [] = []
+_toTreeNodes d t (h:l) = 
+    let nodec = nodeChildren t h -- Node children are separators
         v = nodeValue t h
     in
-    Tree.Node (label d (show h) v) (_toTree d t nodec):_toTree d t l
-_toTree d t ((S h):l) = 
-    let separatorc = [N $ separatorChild t h] -- separator child is a node
+    Tree.Node (label d (show h) v) (_toTreeSeparators d t nodec):_toTreeNodes d t l
+
+_toTreeSeparators :: (Ord c, Show c, Show a) 
+                  => Bool
+                  -> JTree c a 
+                  -> [Sep] 
+                  -> [Tree.Tree String]
+_toTreeSeparators _ _ [] = []    
+_toTreeSeparators d t (h:l) = 
+    let separatorc = [separatorChild t h] -- separator child is a node
         v = separatorValue t h
     in
-    Tree.Node (label d ("<" ++ show h ++ ">") v ) (_toTree d t separatorc):_toTree d t l
+    Tree.Node (label d ("<" ++ show (separatorCluster t h) ++ ">") v ) (_toTreeNodes d t separatorc):_toTreeSeparators d t l
 
 instance (Ord c, Show c, Show a) => Show (JTree c a) where 
     show = Tree.drawTree . toTree False
 
 displayTree b = Tree.drawTree . toTree b
 
-{-
+-- | Display the tree values
+displayTreeValues :: (Show f, Show c) => JTree c f -> IO ()
+displayTreeValues t = 
+  let allValues = treeValues t
+      printAValue (c,NodeValue _ f e) = do 
+        print c 
+        putStrLn "FACTOR"
+        print f 
+        putStrLn "EVIDENCE"
+        print e 
+        putStrLn "------"
 
-Debug functions for tests
+  in 
+  mapM_ printAValue allValues
 
--}
 
---instance Message (Sum Int) String where 
---    newMessage l (NodeValue a b) _ = mconcat (a:b:l)
---
---testTree :: JTree String (Sum Int)
---testTree = let s a= Sum a
---           in
---           addSeparator "ROOT" "RB" "B" .
---           addNode "B" (s 3) (s 3) . 
---           addSeparator "ROOT" "RA" "A"  . 
---           addNode "A" (s 2) (s 2) $ 
---           singletonTree "ROOT" (s 4) (s 5)
---
diff --git a/Bayes/Test.hs b/Bayes/Test.hs
--- a/Bayes/Test.hs
+++ b/Bayes/Test.hs
@@ -10,8 +10,8 @@
 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)
+import Bayes.Factor(testProductProject_prop,testScale_prop,testProjectCommut_prop,testScalarProduct_prop,testProjectionToScalar_prop,testAssocProduct_prop)
+import Bayes.FactorElimination(junctionTreeProperty_prop,junctionTreeAllClusters_prop)
 
 #ifdef LOCAL
 import Bayes.Test.ReferencePatterns(compareAsiaReference,compareCancerReference,comparePokerReference,compareFarmReference)
@@ -30,11 +30,13 @@
                 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
+                testProperty "Test projection to scalar" testProjectionToScalar_prop,
+                testProperty "Test associativity of factor" testAssocProduct_prop
             ]
         , testGroup "Junction Tree" [
                 testProperty "Test the junction tree property" junctionTreeProperty_prop,
-                testCase "Test variable elimination == factor elimination" compareVariableFactor
+                testCase "Test variable elimination == factor elimination" compareVariableFactor,
+                testProperty "Test all clusters are included in the junction tree" junctionTreeAllClusters_prop
             ]
 #ifdef LOCAL
         , testGroup "Reference patterns" [ 
diff --git a/Bayes/VariableElimination.hs b/Bayes/VariableElimination.hs
--- a/Bayes/VariableElimination.hs
+++ b/Bayes/VariableElimination.hs
@@ -11,12 +11,13 @@
  , minDegreeOrder
  , minFillOrder
  , allVariables
+ , marginal
  , EliminationOrder
  ) where
 
 import Bayes
 import Bayes.Factor
-import Data.List(partition,minimumBy,(\\),find)
+import Data.List(partition,minimumBy,(\\),find,foldl')
 import Data.Maybe(fromJust)
 import Data.Function(on)
 import qualified Data.Map as M
@@ -39,87 +40,103 @@
   map createDV s
 
 -- | Used for bucket elimination. Factor are organized by their first DV
-type Buckets f = (EliminationOrder,M.Map DV [f])
+data Buckets f = Buckets !EliminationOrder !(M.Map DV [f])
 
-createBuckets ::  (Graph g, Factor f, Show f) 
-              => BayesianNetwork g f -- ^ Bayesian Network
+createBuckets ::  (Factor f) 
+              => [f] -- ^ Factor to use for computing the marginal one
               -> 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
+createBuckets s e r = 
+  let -- 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) =
+      addDVToBucket (rf, m) dv  =
         let (fk,remaining) = partition (flip containsVariable dv) rf
         in 
         (remaining, M.insert dv fk m)
-      (_,b) = foldr addDVToBucket (s,M.empty) (reverse theOrder)
+      (_,b) = foldl' addDVToBucket (s,M.empty) theOrder
   in
-  (tail theOrder,b)
+  Buckets (tail theOrder) b
 
 -- | Get the factors for a bucket
 getBucket :: DV 
           -> Buckets f 
           -> [f]
-getBucket dv (_,m) = fromJust $ M.lookup dv m
+getBucket dv (Buckets _ m) = fromJust $ M.lookup dv m
 
 -- | Update bucket
-updateBucket :: Factor f => DV -> f -> Buckets f -> Buckets f 
-updateBucket dv f b@(e,m) = 
+updateBucket :: Factor f
+             => DV -- ^ Variable that was eliminated
+             -> f -- ^ New factor resulting from this elimination
+             -> Buckets f 
+             -> Buckets f 
+updateBucket dv f b@(Buckets e m) = 
   if isScalarFactor f 
     then 
-      (tail e,M.insert dv [f] m)
+      Buckets (tail e) (M.insert dv [f] m)
     else
       let b' = removeFromBucket dv b
-          (e',m') = addBucket f b'
+          Buckets e' m' = addBucket b' f 
       in 
-      (tail e',m')
+      Buckets (tail e') m'
 
 -- | Add a factor to the right bucket
-addBucket :: Factor f => f -> Buckets f -> Buckets f
-addBucket f (e,b) = 
+addBucket :: Factor f => Buckets f -> f -> Buckets f
+addBucket (Buckets e b) f = 
   let inBucket = find (f `containsVariable`) e
   in 
   case inBucket of 
-    Nothing -> (e,b)
-    Just bucket -> (e, M.insertWith' (++) bucket [f] b)
+    Nothing -> Buckets e b
+    Just bucket -> Buckets 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) 
+removeFromBucket dv (Buckets e m) = Buckets 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 = 
+marginal :: (Factor f) 
+         => [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
+marginal lf 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
+  let bucket = createBuckets lf p r
       assignmentFactors = map factorFromInstantiation assignment
-      bucket' = foldr addBucket bucket assignmentFactors
-      (_,resultBucket) = foldr marginalizeOneVariable bucket' (reverse p)
+      bucket' = foldl' addBucket bucket assignmentFactors
+      Buckets _ resultBucket = foldl' marginalizeOneVariable bucket' 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 
+  -- We get P(Q , e)
+  resultFactor
  where 
-  marginalizeOneVariable dv currentBucket = 
+  marginalizeOneVariable currentBucket dv   = 
     let fk = getBucket dv currentBucket
         p = factorProduct fk
         f' = factorProjectOut [dv] p
     in
     updateBucket dv f' currentBucket
 
+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 g p r assignment = 
+  let s = allVertexValues g 
+      resultFactor = marginal s p r assignment
+      norm = factorNorm resultFactor
+  in
+  -- We get P(Q | e)
+  resultFactor `factorDivide` norm 
+
 -- | Compute the prior marginal. All the variables in the
 -- elimination order are conditionning variables ( p( . | conditionning variables) )
 priorMarginal :: (Graph g, Factor f, Show f) 
@@ -127,7 +144,13 @@
               -> EliminationOrder -- ^ Ordering of variables to marginalize
               -> EliminationOrder -- ^ Ordering of remaining to keep in result
               -> f
-priorMarginal g ea eb = posteriorMarginal g ea eb []
+priorMarginal g ea eb = 
+  let s = allVertexValues g 
+      resultFactor = marginal s ea eb []
+      norm = factorNorm resultFactor
+  in
+  -- We get P(Q | e)
+  resultFactor `factorDivide` norm    
 
 -- | Compute the interaction graph of the BayesianNetwork
 interactionGraph :: (FoldableWithVertex g,Factor f, UndirectedGraph g')
@@ -139,12 +162,12 @@
   addFactor vertex factor graph = 
     let allvars = factorVariables factor
         edges = [(x,y) | x <- allvars, y <- allvars , x /= y]
-        addNewEdge (va,vb) g = 
+        addNewEdge g (va,vb)  = 
           let g' = addVertex (variableVertex vb) vb . addVertex (variableVertex va) va $ g 
           in
           addEdge (edge (variableVertex va) (variableVertex vb)) () $ g'
     in 
-    foldr addNewEdge graph edges
+    foldl' addNewEdge graph edges
 
 -- | Number of neighbors for a variable in the bayesian network
 nbNeighbors :: UndirectedSG () DV 
@@ -172,16 +195,16 @@
             -> Int 
 degreeOrder g p =
   let  ig = interactionGraph g :: UndirectedSG () DV
-       (_,w) = foldr processVariable (ig,0) p 
+       (_,w) = foldl' processVariable (ig,0) p 
   in 
   w 
  where 
-  addAnEdge (va,vb) g = addEdge (edge va vb) () g
-  processVariable bdv (g,w) = 
+  addAnEdge g (va,vb)  = addEdge (edge va vb) () g
+  processVariable (g,w) bdv  = 
     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)
+        g' = removeVertex (variableVertex bdv) (foldl' addAnEdge g edges)
     in
     if nbNeighbors > w 
       then 
diff --git a/hbayes.cabal b/hbayes.cabal
--- a/hbayes.cabal
+++ b/hbayes.cabal
@@ -7,7 +7,7 @@
 -- The package version. See the Haskell package versioning policy
 -- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for
 -- standards guiding when and how versions should be incremented.
-Version:             0.2
+Version:             0.2.1
 
 -- A short (one-line) description of the package.
 Synopsis:            Inference with Discrete Bayesian Networks
@@ -15,7 +15,7 @@
 -- A longer description of the package.
 Description:  Algorithms for inference with Discrete Bayesian Networks.  
  It is a very preliminary version. It has only been tested on very simple
- examples where it worked. This 0.2 version is using new faster and cleaner algorithms.
+ examples where it worked. This 0.2.1 version is using new faster and cleaner algorithms and correcting a problem with graph triangulation.
 
 -- URL for the project homepage or repository.
 Homepage:            http://www.alpheccar.org