diff --git a/Bayes.hs b/Bayes.hs
--- a/Bayes.hs
+++ b/Bayes.hs
@@ -35,6 +35,7 @@
   , newEdge
   , edgeEndPoints
   , connectedGraph
+  , dag
   -- * SimpleGraph implementation
   -- ** The SimpleGraph type
   , DirectedSG
@@ -63,7 +64,7 @@
 import Control.Monad.State.Strict
 import Control.Monad.Writer.Strict
 import Control.Applicative((<$>))
-import Bayes.Factor
+import Bayes.Factor hiding(isEmpty)
 import Data.Maybe
 import qualified Data.Map as Map
 import qualified Data.Foldable as F
@@ -74,6 +75,7 @@
 import Test.QuickCheck
 import Test.QuickCheck.Arbitrary
 import Data.List(sort,intercalate,nub)
+import Bayes.PrivateTypes hiding(isEmpty)
 
 --import Debug.Trace
 --debug a = trace (show a) a
@@ -120,6 +122,29 @@
 
     foldM createEdge g edges   
 
+-- | Warning : the generated graph is not at all a bayesian network
+-- The variables in the CPT have no reason to correspond to the edges
+-- connected to that CPT.
+-- Only the main variable (first variable) is linked to the right vertex
+instance Factor f => Arbitrary (DirectedSG () f) where
+  arbitrary = do 
+    let createVertex g i = do 
+          let value = fromJust $ factorWithVariables [DV (Vertex i) 2] [0.1,0.9]
+          return $ addVertex (Vertex i) value g
+        createEdge g (va,vb) = do 
+          return $ addEdge (edge va vb) () g 
+
+    nbVertex <- choose (1,8) :: Gen Int
+    
+    g <- foldM createVertex emptyGraph [1..nbVertex]
+
+    let allPairs = [(Vertex x,Vertex y) | x <- [1..nbVertex], y <- [1..nbVertex], x /= y]
+        anEdge (x,y) = arbitrary :: Gen Bool
+
+    edges <- filterM anEdge allPairs
+
+    foldM createEdge g edges
+
 testEdgeRemoval_prop :: DirectedSG String String -> Property
 testEdgeRemoval_prop g = (not . hasNoEdges) g ==> 
   let Just e = someEdge g
@@ -248,7 +273,26 @@
     ingoing :: g a b -> Vertex -> Maybe [Edge]
     outgoing :: g a b -> Vertex -> Maybe [Edge]
 
+-- | Get the root node for the graph
+rootNode :: DirectedGraph g => g a b -> Maybe Vertex
+rootNode g = 
+  let someRoots = filter (isRoot g) . allVertices $ g
+  in 
+  case someRoots of 
+    (h:l) -> Just h 
+    _ -> Nothing
+  where 
+    isRoot g v =
+      case ingoing g v of 
+        Just [] -> True 
+        _ -> False
 
+-- | Check if the graph is a directed Acyclic graph
+dag :: DirectedGraph g => g a b -> Bool 
+dag g = case rootNode g of 
+  Nothing -> isEmpty g 
+  Just r -> dag (removeVertex r g)
+
 -- | Check if the graph is connected
 connectedGraph :: Graph g => g a b -> Bool 
 connectedGraph g = 
@@ -384,12 +428,14 @@
 class FoldableWithVertex g where
   -- | Fold with vertex 
   foldrWithVertex :: (Vertex -> a -> b -> b) -> b -> g c a -> b 
+  foldlWithVertex' :: (b -> Vertex -> a -> b) -> b -> g c a -> b 
 
 instance FoldableWithVertex (SimpleGraph local) where
   foldrWithVertex f s (SP _ vm _) = IM.foldrWithKey (\k (_,v) y -> f (Vertex k) v y) s vm
+  foldlWithVertex' f s (SP _ vm _) = IM.foldlWithKey' (\y k (_,v)  -> f y (Vertex k) v) s vm
 
 _addLabeledVertex vertexName vert@(Vertex v) value (SP em vm name) =
-  let vm' = IM.insertWith noRedundancy v (emptyNeighborhood,value) vm
+  let vm' = IM.insertWith' noRedundancy v (emptyNeighborhood,value) vm
       name' = IM.insert v vertexName name 
   in
   SP em vm' name'
@@ -498,7 +544,7 @@
     else 
       Just . head . M.keys $ em
 
-_addVertex vert@(Vertex v) value (SP em vm nm) = SP em (IM.insertWith noRedundancy v (emptyNeighborhood,value) vm) nm
+_addVertex vert@(Vertex v) value (SP em vm nm) = SP em (IM.insertWith' noRedundancy v (emptyNeighborhood,value) vm) nm
 
 _removeVertex v@(Vertex vertex) g@(SP _ vm _)  = maybe g removeVertexWithValue (IM.lookup vertex vm)
   where
@@ -587,7 +633,7 @@
   tell "\n"
   let r = IM.lookup k nm
   when (isJust r) $ do
-     tell $ fromJust r
+     tell $ "Node " ++ fromJust r
   tell "\n"
   tell $ show v
   tell "\n"
diff --git a/Bayes/Examples.hs b/Bayes/Examples.hs
--- a/Bayes/Examples.hs
+++ b/Bayes/Examples.hs
@@ -69,7 +69,9 @@
 module Bayes.Examples(
    example
  , exampleJunction
+#ifndef LOCAL
  , exampleImport
+#endif
  , exampleDiabete
  , exampleAsia
  , examplePoker
@@ -86,6 +88,8 @@
 import qualified Data.Map as Map
 import System.Directory(getHomeDirectory)
 import System.FilePath((</>))
+
+#ifndef LOCAL
 import Paths_hbayes
 
 -- | Example showing how to import a graph described into
@@ -95,6 +99,7 @@
     path <- getDataFileName "cancer.net"
     r <- importBayesianGraph path
     return (runBN $ fromJust r)
+#endif
 
 -- | Genereic loading functions to load some other
 -- examples from the author's dropbox.
@@ -136,7 +141,7 @@
 
 
 -- | Standard example found in many books about Bayesian Networks.
-example :: (DVSet,SBN CPT)
+example :: ([DV],SBN CPT)
 example = runBN $ do 
     winter <- variable "winter" (t :: Bool)
     sprinkler <- variable "sprinkler" (t :: Bool) 
diff --git a/Bayes/Examples/Tutorial.hs b/Bayes/Examples/Tutorial.hs
--- a/Bayes/Examples/Tutorial.hs
+++ b/Bayes/Examples/Tutorial.hs
@@ -124,16 +124,23 @@
 module Bayes.Examples.Tutorial(
     -- * Tests with the standard network 
       inferencesOnStandardNetwork
+#ifndef LOCAL
     -- * Tests with the cancer network
     , inferencesOnCancerNetwork
+#endif
     , Coma(..)
     , miscTest
+--    , miscDiabete
 	) where 
 
 import Bayes.Factor
 import Bayes
 import Bayes.VariableElimination
+#ifndef LOCAL
 import Bayes.Examples(example, exampleJunction,exampleImport,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)
+#else 
+import Bayes.Examples(example, exampleJunction,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)
+#endif
 import Bayes.FactorElimination
 import Data.Function(on)
 import qualified Data.Map as Map
@@ -141,11 +148,11 @@
 import System.Exit(exitSuccess)
 import qualified Data.List as L((\\))
 
-miscDiabete = do 
-  (varmap,perso) <- exampleDiabete
-  let jtperso = createJunctionTree nodeComparisonForTriangulation perso
-      cho0 = fromJust . Map.lookup "cho_0" $ varmap
-  print $ posterior jtperso cho0
+--miscDiabete = do 
+--  (varmap,perso) <- exampleDiabete
+--  let jtperso = createJunctionTree nodeComparisonForTriangulation perso
+--      cho0 = fromJust . Map.lookup "cho_0" $ varmap
+--  print $ posterior jtperso cho0
 
 miscTest s = do 
   (varmap,perso) <- anyExample s
@@ -171,6 +178,7 @@
 -- from the cancer network.
 data Coma = Present | Absent deriving(Eq,Enum,Bounded)
 
+#ifndef LOCAL
 -- | Inferences with the cancer network
 inferencesOnCancerNetwork = do 
     print "CANCER NETWORK"
@@ -182,12 +190,13 @@
     mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer $ x)) [varA,varB,varC,varD,varE]
 
     print "UPDATED EVIDENCE : Coma present"
-    let jtcancer' = updateEvidence [varD =: Present] jtcancer 
+    let jtcancer' = changeEvidence [varD =: Present] jtcancer 
     mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varD,varE]
 
     print "UPDATED EVIDENCE : Coma absent"
-    let jtcancer' = updateEvidence [varD =: Absent] jtcancer 
+    let jtcancer' = changeEvidence [varD =: Absent] jtcancer 
     mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varD,varE]
+#endif 
 
 -- | Inferences with the standard network
 inferencesOnStandardNetwork = do
@@ -221,20 +230,20 @@
     print m
     putStrLn ""
 
-    let jt' = updateEvidence [wet =: True] jt 
+    let jt' = changeEvidence [wet =: True] jt 
 
     print "Posterior Marginal : probability of rain if grass wet"
     let m = posterior jt' rain
     print m
     putStrLn ""
 
-    let jt'' = clearEvidence jt'
+    let jt'' = changeEvidence [] jt'
     print "Prior Marginal : probability of rain"
     let m = posterior jt rain
     print m
     putStrLn ""
 
-    let jt3 = updateEvidence [wet =: True, sprinkler =: True] jt'
+    let jt3 = changeEvidence [wet =: True, sprinkler =: True] jt'
 
     print "Posterior Marginal : probability of rain if grass wet and sprinkler used"
     let m = posterior jt3 rain
diff --git a/Bayes/Factor.hs b/Bayes/Factor.hs
--- a/Bayes/Factor.hs
+++ b/Bayes/Factor.hs
@@ -16,10 +16,9 @@
  -- * Implementation
  , Vertex(..)
  -- ** Discrete variables and instantiations
- , DV(..)
- , DVSet(..)
+ , DV
+ --, DVSet(..)
  , DVI
- , DVISet(..)
  , setDVValue
  , instantiationValue
  , instantiationVariable
@@ -40,7 +39,7 @@
 
 import qualified Data.Vector.Unboxed as V
 import Data.Vector.Unboxed((!))
-import Data.Maybe(fromJust,mapMaybe)
+import Data.Maybe(fromJust,mapMaybe,isJust)
 import qualified Data.List as L
 import Text.PrettyPrint.Boxes hiding((//))
 import Test.QuickCheck
@@ -48,79 +47,19 @@
 import qualified Data.IntMap as IM
 import Control.Monad
 import System.Random(Random)
+import Data.List(partition)
+import Bayes.PrivateTypes
 
 --import Debug.Trace
 
 --debug a = trace ("\nDEBUG\n" ++ show a ++ "\n") a
 
--- | Vertex type used to identify a vertex in a graph
-newtype Vertex = Vertex {vertexId :: Int} deriving(Eq,Ord)
 
-instance Show Vertex where 
-    show (Vertex v) = "v" ++ show v
-
--- | A Set of variables used in a factor. s is the set and a the variable
-class Set s where
-    -- | Empty set
-    emptySet :: s a
-    -- | Union of two sets
-    union :: Eq a => s a -> s a -> s a
-    -- | Intersection of two sets
-    intersection :: Eq a => s a -> s a -> s a
-    -- | Difference of two sets
-    difference :: Eq a => s a -> s a -> s a
-    -- | Check if the set is empty
-    isEmpty :: s a -> Bool
-    -- | Check if an element is member of the set
-    isElem :: Eq a => a -> s a -> Bool
-    -- | Add an element to the set
-    addElem :: Eq a => a -> s a -> s a
-    -- | Number of elements in the set
-    nbElements :: s a -> Int
-
-    -- | Check if a set is subset of another one
-    subset :: Eq a => s a -> s a -> Bool
-
-    -- | Check set equality
-    equal :: Eq a => s a -> s a -> Bool
-    equal sa sb = (sa `subset` sb) && (sb `subset` sa)
-
-instance Set [] where
-    emptySet = []
-    union = L.union
-    intersection = L.intersect
-    difference a b = a L.\\ b
-    isEmpty [] = True 
-    isEmpty _ = False
-    isElem = L.elem 
-    addElem a l = if a `elem` l then l else a:l
-    nbElements = length
-    subset sa sb = all (`elem` sb) sa
-
--- | A discrete variable has a number of levels which is required to size the factors
-class BayesianDiscreteVariable v where
-    dimension :: v -> Int 
-
-
 -- | A vertex associated to another value (variable dimension, variable value ...)
 class LabeledVertex l where
     variableVertex :: l -> Vertex
 
--- | A discrete variable
-data DV = DV !Vertex !Int deriving(Eq,Ord)
 
--- | A set of discrete variables
-type DVSet = [DV]
-
-instance Show DV where
-    show (DV v d) = show v ++ "(" ++ show d ++ ")"
-
--- | Discrete Variable instantiation. A variable and its value
-data DVI a = DVI DV !a deriving(Eq)
-
-instance Show a => Show (DVI a) where 
-   show (DVI (DV v _) i) = show v ++ "=" ++ show i
-
 -- | Convert a variable instantation to a factor
 -- Useful to create evidence factors
 factorFromInstantiation :: Factor f => DVI Int -> f
@@ -129,88 +68,16 @@
     in
     fromJust . factorWithVariables [dv] . map (setValue) $ [0..dimension dv-1]
 
--- | A set of variable instantiations
-type DVISet a = [DVI a]
 
-instance BayesianDiscreteVariable DV where
-    dimension (DV _ d) = d
 
--- | Create a discrete variable instantiation for a given discrete variable
-setDVValue :: DV -> a -> DVI a
-setDVValue v a = DVI v a
 
-getMinBound :: Bounded a => a -> a 
-getMinBound _ = minBound
-
--- | Create a variable instantiation using values from
--- an enumeration
-(=:) :: (Bounded b, Enum b) => DV -> b -> DVI Int 
-(=:) a b = setDVValue a (fromEnum b - fromEnum (getMinBound b))
-
--- | Extract value of the instantiation
-instantiationValue (DVI _ v) = v
-
--- | Discrete variable from the instantiation
-instantiationVariable (DVI dv _) = dv
-
 instance LabeledVertex (DVI a) where
     variableVertex (DVI v _) = variableVertex v
 
 instance LabeledVertex DV where
     variableVertex (DV v _) = v
 
--- | Extend indexing to full variable set using a bool
--- list and a default value
--- For instance [True, False, True, False] 5 [2,3] ---> [2,5,3,5]
-extend :: [Bool] -> a -> [a] -> [a]
-extend [] _ l = l
-extend (h:t) d [] = d:extend t d []
-extend (False:t) d l = d:extend t d l
-extend (True:t) d (h:l') = h:extend t d l'
 
--- | Inner loop function using full indices for full variables
-type InnerLoop a = [Int] -> a
-
--- | Outer loop function using result from inner loop
--- and outer vars indices
-type OuterLoop a b = [Int] -> [a] -> b
-
--- | Iter on outer var and inner var
--- Inner body is called with indiced for full set
--- Outer body is called with indices for outer set
-forSubA :: DVSet -- ^ All variables
-        -> DVSet -- ^ Outer variables
-        -> (DVSet -> [Int] -> [a]) -- ^ Inner loop body
-        -> OuterLoop a b -- ^ Outer loop function
-        -> [b]
-forSubA allvars outervars inner outer = 
-    let sCode s e = if (e `isElem` s) then True else False
-        selection s = map (sCode s) allvars
-        computeOuter iouter =
-            let outerIdx =  extend (selection outervars) 0 iouter
-                innerValues = inner allvars outerIdx
-            in 
-            outer iouter innerValues
-    in
-    map computeOuter (forAllIndices outervars)
-
--- | Use indices with full variable set
-forSubB :: DVSet -- ^ Inner vars 
-        -> InnerLoop a -- ^ Inner loop function
-        -> DVSet -- ^ All vars
-        -> [Int] -- ^ Outer indices
-        -> [a]
-forSubB innervars f allvars  outerIdx  = 
-        let sCode s e = if (e `isElem` s) then True else False
-            selection s = map (sCode s) allvars
-            computeInner iinner =
-                let innerIdx = extend (selection innervars) 0 iinner
-                    idx = zipWith (+) outerIdx innerIdx
-                    in 
-                    f idx
-        in
-        map computeInner (forAllIndices innervars)
-
 -- | Norm the factor
 normedFactor :: Factor f => f -> f 
 normedFactor f = factorDivide f (factorNorm f)
@@ -218,7 +85,7 @@
 -- | A factor as used in graphical model
 -- It may or not be a probability distribution. So it has no reason to be
 -- normalized to 1
-class FactorPrivate f => Factor f where
+class Factor f where
     -- | When all variables of a factor have been summed out, we have a scalar
     isScalarFactor :: f -> Bool 
     -- | An empty factor with no variable and no values
@@ -226,20 +93,24 @@
     -- | Check if a given discrete variable is contained in a factor
     containsVariable :: f -> DV  -> Bool
     -- | Give the set of discrete variables used by the factor
-    factorVariables :: f -> DVSet    
+    factorVariables :: f -> [DV]    
     -- | Return A in P(A | C D ...). It is making sense only if the factor is a conditional propbability
     -- table. It must always be in the vertex corresponding to A in the bayesian graph
     factorMainVariable :: f -> DV
-    factorMainVariable = head . factorVariables
+    factorMainVariable f = let vars = factorVariables f 
+      in
+      case vars of 
+        [] -> error "Can't get the main variable of a scalar factor"
+        (h:_) -> h 
     -- | Create a new factors with given set of variables and a list of value
     -- for initialization. The creation may fail if the number of values is not
     -- coherent with the variables and their levels.
     -- For boolean variables ABC, the value must be given in order
     -- FFF, FFT, FTF, FTT ...
-    factorWithVariables :: DVSet -> [Double] -> Maybe f
+    factorWithVariables :: [DV] -> [Double] -> Maybe f
     -- | Value of factor for a given set of variable instantitation.
     -- The variable instantion is like a multi-dimensional index.
-    factorValue :: f -> DVISet Int -> Double
+    factorValue :: f -> [DVI Int] -> Double
     -- | Position of a discrete variable in te factor (p(AB) is differennt from p(BA) since values
     -- are not organized in same order in memory)
     variablePosition :: f -> DV -> Maybe Int
@@ -258,107 +129,42 @@
 
     -- | Create an evidence factor from an instantiation.
     -- If the instantiation is empty then we get nothing
-    evidenceFrom :: DVISet Int -> Maybe f
+    evidenceFrom :: [DVI Int] -> Maybe f
     
 
     -- | Divide all the factor values
     factorDivide :: f -> Double -> f
     factorDivide f d = (1.0 / d) `factorScale` f 
 
+    factorToList :: f -> [Double]
+
     -- | Multiply factors. 
     factorProduct :: [f] -> f
-    factorProduct [] = factorFromScalar 1.0
-    factorProduct l = 
-        let allVars = L.foldl1' union . map factorVariables $ l
-        in 
-        if L.null allVars 
-            then 
-                factorFromScalar (product . map factorNorm $ l) 
-            else
-                let getFactorValueAtIndex i factor = factorValuePrivate factor (reorder i factor)
-                    instantiationProduct instantiation = product . map (getFactorValueAtIndex instantiation) $ l
-                    values = [instantiationProduct x | x <- forAllInstantiations allVars]
-                in 
-                fromJust $ factorWithVariables allVars values
 
     -- | Project out a factor. The variable in the DVSet are summed out
-    factorProjectOut :: DVSet -> f -> f
-    factorProjectOut s f = 
-        let alls = factorVariables f
-            s' = alls `difference` s
-        in 
-        if null s'
-            then 
-                factorFromScalar (factorNorm f)
-            else
-                let dstValues = forSubA alls s' 
-                                   (forSubB s $ factorValuePrivate f)
-                                   (\i c -> sum c)
-                in 
-                fromJust $ factorWithVariables s' dstValues
+    factorProjectOut :: [DV] -> f -> f
+
     -- | Project to. The variable are kept and other variables are removed
-    factorProjectTo :: DVSet -> f -> f 
+    factorProjectTo :: [DV] -> f -> f 
     factorProjectTo s f = 
         let alls = factorVariables f 
             s' = alls `difference` s 
         in 
         factorProjectOut s' f
 
--- | Used internaly when we know the position of a variable in the factor
--- then we can identify the variable with an int. May be a bit faster for some
--- algorithms
-class FactorPrivate f where
-    factorValuePrivate :: f -> [Int] -> Double
-
--- | Return all the index (position in the factor) for a DV
-allValues :: DV -> [Int]
-allValues (DV _ i) = [0..i-1]
-
--- | Generate all indexes for a set of variables
-forAllIndices :: DVSet -> [[Int]]
-forAllIndices = mapM allValues
-
--- | Generate all instantiations of variables
-forAllInstantiations :: DVSet -> [DVISet Int]
-forAllInstantiations = mapM oneInstantiation
- where
-    oneInstantiation v@(DV vertex _) = map (setDVValue v) . allValues $ v
-
 -- | Change the layout of values in the
 -- factor to correspond to a new variable order
-changeVariableOrder :: DVSet -- ^ Old order
-                    -> DVSet -- ^ New order 
+-- Used to import external files
+changeVariableOrder :: DVSet s -- ^ Old order
+                    -> DVSet s' -- ^ New order 
                     -> [Double] -- ^ Old values
                     -> [Double] -- ^ New values
-changeVariableOrder oldOrder newOrder oldValues =
+changeVariableOrder (DVSet oldOrder) newOrder oldValues =
     let oldFactor = fromJust $ factorWithVariables oldOrder oldValues :: CPT
     in
     [factorValue oldFactor i | i <- forAllInstantiations newOrder]
 
 
--- | Order the variable to get a multiindex which is
--- making sense in the CPT. Only the subset in CPT is selectionned and reordered
-reorder :: Factor f => DVISet Int -> f  -> [Int]
-reorder i f = 
-    let nbDestVars = nbElements . factorVariables $ f
-        v = V.replicate nbDestVars 0
-        asDV v = DV v 0
-        vectorPair bdvi = do 
-            pos <- variablePosition f . asDV . variableVertex $ bdvi
-            let value = instantiationValue bdvi
-            return (pos, value)
-        allPos = mapMaybe vectorPair i
-    in
-    let testError = maybe False (const True) $ do 
-        guard $ length allPos == nbDestVars
-        guard $ and . map ( (< nbDestVars) . fst)  $ allPos
-        return ()
-    in
-    case testError of
-      False -> error ("reorder has not set all destination indexes ! allpos = " ++ show allPos ++ " nbDestVars = " ++ show nbDestVars ++ "\n" ) 
-      True -> V.toList $ v V.// allPos
-
-
 -- | Mainly used for conditional probability table like p(A B | C D E) but the normalization to 1
 -- is not imposed. And the conditionned variables are not different from the conditionning ones.
 -- The dimensions for each variables are listed.
@@ -368,11 +174,11 @@
 -- the knowledge of the dependencies is.
 -- So, this same structure is used for a probability too (conditioned on nothing)
 data CPT = CPT {
-           dimensions :: DVSet -- ^ Dimensions for all variables
-         , mapping :: IM.IntMap Int -- ^ Mapping from vertex number to position in dimensions
-         , values :: V.Vector Double -- ^ Table of values
+           dimensions :: ![DV] -- ^ Dimensions for all variables
+         , mapping :: !(IM.IntMap Int) -- ^ Mapping from vertex number to position in dimensions
+         , values :: !(V.Vector Double) -- ^ Table of values
          }
-         | Scalar Double
+         | Scalar !Double
 
 debugCPT (Scalar d) = do 
    putStrLn "SCALAR CPT"
@@ -496,9 +302,11 @@
 -- meaning of the variables according to their position.
 isomorphicFactor :: Factor f => f -> f -> Bool
 isomorphicFactor fa fb = maybe False (const True) $ do 
-    let va = factorVariables fa 
-        vb = factorVariables fb 
-    guard (va `equal` vb)
+    let sa = factorVariables fa 
+        sb = factorVariables fb 
+        va = DVSet sa 
+        vb = DVSet sb
+    guard (sa `equal` sb)
     guard (factorDimension fa == factorDimension fb)
     guard $ and [factorValue fa ia `nearlyEqual` factorValue fb ia | ia <- forAllInstantiations va]
     return ()
@@ -508,22 +316,25 @@
 Following functions are used to typeset the factor when displaying it
 
 -}
-vname :: Int -> Int -> Box
-vname vc i = text $ "v" ++ show vc ++ "=" ++ show i
+-- | Display a variable name and its size
+vname :: Int -> DVI Int -> Box
+vname vc i = text $ "v" ++ show vc ++ "=" ++ show (instantiationValue i)
 
-dispFactor :: FactorPrivate f => f -> DV -> [Int] -> DVSet -> Box
+dispFactor :: Factor f => f -> DV -> [DVI Int] -> [DV] -> Box
 dispFactor cpt h c [] = 
-    let dstIndexes = allValues h
+    let dstIndexes = allInstantiationsForOneVariable h
         dependentIndexes =  reverse c
         factorValueAtPosition p = 
-            let v = factorValuePrivate cpt p
+            let v = factorValue cpt p
             in
             text . show  $ v
     in
     vsep 0 center1 . map (factorValueAtPosition . (:dependentIndexes)) $ dstIndexes
 
 dispFactor cpt dst c (h@(DV (Vertex vc) i):l) = 
-    hsep 1 top . map (\i -> vcat center1 [vname vc i,dispFactor cpt dst (i:c) l])  $ (allValues h)
+    let allInst = allInstantiationsForOneVariable h
+    in
+    hsep 1 top . map (\i -> vcat center1 [vname vc i,dispFactor cpt dst (i:c) l])  $ allInst
 
 instance Show CPT where
     show (Scalar v) = "\nScalar Factor:\n" ++ show v
@@ -532,11 +343,17 @@
     show c@(CPT d _ v) = 
         let h@(DV (Vertex vc) _) = head d
             table = dispFactor c h [] (tail d)
-            dstColumn = vcat center1 $ replicate (length d - 1) (text "") ++ map (vname vc) (allValues h)
+            dstIndexes = map head (forAllInstantiations . DVSet $ [h])
+            -- In P(A | B ...), the dst column is containing the possible values for the
+            -- variables A with a header made of space to be aligned with the other part of the table.
+            -- In the other part of the table, this header is containing the variable values for the other varibles
+            dstColumn = vcat center1 $ replicate (length d - 1) (text "") ++ map (vname vc) dstIndexes
         in
         "\n" ++ show d ++ "\n" ++ render (hsep 1 top [dstColumn,table])
 
 instance Factor CPT where
+    factorToList (Scalar v) = [v]
+    factorToList (CPT _ _ v) = V.toList v
     emptyFactor = emptyCPT
     isScalarFactor (Scalar _) = True
     isScalarFactor _ = False
@@ -548,55 +365,111 @@
     factorWithVariables = createCPTWithDims
     factorVariables (CPT v _ _) = v
     factorVariables (Scalar _) = []
-    factorNorm f@(CPT _ _ _) = sum [ factorValuePrivate f x | x <- forAllIndices (factorVariables f)]
+    factorNorm f@(CPT d _ vals) = 
+        let vars = DVSet d
+            strides = indexStrides vars
+        in
+        sum [ vals!(indexPosition strides x) | x <- indicesForDomain vars]
     factorNorm (Scalar v) = v
     variablePosition (CPT _ m _) (DV (Vertex i) _) = IM.lookup i m
     variablePosition (Scalar _) _ = Nothing
     factorScale s (Scalar v) = Scalar (s*v)
-    factorScale s f = 
-        let newValues = map (s *) [ factorValuePrivate f x | x <- forAllIndices (factorVariables f)]
+    factorScale s f@(CPT d _ vals) = 
+        let vars = DVSet d
+            strides = indexStrides vars
+            newValues = map (s *) [ vals!(indexPosition strides x) | x <- indicesForDomain vars]
         in 
         fromJust $ factorWithVariables (factorVariables f) newValues
     factorValue (Scalar v) _ = v 
-    factorValue f i = 
-        let multiIndex = reorder i f
+    factorValue f@(CPT d _ v) i = 
+        let vars = DVSet d
+            (dv,pos) = instantiationDetails i
+            strides = indexStridesFor vars dv
         in 
-        factorValuePrivate f multiIndex
+        v!(indexPosition strides pos)
     evidenceFrom [] = Nothing 
     evidenceFrom l = 
-        let index = map instantiationValue l 
-            variables = map instantiationVariable l
+        let (variables,index) = instantiationDetails l
+            DVSet nakedVars = variables
             setValueForIndex i = if i == index then 1.0 else 0.0 
         in
-        factorWithVariables variables . map setValueForIndex $ forAllIndices variables
-
-instance FactorPrivate CPT where
-    factorValuePrivate = getCPTValue
+        factorWithVariables nakedVars . map setValueForIndex $ indicesForDomain variables
+    factorProduct [] = factorFromScalar 1.0
+    factorProduct l = 
+        let allVars = DVSet $ L.foldl1' union . map factorVariables $ l
+            DVSet nakedVars = allVars
+            (scalars,cpts) = partition isScalarFactor l
+            stridesFromCPT (CPT d _ _) = indexStridesFor (DVSet d) allVars
+            ps = product . map (flip factorValue undefined) $ scalars
+            cptsStrides = map stridesFromCPT cpts
+        in 
+        if L.null cpts 
+            then 
+                factorFromScalar ps
+            else
+                let getFactorValueAtIndex i (strides,factor@(CPT _ _ vals)) = vals!(indexPosition strides i)
+                    instantiationProduct instantiation = product . map (getFactorValueAtIndex instantiation) $ (zip cptsStrides cpts)
+                    values = [ps * instantiationProduct x | x <- indicesForDomain allVars]
+                in 
+                values `seq` fromJust $ factorWithVariables nakedVars values
+    factorProjectOut _ f@(Scalar v) = f
+    factorProjectOut s f@(CPT d _ v) = 
+        let variablesToSum = s
+            variablesToKeep = d `difference` s 
+            keepSet = DVSet variablesToKeep
+            sumSet = DVSet variablesToSum 
+            strides = indexStridesFor (DVSet d) (DVSet $ variablesToKeep ++ variablesToSum)
 
+            values = do 
+                  keepIndex <- indicesForDomain keepSet 
+                  let l = do
+                        sumIndex <- indicesForDomain sumSet 
+                        return $ v!(indexPosition strides $ combineIndex strides keepIndex sumIndex)
+                  return (sum l)
+        in
+        values `seq` fromJust $ factorWithVariables variablesToKeep values
+        
+-- | Used to combined the keep and sum indices in the factor projection
+combineIndex :: Strides s'' -> [Index s] -> [Index s'] -> [Index s''] 
+combineIndex _ la lb = map (Index . fromIndex) la ++ map (Index .fromIndex) lb
 
 -- | An empty CPT
 emptyCPT :: CPT
 emptyCPT = CPT [] IM.empty V.empty
 
--- | Convertion of a multiindex to its
--- position inside of the data vector of a 'CPT'
-indexPosition :: DVSet -> [Int] -> Int
-indexPosition [] _ = 0
-indexPosition d pos = 
-    let dim = map dimension d
+newtype Strides s = Strides [Int] deriving(Eq,Show)
+
+-- | Generate strides to read the first CPT using an index having meaning in the second CPT
+indexStridesFor :: DVSet s -- ^ DVSet to be read
+                -> DVSet s' -- ^ DVSet to interpret the index
+                -> Strides s'
+indexStridesFor dr@(DVSet drvars) di@(DVSet divars) =
+    let Strides originStrides = indexStrides dr
+        reference = zip drvars originStrides 
+        getNewStrides dv = maybe 0 id (lookup dv reference)
+    in 
+    Strides $ map getNewStrides divars
+    
+
+-- | Generate the strides to read a given factor using a multiindex
+-- using the same order as the factor variables
+indexStrides :: DVSet s -> Strides s
+indexStrides d@(DVSet dvars)  = 
+    let dim = map dimension dvars
         pos' = scanr (*) (1::Int) (tail dim)
-        c = sum . map (\(x,y) -> x * y) $ (zip pos' pos)
     in 
-    c
+    Strides pos'
+-- | Convertion of a multiindex to its
+-- position inside of the data vector of a 'CPT'
+indexPosition :: Strides s -> [Index s] -> Int
+{-# INLINE indexPosition #-}
+indexPositions _ []  = 0
+indexPosition (Strides pos') pos = sum . map (\(x,y) -> x * fromIndex y) $ (zip pos' pos)
 
--- | Get the value at a given position. Positions are starting at zero
-getCPTValue :: CPT -> [Int] -> Double
-getCPTValue (Scalar v) _ = v
-getCPTValue cpt@(CPT d _ v) pos = v!(indexPosition d pos)
 
 -- | Create a CPT given some dimensions and a list of Doubles.
 -- Returns nothing is the length are not coherents.
-createCPTWithDims :: DVSet -> [Double] -> Maybe CPT
+createCPTWithDims :: [DV] -> [Double] -> Maybe CPT
 createCPTWithDims dims values = 
     let createDVIndex i (DV (Vertex v) _)  = (v,i)
         m = IM.fromList . zipWith createDVIndex ([0,1..]::[Int]) $ dims
diff --git a/Bayes/FactorElimination.hs b/Bayes/FactorElimination.hs
--- a/Bayes/FactorElimination.hs
+++ b/Bayes/FactorElimination.hs
@@ -9,37 +9,43 @@
     -- * Triangulation
     , nodeComparisonForTriangulation
     , numberOfAddedEdges
+    , weight
+    , weightedEdges
     , triangulate
     -- * Junction tree
-    , minimumSpanningTree
     , createClusterGraph
     , Cluster
     , createJunctionTree
+    , createUninitializedJunctionTree
     , JunctionTree
     -- * Shenoy-Shafer message passing
     , collect 
     , distribute
     , posterior 
     -- * Evidence
-    , clearEvidence
-    , updateEvidence
+    , changeEvidence
     -- * Test 
     , junctionTreeProperty_prop
-    , createVerticesJunctionTree
     , VertexCluster
+    -- * For debug 
+    , junctionTreeProperty
+    , maximumSpanningTree
+    , fromVertexCluster
     ) where
 
 import Bayes
 import qualified Data.Foldable as F
 import Data.Maybe(fromJust,mapMaybe,isJust)
-import Control.Monad(mapM)
+import Control.Monad(mapM,guard)
 import Bayes.Factor hiding (isEmpty)
 import Data.Function(on)
-import Data.List(minimumBy,maximumBy,inits)
+import Data.List(minimumBy,maximumBy,inits,foldl')
 import qualified Data.Set as Set
 import qualified Data.Map as Map
 import qualified Data.Functor as Functor
 import qualified Data.Tree as T 
+import Bayes.FactorElimination.JTree
+import Control.Applicative((<$>))
 
 import Test.QuickCheck hiding ((.||.), collect)
 import Test.QuickCheck.Arbitrary
@@ -63,6 +69,15 @@
     in 
     length [edge x y | x <- nodes, y <- nodes, x /= y, not (isLinkedWithAnEdge g x y)]
 
+weightedEdges :: (UndirectedGraph g, Factor f) 
+              => g a f 
+              -> Vertex 
+              -> Int 
+weightedEdges g v = 
+    let nodes = fromJust $ neighbors g v
+    in 
+    sum [weight g x * weight g y | x <- nodes, y <- nodes, x /= y, not (isLinkedWithAnEdge g x y)]
+
 -- | Weight of a node
 weight :: (UndirectedGraph g, Factor f)
        => g a f 
@@ -101,6 +116,13 @@
 instance Show VertexCluster where 
     show (VertexCluster s) = show . Set.toList $ s
 
+
+
+instance IsCluster Cluster where 
+  overlappingEvidence (Cluster c) e = filter (\x -> Set.member (instantiationVariable x) c) e
+  clusterVariables (Cluster s) = Set.toList s
+  mkSeparator (Cluster sa) (Cluster sb) = Cluster $ Set.intersection sa sb
+
 -- | Triangulate a graph using a cost function
 -- The result is the triangulated graph and the list of clusters
 -- which may not be maximal.
@@ -153,155 +175,6 @@
       (Nothing,_) -> checkIfMaximal (current:reversedPrefix) (head suffix) (tail suffix)
       (Just r,l) -> checkIfMaximal (r:reverse l) (head suffix) (tail suffix)
 
-
--- | Create the cluster graph
-createClusterGraph :: UndirectedGraph g
-                   => [VertexCluster] 
-                   -> g Int VertexCluster
-createClusterGraph c =
-  let numberedClusters = zip c (map Vertex [0..])
-      addCluster (c,v) g = addVertex v c g
-      graphWithoutEdges = foldr addCluster emptyGraph numberedClusters
-      separatorSize ca cb = Set.size $ Set.intersection (fromVertexCluster ca) (fromVertexCluster cb)
-      allEdges = [(cx,cy) | cx <- numberedClusters, cy <- numberedClusters, cx /= cy]
-      addClusterEdge ((ca,va),(cb,vb)) g = addEdge (edge va vb) (separatorSize ca cb) g
-  in 
-  foldr addClusterEdge graphWithoutEdges allEdges
-
-
-{-
-
-Minimum spanning tree using Prim's algorithm
-  
--}
-
--- | Tree with values on edges
-data Tree b a = Node a [(b,Tree b a)] deriving(Eq)
-
-{-
-
-Implementation of show for the tree
- 
--}
-standardHaskellTree :: (Show f, Show b) => Tree b (JTNodeValue f) -> T.Tree String 
-standardHaskellTree n@(Node a []) = T.Node (show $ nodeCluster n) []
-standardHaskellTree n@(Node a l) = T.Node (show $ nodeCluster n) (map (standardHaskellTree  . snd) l)
-
-standardVertexTree :: Tree () VertexCluster -> T.Tree String 
-standardVertexTree n@(Node a []) = T.Node (show a) []
-standardVertexTree n@(Node a l) = T.Node (show a) (map (standardVertexTree  . snd) l)
-  
-showFactorsAndEdges :: (Show f, Show b) => Tree b (JTNodeValue f) -> (String -> String) 
-showFactorsAndEdges  n@(Node a []) = (++ show (nodeValueFactor a))
-showFactorsAndEdges  n@(Node a l) = foldl1 (.) (map factorAndEdge l) . (++ show (nodeValueFactor a)) 
-  where 
-    factorAndEdge (s,t) = showFactorsAndEdges t . (++ show s) 
-
-instance (Show f ,Show b)=> Show (Tree b (JTNodeValue f)) where 
-  show t = "JUNCTION TREE\n" ++ T.drawTree (standardHaskellTree t) ++ "\n" ++ showFactorsAndEdges t "" ++ "\n------\n"
-
-instance Show (Tree () VertexCluster) where 
-  show t = "JUNCTION TREE\n" ++ T.drawTree (standardVertexTree t) ++ "\n"
-
-instance Functor.Functor (Tree b) where 
-  fmap f (Node a l) = Node (f a) (map (mapEdge f) l)
-    where 
-      mapEdge f (e,c) = (e, fmap f c)
-
--- | Expand a tree (encoded as a list of edges)
--- by adding vertices and keeping track of the vertices which have
--- already been added.
--- The selection of where to connect the new vertices is based upon cost of the new edges
-expand :: UndirectedGraph g 
-       => g Int f 
-       -> [Edge] -- ^ List of edges
-       -> [Vertex] -- ^ Vertices in Tree
-       -> [Vertex] -- ^ Vertices to add
-       -> [Edge] -- ^ Updated sets and edge list
-expand g theEdges inTree remaining | null remaining = theEdges
-                                   | otherwise = 
-                                        let (treeVertex,outVertex) = maximumBy (compare `on` (edgeCost g)) $ [(vin,vout) | vin <- inTree, vout <-remaining,isLinkedWithAnEdge g vin vout]
-                                        in 
-                                        expand g (edge treeVertex outVertex : theEdges) (outVertex : inTree)
-                                          (filter (/= outVertex) remaining)
-
-  where 
-    edgeCost g (va,vb) = fromJust $ edgeValue g (edge va vb)
-
-leaf x = Node x []
-treeEdge c b = (c,b)
-
--- | Create a tree based upon a description with edges
-createTreeFromMap :: Vertex -- ^ Root vertex
-                  -> Map.Map Vertex [Vertex] -- ^ Tree branches
-                  -> Tree () Vertex 
-createTreeFromMap root m = 
-  let growTree m t@(Node a _) | Map.null m = t
-                              | otherwise = 
-                                    case Map.lookup a m of 
-                                      Nothing -> t 
-                                      Just l -> Node a . map (treeEdge () . growTree m . leaf) $ l
-  in
-  growTree m (leaf root)
-                   
--- | Implementing the Prim's algorithm for minimum spanning tree
-minimumSpanningTree :: UndirectedGraph g 
-                    => g Int f 
-                    -> Tree () f 
-minimumSpanningTree g = 
-  let startRoot = fromJust $ someVertex g 
-      remainingVertices = filter (/= startRoot) (allVertices g)
-      foundEdges = expand g [] [startRoot] remainingVertices
-      m = Map.fromListWith (++) . map ((\(a,b) -> (a,[b])) . edgeEndPoints) $ foundEdges
-      theTree = createTreeFromMap startRoot m
-  in 
-  Functor.fmap (fromJust . vertexValue g) theTree
-      
-   
-{-
-
-Junction tree algorithm
-
--}
-
--- | Check if all variables of a factor are included in a cluster
-vertexClusterIsContainingFactor :: Factor f => VertexCluster -> f -> Bool 
-vertexClusterIsContainingFactor c f = 
-  let factorVars = Set.fromList . map variableVertex . factorVariables $ f
-  in 
-  Set.isSubsetOf factorVars (fromVertexCluster c)
-
--- | Check if all variables of a factor are included in a cluster
-clusterIsContainingVariable :: DV -> Cluster  -> Bool 
-clusterIsContainingVariable v c  =  
-  Set.member v (Set.fromList $ fromCluster c)
-
--- | Separator which can be in 3 state depending how many messages have passed through it
-data Separator f = NoMessage !Cluster
-                 | Collect !Cluster !f 
-                 | Distribute !Cluster !f !f -- Upward and downward message
-                 deriving(Eq)
-
-instance Show f => Show (Separator f) where 
-  show (NoMessage c) = "NoMessage: " ++ show c 
-  show (Collect c u) = "Collect: " ++ show c ++ "\n" ++ "\n <----- \n" ++ show u ++ "\n"
-  show (Distribute c u d) = "Distribute: " ++ show c ++ "\n <----- \n" ++ show u ++ "\n" ++ " -----> \n" ++ show d ++ "\n"
-
-
--- | Evidence if some is used for the node
-type Evidence f = f
-
--- | Evidence for cluster, factor for cluster
-data JTNodeValue f = JTNodeValue !Cluster !(Evidence f) !f deriving(Eq,Show)
-
--- | Cluster of discrete variables.
--- Discrete variables instead of vertices are needed because the
--- factor are using 'DV' and we need to find
--- which factors must be contained in a given cluster.
-newtype Cluster = Cluster (Set.Set DV) deriving(Eq,Show)
-
-fromCluster (Cluster s) = Set.toList s 
-
 -- | Convert the clusters from vertex to 'DV' clusters
 vertexClusterToCluster :: (Factor f , Graph g)
                        => g e f 
@@ -313,150 +186,115 @@
   in 
   Cluster . Set.fromList $ variables
 
--- | Vertices contained in a cluster
-clusterVertices :: VertexCluster -> [Vertex]
-clusterVertices = Set.toList . fromVertexCluster
 
--- | Find all factors contained in a cluster
-findFactorsForCluster :: (Factor f , Graph g)
-                      => BayesianNetwork g f
-                      -> VertexCluster
-                      -> [f]
-findFactorsForCluster g c = 
-  filter (vertexClusterIsContainingFactor c) . mapMaybe (vertexValue g) . clusterVertices $ c
-
--- | The junction tree
-type JunctionTree f = Tree (Separator f) (JTNodeValue f)
-
--- | Get the potential for a cluster
-mkNodePotential :: (Graph g, Factor f, Show f)
-                => BayesianNetwork g f 
-                -> VertexCluster 
-                -> Set.Set Vertex
-                -> (JTNodeValue f, Set.Set Vertex)
-mkNodePotential g c set =  
-  let -- Factor found in a cluster but they may already be used in another cluster
-      foundFactors = findFactorsForCluster g c
-      -- Get the vertices for the factor
-      vertexForFactors = map (variableVertex . factorMainVariable) foundFactors 
-      -- Keep only the factors which are not already used
-      isNotUsed (v,f) = Set.member v set
-      factorsNotYetUsed = filter isNotUsed (zip vertexForFactors foundFactors)
-      set' = Set.difference set (Set.fromList $ map fst factorsNotYetUsed)
-      factorsToUse = map snd factorsNotYetUsed
-    
-      potential = factorProduct factorsToUse
-  in 
-  (JTNodeValue (vertexClusterToCluster g c) (factorFromScalar 1.0) potential, set')
-
--- | Generate the evidence potential for a given cluster
-evidenceForCluster :: Factor f => DVISet Int -> Cluster -> Maybe (Evidence f)
-evidenceForCluster assignments cluster@(Cluster c) = 
-  let c' = Set.fromList (map instantiationVariable assignments) 
-      common = Set.intersection c' c 
-      selectedVariables = filter (\c -> Set.member (instantiationVariable c) common) assignments
+-- | Create the cluster graph
+createClusterGraph :: (UndirectedGraph g, Factor f, Graph g')
+                   => g' e f
+                   -> [VertexCluster] 
+                   -> g Int Cluster
+createClusterGraph bn c =
+  let numberedClusters = zip c (map Vertex [0..])
+      addCluster g (c,v)  = addVertex v (vertexClusterToCluster bn c) g
+      graphWithoutEdges = foldl' addCluster emptyGraph numberedClusters
+      separatorSize ca cb = Set.size $ Set.intersection (fromVertexCluster ca) (fromVertexCluster cb)
+      allEdges = [(cx,cy) | cx <- numberedClusters, cy <- numberedClusters, cx /= cy]
+      addClusterEdge g ((ca,va),(cb,vb)) = addEdge (edge va vb) (separatorSize ca cb) g
   in 
-  evidenceFrom selectedVariables
+  foldl' addClusterEdge graphWithoutEdges allEdges
 
 
--- | Get the cluster for a node
-nodeCluster :: Tree a (JTNodeValue f) -> Cluster 
-nodeCluster (Node (JTNodeValue c _ _ ) _) = c 
-
-emptyCluster :: Cluster 
-emptyCluster = Cluster Set.empty
-
-nodeValueFactor (JTNodeValue _ _ f ) = f
-nodeValueEvidence (JTNodeValue _ e _) = e
-
-nodeValueWithNewEvidence (JTNodeValue a e b) e' = JTNodeValue a e' b
-clearNodeValueEvidence (JTNodeValue a _ b)  = JTNodeValue a (factorFromScalar 1.0) b
-
--- | Get the cluster for a separator
-separatorCluster :: Separator f -> Cluster 
-separatorCluster (NoMessage c) = c
-separatorCluster (Collect c _) = c 
-separatorCluster (Distribute c _ _) = c 
-
+{-
 
-upMessage (Distribute _ u _) = Just u 
-upMessage (Collect _ u ) = Just u 
-upMessage _ = Nothing 
+Maximum spanning tree using Prim's algorithm
+  
+-}
 
-downMessage (Distribute _ _ d) = Just d 
-downMessage _ = Nothing 
+-- | Get all possible edges between the leaves and the remaining nodes
+possibilities :: (Ord c , UndirectedGraph g) 
+              => g Int c -- ^ Original graph to get the edge value 
+              -> JTree c (Vertex,f) -- ^ Tree to get the vertex for a leaf
+              -> [Vertex] -- ^ Vertices to add to the tree
+              -> [c] -- ^ List of leaves
+              -> [(Vertex,c,Int)] -- ^ Found edge to add
+possibilities g currentT remaining leavesClusters = do 
+  rv <- remaining
+  lv <- leavesClusters
+  let NodeValue (lvVertex,lvCluster) _ = nodeValue currentT lv
+  guard (isLinkedWithAnEdge g rv lvVertex)
+  let ev = fromJust $ edgeValue g (edge rv lvVertex)
+  return $ (rv,lv,ev)
 
-computeSeparatorCluster :: (Factor f, Graph g) 
-                        => BayesianNetwork g f 
-                        -> VertexCluster 
-                        -> VertexCluster
-                        -> Cluster
-computeSeparatorCluster g parent child = 
-  let theNodeCluster (Node c _) = c 
-      childVertices = fromVertexCluster child 
-      parentVertices = fromVertexCluster parent
-      separatorVertices = VertexCluster $ Set.intersection childVertices parentVertices
+-- | Find the max edge to add to the tree
+findMax :: (UndirectedGraph g, Ord c, Factor f)
+        => g Int c -- ^ Graph
+        -> [Vertex] -- ^ Nodes to add 
+        -> JTree c (Vertex,f)
+        -> ([Vertex],(Vertex,c),c) 
+findMax g remaining currentT = 
+  let leavesClusters = treeNodes currentT
+      edgeValue (_,_,e) = e
+      (rf,lf,ef) = maximumBy (compare `on` edgeValue) (possibilities g currentT remaining leavesClusters)
+      remaining' = filter (/= rf) remaining 
+      foundCluster = fromJust $ vertexValue g rf
   in
-  vertexClusterToCluster g  separatorVertices
+  (remaining', (rf, foundCluster), lf)
 
-dfs :: (n -> n -> e -> e') -- Parent, child node and their egde
-    -> (n -> a -> (n', a)) -- Node and current value -> new value and new nod
-    -> Tree e n  -- Tree to traverse
-    -> a -- Start value
-    -> (Tree e' n', a) -- New tree and new value
-dfs edgef nodef n@(Node nodevalue []) current = 
-  let (newnodevalue, newval) = nodef nodevalue current 
-  in 
-  (Node newnodevalue [],newval) 
-dfs edgef nodef n@(Node nodevalue children) current =
-  let (newnodevalue, newval) = nodef nodevalue current 
-      applyEdgeFunction (e,Node childvalue _) = edgef nodevalue childvalue e
-      applyToChildren childrenNode val = dfs edgef nodef childrenNode val
-      edges' = map applyEdgeFunction children
-      recurseOnChildren s r [] = (s,reverse r)
-      recurseOnChildren s r (a:l) = 
-        let (a',s') = applyToChildren a s
-        in 
-        recurseOnChildren s' (a':r) l 
-      (lastval,newSubTrees) = recurseOnChildren newval [] (map snd children)
-  in 
-  (Node newnodevalue (zip edges' newSubTrees),lastval)
+removeVertices :: JTree c (Vertex,f) -> JTree c f
+removeVertices t = t { nodeValueMap = Map.map removeVertexFromNode (nodeValueMap t)
+                     , separatorValueMap = Map.map removeVertexFromSeparator (separatorValueMap t)
+                     }
+ where 
+  removeVertexFromNode (NodeValue (_,f) (_,e)) = NodeValue f e 
+  removeVertexFromSeparator (SeparatorValue (_,u) (Just (_,d))) = SeparatorValue u (Just d)
+  removeVertexFromSeparator (SeparatorValue (_,u) Nothing) = SeparatorValue u Nothing 
+  removeVertexFromSeparator EmptySeparator = EmptySeparator
 
-setFactorEdgeUpdate :: (Graph g, Factor f) 
-                    => BayesianNetwork g f 
-                    -> VertexCluster 
-                    -> VertexCluster
-                    -> () 
-                    -> Separator f
-setFactorEdgeUpdate g parentvalue childvalue _ = NoMessage $ computeSeparatorCluster g parentvalue childvalue 
+-- | Implementing the Prim's algorithm for minimum spanning tree
+maximumSpanningTree :: (UndirectedGraph g, IsCluster c, Factor f, Ord c) 
+                    => g Int c 
+                    -> JTree c f
+maximumSpanningTree g = 
+    let rootNodeVertex = fromJust $ someVertex g 
+        rootNodeValue = fromJust $ vertexValue g rootNodeVertex
+        unitFactor = factorFromScalar 1.0 
+        startTree = singletonTree rootNodeValue (rootNodeVertex,unitFactor) (rootNodeVertex,unitFactor) 
+        remainingVertices = filter (/= rootNodeVertex) (allVertices g) 
+    in 
+    removeVertices $ buildTree g remainingVertices startTree 
 
-setFactorNodeUpdate :: (Graph g, Factor f, Show f) 
-                    => BayesianNetwork g f 
-                    -> VertexCluster
-                    -> Set.Set Vertex 
-                    -> (JTNodeValue f, Set.Set Vertex)
-setFactorNodeUpdate g nodeValue set = mkNodePotential g nodeValue set
+buildTree :: (UndirectedGraph g , IsCluster c, Factor f, Ord c)
+          => g Int c 
+          -> [Vertex]
+          -> JTree c (Vertex,f) 
+          -> JTree c (Vertex,f) 
+buildTree g [] currentT = currentT 
+buildTree g l@(h:t) currentT = 
+    let unitFactor = factorFromScalar 1.0
+        (l',(foundElemVertex,foundElemValue),leaf) = findMax g l currentT
+        sep = mkSeparator foundElemValue leaf
+        newTree = addSeparator leaf sep foundElemValue . 
+                  addNode foundElemValue (foundElemVertex,unitFactor) (foundElemVertex,unitFactor) $ currentT
+    in 
+    buildTree g l' newTree
+   
+{-
 
--- | Set a factor for a node
-setFactors :: (Graph g, Factor f, Show f)
-           => BayesianNetwork g f -- ^ Bayesian graph
-           -> Tree () VertexCluster  -- ^ Cluster tree with no factors
-           -> Set.Set Vertex
-           -> (JunctionTree f,Set.Set Vertex) -- ^ Initialized junction tree
-setFactors g = dfs (setFactorEdgeUpdate g) (setFactorNodeUpdate g) 
+Junction tree algorithm
 
+-}
 
+
 -- | Create a junction tree with only the clusters and no factors
-createVerticesJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g)
-                           => (UndirectedSG () b -> Vertex -> Vertex -> Ordering) -- ^ Weight function on the moral graph
-                           -> g () b -- ^ Input directed graph
-                           -> Tree () VertexCluster -- ^ Junction tree
-createVerticesJunctionTree cmp g =  
+createUninitializedJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f)
+                                => (UndirectedSG () f -> Vertex -> Vertex -> Ordering) -- ^ Weight function on the moral graph
+                                -> g () f -- ^ Input directed graph
+                                -> JunctionTree f -- ^ Junction tree
+createUninitializedJunctionTree cmp g =  
   let theMoralGraph = moralGraph g
       (clusters,_) = triangulate (cmp theMoralGraph) theMoralGraph
-      g'' = createClusterGraph clusters :: UndirectedSG Int VertexCluster
+      g'' = createClusterGraph g clusters :: UndirectedSG Int Cluster
   in 
-  minimumSpanningTree g''
+  maximumSpanningTree g''
 
 -- | Create a function tree
 createJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f, Show f)
@@ -464,137 +302,63 @@
                   -> BayesianNetwork g f -- ^ Input directed graph
                   -> JunctionTree f -- ^ Junction tree
 createJunctionTree cmp g = 
-  let cTree = createVerticesJunctionTree cmp g 
-      factorSet = Set.fromList (allVertices g) -- Tracking of factors which have not yet been put in the junction tree
+  let cTree = createUninitializedJunctionTree cmp g 
       -- A vertex is linked with a factor so vertex is used as the identifier
-      (newTree, _) = setFactors g cTree factorSet
-  in 
-  distribute Nothing . collect $ newTree
-
-
-collectMessages :: Factor f => (Separator f , JunctionTree f) -> (Separator f , JunctionTree f)
-collectMessages (separator, Node nc []) = 
-  let sc = separatorCluster separator
-      newPotential = factorProduct [nodeValueFactor nc,nodeValueEvidence nc] 
-      newMessage = factorProjectTo (fromCluster sc) newPotential
-  in
-  (Collect sc newMessage, Node nc []) -- Copy node factor to node current potential
-collectMessages (separator,(Node nc l)) = 
-  let sc = separatorCluster separator
-      messagesFromSubTrees = map collectMessages l 
-      newPotential = factorProduct (nodeValueEvidence nc:nodeValueFactor nc:(mapMaybe (upMessage . fst) messagesFromSubTrees))
-      newMessage = factorProjectTo (fromCluster sc) newPotential 
-  in 
-  (Collect sc newMessage, Node nc messagesFromSubTrees)
-
--- | Collect phase of the junction tree
-collect :: Factor f => JunctionTree f -> JunctionTree f 
-collect t = let (_,t') = collectMessages (NoMessage emptyCluster, t) in t'
-
-notSameCluster a b = nodeCluster a /= nodeCluster b 
-
--- | Distribute phase of the junction tree
-distribute :: Factor f => Maybe (Separator f) -> JunctionTree f -> JunctionTree f 
-distribute down n@(Node nc []) = n
-distribute down (Node nc l) = 
-  let receivedDownMessage = if isJust down then fromJust . downMessage . fromJust $ down else factorFromScalar 1.0
-      getUpMessage (edge,c) = upMessage edge 
-      upMessagesForSendingTo i = fromJust . mapM getUpMessage . filter ((i `notSameCluster`) . snd) $ l
-      newPotential i = factorProduct (nodeValueFactor nc:nodeValueEvidence nc:receivedDownMessage:upMessagesForSendingTo i)
-      newMessage sc i = factorProjectTo (fromCluster sc) (newPotential i)
-      distributeMessage s@(Collect sc dm,i) = 
-        let newSeparator = Distribute sc dm (newMessage sc i)
-        in 
-        (newSeparator,distribute (Just newSeparator) i)
-      distributeMessage _ = error "Distribute message can only update a collect phase message"
-      subTrees = map distributeMessage l
+      newTree = setFactors g cTree
   in 
-  Node nc subTrees
-
--- | Depth first search in  tree
-findInTree :: (Tree edge a -> Bool) -> Maybe edge -> Tree edge a -> Maybe (Maybe edge,Tree edge a)
-findInTree cmp e n@(Node a []) = if (cmp n) then Just (e,n) else Nothing 
-findInTree cmp e n@(Node a l) = 
-  let findSome [] = Nothing
-      findSome ((e',h):t) = 
-        case findInTree cmp (Just e') h of 
-          Nothing -> findSome t 
-          Just r -> Just r
-  in
-  case cmp n of 
-    True -> Just (e,n) 
-    False -> findSome l
+  distribute . collect $ newTree
 
 
 -- | Compute the marginal posterior (if some evidence is set on the junction tree)
   -- otherwise compute just the marginal prior.
 posterior :: Factor f => JunctionTree f -> DV -> Maybe f
-posterior t v = do 
-  (maybeEdge,Node n l) <- findInTree (clusterIsContainingVariable v . nodeCluster) Nothing t
-  let receivedDownMessage = maybe (factorFromScalar 1.0) id $ 
-                               do
-                                 e <- maybeEdge 
-                                 downMessage e
-      upMessages = fromJust . mapM (upMessage . fst) $ l
-      p = factorProduct (receivedDownMessage:nodeValueEvidence n:nodeValueFactor n:upMessages)
-  return $ normedFactor $ factorProjectTo [v] p 
+posterior t v = 
+  case snd $ traverseTree (findClusterFor v) Nothing t of 
+    Nothing -> Nothing
+    Just c -> let NodeValue f e = nodeValue t c 
+                  d = maybe (factorFromScalar 1.0) id $ downMessage t =<< (nodeParent t c)
+                  u = map (upMessage t) (nodeChildren t c)
+                  unNormalized = factorProjectTo [v] (factorProduct (f:e:d:u))
+              in 
+              Just $ factorDivide unNormalized (factorNorm unNormalized)
 
--- | Apply some evidence modifications in the tree
-applyEvidenceWith :: (JunctionTree f -> JunctionTree f) -- ^ Node modification function. Only change node value. Not the children
-                  -> JunctionTree f -- ^ Input tree
-                  -> JunctionTree f
-applyEvidenceWith nodeChange n@(Node _ []) = nodeChange n 
-applyEvidenceWith nodeChange n@(Node _ l) =
-  let Node n' l' = nodeChange n 
-      changeChildren (e,c) = (e,applyEvidenceWith nodeChange c)
-  in 
-  Node n' (map changeChildren l')
+-- | Find a cluster containing the variable
+findClusterFor :: DV 
+               -> Maybe Cluster
+               -> Cluster -- ^ Current cluster
+               -> NodeValue f -- ^ Current value
+               -> Action (Maybe Cluster) (NodeValue f)
+findClusterFor dv s c@(Cluster sc) v = 
+  case Set.member dv sc of 
+    False -> Skip s 
+    True -> Stop (Just c)
 
--- | Change the evidence for a node
-evidenceWith :: Factor f => DVISet Int -> JunctionTree f -> JunctionTree f
-evidenceWith assignments t@(Node n l) = 
-  let n' = case evidenceForCluster assignments (nodeCluster t) of 
-             Nothing -> n 
-             Just e' -> nodeValueWithNewEvidence n e'
-  in 
-  Node n' l
 
--- | Remove the evidence for a node
-clearNodeEvidence (Node n l) = Node (clearNodeValueEvidence n) l 
-
--- | Remove evidence in the junction tree
-clearEvidence :: Factor f => JunctionTree f -> JunctionTree f
-clearEvidence = distribute Nothing . collect . applyEvidenceWith (clearNodeEvidence)
-
--- | Update evidence in the tree
-updateEvidence :: Factor f => DVISet Int -> JunctionTree f -> JunctionTree f
-updateEvidence assignments = distribute Nothing . collect . applyEvidenceWith (evidenceWith assignments)
-
--- | Used to implement quickcheck.
--- The junction tree property is the property that CA intersection CB is included in all clusters in the path
--- from CA to CB.
-junctionTreeProperty :: [VertexCluster] -> Tree () VertexCluster -> Bool
-junctionTreeProperty path (Node _ []) = True 
-junctionTreeProperty path (Node c l) = 
-  let children = map snd l 
-  in
-  checkPath c (reverse path) && all (junctionTreeProperty (c:path)) children 
-
-junctionTreeProperty_prop :: DirectedSG () String -> Property 
+junctionTreeProperty_prop :: DirectedSG () CPT -> Property 
 junctionTreeProperty_prop g = (not . isEmpty) g && (not . hasNoEdges) g && connectedGraph g ==> 
   let cmp ug = (compare `on` (numberOfAddedEdges ug))
+      t = createUninitializedJunctionTree cmp g
   in
-  junctionTreeProperty [] (createVerticesJunctionTree cmp g)
+  junctionTreeProperty t [] (root t)
 
--- | Check that the intersection of C with any parent in included in any cluster between the parent and C.
-checkPath :: VertexCluster -> [VertexCluster] -> Bool 
-checkPath c l = 
-  let parentSets = map fromVertexCluster l
-      allIntersections = map (Set.intersection (fromVertexCluster c)) parentSets
-      pathsToEachParent = tail . inits $ parentSets
-      isSubsetOfAllParents i parents = all (Set.isSubsetOf i) parents
+junctionTreeProperty :: JTree Cluster CPT -> [Cluster] -> Cluster -> Bool
+junctionTreeProperty t path c = 
+  let cl = map (separatorChild t) . nodeChildren t $ c
+  in
+  checkPath c path && all (junctionTreeProperty t (c:path)) cl 
+
+
+-- | Check that the intersection of C with any parent in included in all cluster between the parent and C.
+checkPath :: Cluster -> [Cluster] -> Bool 
+checkPath _ [] = True
+checkPath (Cluster c) l = 
+  let clusterSet (Cluster s) = s -- x
+      parentSets = map clusterSet l -- Example a b c d where a is the root
+      allIntersectionsWithParents = map (Set.intersection c) parentSets -- a ^ x, b ^ x , c ^ x , d ^ x
+      pathsToEachParent = tail . inits $ parentSets -- a, ab, abc, abcd
+      isSubsetOfAllParents i path = all (Set.isSubsetOf i) path
   in    
-  and $ zipWith isSubsetOfAllParents allIntersections pathsToEachParent
+  and $ zipWith isSubsetOfAllParents allIntersectionsWithParents pathsToEachParent
 {-
 
 Moral graph
@@ -625,8 +389,8 @@
     (originGraph',dstGraph')
 
 -- | Add the missing parent links
-addMissingLinks :: DirectedGraph g => Vertex -> b -> g () b -> g () b
-addMissingLinks v _ g = 
+addMissingLinks :: DirectedGraph g => g () b -> Vertex -> b -> g () b
+addMissingLinks g v _ = 
     let (_,g') = connectAllNodesWith g g (\e m -> m) (\e m -> addEdge e () m) (parents g v)
     in 
     g'
@@ -637,11 +401,11 @@
                     => g  () b 
                     -> g' () b 
 convertToUndirected m = 
-    let addVertexWithLabel v dat g = 
+    let addVertexWithLabel g v dat  = 
            let theName = fromJust $ vertexLabel m v
            in 
            addLabeledVertex theName v dat g
-        newDiscreteGraph = foldrWithVertex addVertexWithLabel emptyGraph m
+        newDiscreteGraph = foldlWithVertex' addVertexWithLabel emptyGraph m
         addEmptyEdge edge g = addEdge edge () g
     in 
     foldr addEmptyEdge newDiscreteGraph . allEdges $ m
@@ -651,4 +415,4 @@
 moralGraph :: (NamedGraph g, FoldableWithVertex g, DirectedGraph g) 
            => g () b -> UndirectedSG () b 
 moralGraph g = 
-    convertToUndirected  . foldrWithVertex addMissingLinks g $ g
+    convertToUndirected  . foldlWithVertex' addMissingLinks g $ g
diff --git a/Bayes/FactorElimination/JTree.hs b/Bayes/FactorElimination/JTree.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/FactorElimination/JTree.hs
@@ -0,0 +1,486 @@
+{- | Junction Trees 
+
+The Tree data structures are not working very well with message passing algorithms. So, junction trees are using
+a different representation
+
+-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+module Bayes.FactorElimination.JTree(
+      IsCluster(..)
+    , Cluster(..)
+    , JTree(..)
+    , JunctionTree(..)
+    , setFactors
+    , distribute 
+    , collect
+    , fromCluster
+    , changeEvidence
+    , nodeIsMemberOfTree
+    , singletonTree
+    , addNode 
+    , addSeparator
+    , leaves
+    , nodeValue
+    , NodeValue(..)
+    , SeparatorValue(..)
+    , downMessage
+    , upMessage 
+    , nodeParent 
+    , nodeChildren
+    , traverseTree
+    , separatorChild
+    , treeNodes
+    , Action(..)
+    ) where 
+
+import qualified Data.Map as Map
+import qualified Data.Tree as Tree
+import Data.Maybe(fromJust,mapMaybe)
+import qualified Data.Set as Set
+import Data.Monoid
+import Data.List((\\), intersect,partition, foldl')
+import Bayes.PrivateTypes 
+import Bayes.Factor
+import Bayes
+
+import Debug.Trace 
+debug s a = trace (s ++ " " ++ show a ++ "\n") a
+
+type UpMessage a = a 
+type DownMessage a = Maybe a
+
+-- | Separator value
+data SeparatorValue a = SeparatorValue !(UpMessage a) !(DownMessage a)
+                      | EmptySeparator -- ^ Use to track the progress in the collect phase
+                      deriving(Eq)
+
+instance Show a => Show (SeparatorValue a) where 
+    show EmptySeparator = ""
+    show (SeparatorValue u Nothing) = "u(" ++ show u ++ ")"
+    show (SeparatorValue u (Just d)) = "u(" ++ show u ++ ") d(" ++ show d ++ ")"
+
+type FactorValue a = a 
+type EvidenceValue a = a
+
+-- | Node value
+data NodeValue a = NodeValue !(FactorValue a) !(EvidenceValue a) deriving(Eq)
+
+instance Show a => Show (NodeValue a) where 
+    show (NodeValue f e) = "f(" ++ show f ++ ") e(" ++ show e ++ ")"
+
+-- | Junction tree.
+-- 'c' is the node / separator identifier (for instance a set of 'DV')
+-- a are the values for a node or separator
+data JTree  c f = JTree {  root :: !c
+                        -- | Leaves of the tree
+                        ,  leavesSet :: !(Set.Set c)
+                        -- | The children of a node are separators
+                        ,  childrenMap :: !(Map.Map c [c])
+                        -- | Parent of a node
+                        ,  parentMap :: !(Map.Map c c)
+                        -- | Parent of a separator
+                        ,  separatorParentMap :: !(Map.Map c c)
+                        -- | The child of a seperator is a node
+                        ,  separatorChildMap :: !(Map.Map c c)
+                        -- | Values for nodes and seperators
+                        ,  nodeValueMap :: !(Map.Map c (NodeValue f))
+                        ,  separatorValueMap :: !(Map.Map c (SeparatorValue f))
+                        } deriving(Eq)
+
+-- | Create a singleton tree with just one root node
+singletonTree r factorValue evidenceValue = 
+    let t = JTree r Set.empty Map.empty Map.empty Map.empty Map.empty Map.empty Map.empty
+    in 
+    addNode r factorValue evidenceValue t
+
+-- | Leaves of the tree
+leaves :: JTree c a -> [c]
+leaves = Set.toList . leavesSet
+
+-- | All nodes of the tree
+treeNodes :: JTree c a -> [c]
+treeNodes = Map.keys . nodeValueMap
+
+-- | Value of a node
+nodeValue :: Ord c => JTree c a -> c -> NodeValue a 
+nodeValue t e = fromJust $ Map.lookup e (nodeValueMap t)
+
+-- | Change the value of a node
+setNodeValue :: Ord c => c -> NodeValue a -> JTree c a -> JTree c a
+setNodeValue c v t = t {nodeValueMap = Map.insert c v (nodeValueMap t)} 
+
+-- | Parent of a node
+nodeParent :: Ord c => JTree c a -> c -> Maybe c 
+nodeParent t e = Map.lookup e (parentMap t)
+
+-- | Value of a node
+separatorValue :: Ord c => JTree c a -> c -> SeparatorValue a 
+separatorValue t e = fromJust $ Map.lookup e (separatorValueMap t)
+
+-- | Parent of a separator
+separatorParent :: Ord c => JTree c a -> c -> c 
+separatorParent t e = fromJust $ Map.lookup e (separatorParentMap t)
+
+-- | UpMessage for a separator node
+upMessage :: Ord c => JTree c a -> c -> a
+upMessage t c = case separatorValue t c of 
+                  SeparatorValue up _ -> up 
+                  _ -> error "Trying to get an up message on an empty seperator ! Should never occur !"
+
+-- | DownMessage for a separator node
+downMessage :: Ord c => JTree c a -> c -> Maybe a 
+downMessage t c = case separatorValue t c of 
+     SeparatorValue _ (Just down) -> Just down 
+     SeparatorValue _ Nothing -> Nothing
+     _ -> error "Trying to get a down message on an empty separator ! Should never occur !"
+
+-- | Return the separator childrens of a node
+nodeChildren :: Ord c => JTree c a -> c -> [c]
+nodeChildren t e = maybe [] id $ Map.lookup e (childrenMap t)
+
+-- | Return the child of a separator
+separatorChild :: Ord c => JTree c a -> c -> c 
+separatorChild t e = fromJust $ Map.lookup e (separatorChildMap t)
+
+-- | Check if a node is member of the tree
+nodeIsMemberOfTree :: Ord c => c -> JTree c a -> Bool 
+nodeIsMemberOfTree c t = Map.member c (nodeValueMap t)
+
+-- | Add a separator between two nodes.
+-- The nodes MUST already be in the tree
+addSeparator :: (Ord c) 
+             => c -- ^ Origin node 
+             -> c -- ^ Separator
+             -> c -- ^ Destination node 
+             -> JTree c a -- ^ Current tree 
+             -> JTree c a -- ^ Modified tree 
+addSeparator node sep dest t = 
+    t { childrenMap = Map.insertWith' (++) node [sep] (childrenMap t)
+      , separatorChildMap = Map.insert sep dest (separatorChildMap t)
+      , separatorValueMap = Map.insert sep EmptySeparator (separatorValueMap t)
+      , leavesSet = Set.delete node (leavesSet t) 
+      , parentMap = Map.insert dest sep (parentMap t)
+      , separatorParentMap = Map.insert sep node (separatorParentMap t)
+      }
+
+-- | Add a new node
+addNode :: (Ord c) 
+        => c -- ^ Node
+        -> a -- ^ Factor value 
+        -> a -- ^ Evidence value
+        -> JTree c a 
+        -> JTree c a 
+addNode node factorValue evidenceValue t = 
+   t { nodeValueMap = Map.insert node (NodeValue factorValue evidenceValue) (nodeValueMap t)
+     , leavesSet = Set.insert node (leavesSet t) 
+     }
+
+-- | Update the up message of a separator
+updateUpMessage :: Ord c 
+                => Maybe c -- ^ Separator node to update (if any : none for root node)
+                -> a -- ^ New value
+                -> JTree c a -- ^ Old tree
+                -> JTree c a
+updateUpMessage Nothing _ t = t
+updateUpMessage (Just sep) newval t = 
+    let newSepValue =  case separatorValue t sep of 
+                         EmptySeparator -> SeparatorValue newval Nothing
+                         SeparatorValue up down -> SeparatorValue newval down 
+    in 
+    t {separatorValueMap = Map.insert sep newSepValue (separatorValueMap t)}
+
+-- | Update the down message of a separator
+updateDownMessage :: Ord c 
+                  => c -- ^ Separator node to update
+                  -> a -- ^ New value
+                  -> JTree c a -- ^ Old tree
+                  -> JTree c a
+updateDownMessage sep newval t = 
+    let newSepValue = case separatorValue t sep of 
+                        EmptySeparator -> error "Can't set a down message on an empty separator"
+                        SeparatorValue up _ -> SeparatorValue up (Just newval)
+    in 
+    t {separatorValueMap = Map.insert sep newSepValue (separatorValueMap t)}
+
+{-
+
+Message passing algorithms
+    
+-}
+
+-- | Functions used to generate new messages
+class Message f c | f -> c where
+    -- | Generate a new message from the received ones
+    newMessage :: [f] -> NodeValue f -> c -> f 
+
+
+-- | Check that a separator is initialized
+separatorInitialized :: SeparatorValue a -> Bool
+separatorInitialized EmptySeparator = False 
+separatorInitialized _ = True
+
+allSeparatorsHaveReceivedAMessage :: Ord c
+                                  => JTree c a -- ^ Tree
+                                  -> [c] -- ^ Separators
+                                  -> Bool 
+allSeparatorsHaveReceivedAMessage t seps = 
+  all separatorInitialized . map (separatorValue t) $ seps
+
+-- | Update the up separator by sending a message
+-- But only if all the down separators have received a message
+updateUpSeparator :: (Message a c, Ord c) 
+                  => JTree c a 
+                  -> c -- ^ Node generating the new upMessage
+                  -> JTree c a 
+updateUpSeparator t h  = 
+    let seps = nodeChildren t h
+    in
+    case allSeparatorsHaveReceivedAMessage t seps of 
+      False -> t 
+      True -> let incomingMessages = map (upMessage t) seps
+                  currentValue = nodeValue t h
+                  destinationNode = nodeParent t h
+              in 
+              case destinationNode of 
+                Nothing -> t -- When root
+                Just p -> let generatedMessage = newMessage incomingMessages currentValue p
+                          in 
+                          updateUpMessage destinationNode generatedMessage t
+
+-- | Update the down separator by sending a message
+updateDownSeparator :: (Message a c, Ord c) 
+                    => c -- ^ Node generating the message 
+                    -> JTree c a 
+                    -> c -- ^ Child receiving the message
+                    -> JTree c a 
+updateDownSeparator node t child  = 
+    let incomingMessagesFromBelow = map (upMessage t) (nodeChildren t node \\ [child])
+        messageFromAbove = downMessage t =<< (nodeParent t node)
+        incomingMessages = maybe incomingMessagesFromBelow (\x -> x:incomingMessagesFromBelow) messageFromAbove
+        currentValue = nodeValue t node
+        generatedMessage = newMessage incomingMessages currentValue child
+    in 
+    updateDownMessage child generatedMessage t
+
+unique :: Ord c => [c] -> [c]
+unique = Set.toList . Set.fromList
+
+data TraversalState = ACluster | ASeparator
+
+-- | Collect message taking into account that the tree depth may be different for different leaves.
+collect :: (Ord c, Message a c) 
+        => JTree c a 
+        -> JTree c a
+collect t = _collect ACluster (leaves t) t
+
+_collect :: (Ord c, Message a c) 
+         => TraversalState -- ^ Node processing phase or separator processing phase
+         -> [c]
+         -> JTree c a -- ^ Tree
+         -> JTree c a -- ^ Modified tree 
+_collect _ [] t = t
+_collect ACluster l t = 
+    let newTree = foldl' updateUpSeparator t l
+    in 
+    _collect ASeparator (mapMaybe (nodeParent t) l) newTree
+_collect ASeparator l t = _collect ACluster (unique . map (separatorParent t) $ l) t
+    
+distribute :: (Ord c, Message a c)
+           => JTree c a 
+           -> JTree c a
+distribute t = _distribute ACluster t (root t) 
+
+_distribute :: (Ord c, Message a c)
+            => TraversalState -- ^ True if node
+            -> JTree c a 
+            -> c -- ^ Destination of the distribute
+            -> JTree c a 
+_distribute ACluster t node = 
+    let children = nodeChildren t node
+        newTree = foldl' (updateDownSeparator node) t $ children
+    in
+    foldl' (_distribute ASeparator) newTree children
+_distribute ASeparator t node = _distribute ACluster t (separatorChild t node)
+
+{-
+
+Factors and evidence modifications
+
+-}
+
+-- | This class is used to check if evidence or a factor is relevant
+-- for a cluster
+class IsCluster c where 
+  overlappingEvidence :: c -> [DVI Int] -> [DVI Int]
+  clusterVariables :: c -> [DV]
+  mkSeparator :: c -> c -> c
+
+instance IsCluster [DV] where 
+  overlappingEvidence c e = filter (\x -> instantiationVariable x `elem` c) e
+  clusterVariables = id
+  mkSeparator = intersect
+
+data Action s a = Skip !s 
+                | ModifyAndStop !s !a
+                | Modify !s !a
+                | Stop !s
+
+-- | Traverse a tree and modify it
+traverseTree :: Ord c 
+             => (s -> c -> NodeValue f -> Action s (NodeValue f)) -- ^ Modification function
+             -> s -- ^ Current state
+             -> JTree c f -- ^ Input tree
+             -> (JTree c f,s)
+traverseTree action state t = _traverseTree True action (t,state) (root t)
+
+_traverseTree False action (t,state) current  = _traverseTree True action (t,state)  (separatorChild t current) 
+_traverseTree True action (t,state) current = 
+  case action state current (nodeValue t current) of 
+     Stop newState -> (t,newState)
+     ModifyAndStop _ newValue -> (setNodeValue current newValue t, state) 
+     Skip newState -> foldl' (_traverseTree False action) (t,newState) (nodeChildren t current)
+     Modify newState newValue -> 
+         let newTree = setNodeValue current newValue t 
+         in 
+         foldl' (_traverseTree False action) (newTree,newState) (nodeChildren newTree current)
+
+mapWithCluster :: Ord c 
+               => (c -> NodeValue f -> NodeValue f)
+               -> JTree c f 
+               -> JTree c f        
+mapWithCluster f t = t {nodeValueMap = Map.mapWithKey f (nodeValueMap t)}
+
+-- | Set the factors in the tree 
+setFactors :: (Graph g, Factor f, Show f, IsCluster c, Ord c)
+           => BayesianNetwork g f 
+           -> JTree c f 
+           -> JTree c f
+setFactors g t = 
+  let factors = allVertexValues g 
+  in
+  fst . traverseTree updateFactor factors $ t 
+
+
+-- | Update factors in a cluster
+updateFactor :: (Factor f, IsCluster c) 
+             => [f] -- ^ Remaining list of factors to attribute
+             -> c -- ^ Current cluster
+             -> NodeValue f -- ^ Current value
+             -> Action [f] (NodeValue f)
+updateFactor lf c (NodeValue _ evidence) | null lf = Stop lf
+                                         | otherwise =
+  let isFactorIncluded l = all (`elem` clusterVariables c) (factorVariables l)
+      (attributedFactors,remainingFactors) = partition isFactorIncluded lf 
+  in 
+  Modify remainingFactors  (NodeValue (factorProduct attributedFactors) evidence)
+
+-- | Change evidence in the network
+changeEvidence :: (IsCluster c, Ord c, Factor f, Message f c)
+               => [DVI Int] -- ^ Evidence
+               -> JTree c f 
+               -> JTree c f 
+changeEvidence e t =  distribute . 
+                      collect . fst .
+                      traverseTree changeNodeEvidence e $ 
+                      t { separatorValueMap = Map.map (const EmptySeparator) (separatorValueMap t)}
+
+changeNodeEvidence :: (IsCluster c, Factor f) 
+                   => [DVI Int] -- ^ Evidence
+                   -> c  -- ^ Current cluster
+                   -> NodeValue f -- ^ Current value
+                   -> Action [DVI Int] (NodeValue f)
+changeNodeEvidence [] c v = Stop []
+changeNodeEvidence e c (NodeValue f olde) = 
+  let oe = overlappingEvidence c e
+      ns = e \\ oe
+      newEvidence = factorProduct $ map factorFromInstantiation oe
+  in 
+  Modify ns (NodeValue f newEvidence)
+
+
+-- | Cluster of discrete variables.
+-- Discrete variables instead of vertices are needed because the
+-- factor are using 'DV' and we need to find
+-- which factors must be contained in a given cluster.
+newtype Cluster = Cluster (Set.Set DV) deriving(Eq,Ord)
+
+instance Show Cluster where 
+  show (Cluster s) = show . Set.toList $ s
+
+fromCluster (Cluster s) = Set.toList s 
+
+
+instance Factor f => Message f Cluster where 
+  newMessage input (NodeValue f e) dv = factorProjectTo (fromCluster dv) (factorProduct (f:e:input))
+
+
+type JunctionTree f = JTree Cluster f
+
+{-
+
+Implement the show function to see the structure of the tree
+(without the values)
+    
+-}
+
+data NodeKind c = N !c | S !c
+
+label True c a = c ++ "=" ++ show a 
+label False c _ = c
+
+-- | Convert the JTree into a tree of string
+-- using the cluster.
+toTree :: (Ord c, Show c, Show a) 
+       => Bool -- ^ True if the data must be displayed
+       -> JTree c a 
+       -> Tree.Tree String
+toTree d t = 
+    let r = root t
+        v = nodeValue t r
+        nodec = map S (nodeChildren t r)
+    in 
+    Tree.Node (label d (show r) v) (_toTree d t nodec)
+
+
+_toTree :: (Ord c, Show c, Show a) 
+        => Bool
+        -> JTree c a 
+        -> [NodeKind c] 
+        -> [Tree.Tree String]
+_toTree _ _ [] = []
+_toTree d t ((N h):l) = 
+    let nodec = map S (nodeChildren t h) -- Node children are separators
+        v = nodeValue t h
+    in
+    Tree.Node (label d (show h) v) (_toTree d t nodec):_toTree d t l
+_toTree d t ((S h):l) = 
+    let separatorc = [N $ separatorChild t h] -- separator child is a node
+        v = separatorValue t h
+    in
+    Tree.Node (label d ("<" ++ show h ++ ">") v ) (_toTree d t separatorc):_toTree d t l
+
+instance (Ord c, Show c, Show a) => Show (JTree c a) where 
+    show = Tree.drawTree . toTree False
+
+displayTree b = Tree.drawTree . toTree b
+
+{-
+
+Debug functions for tests
+
+-}
+
+--instance Message (Sum Int) String where 
+--    newMessage l (NodeValue a b) _ = mconcat (a:b:l)
+--
+--testTree :: JTree String (Sum Int)
+--testTree = let s a= Sum a
+--           in
+--           addSeparator "ROOT" "RB" "B" .
+--           addNode "B" (s 3) (s 3) . 
+--           addSeparator "ROOT" "RA" "A"  . 
+--           addNode "A" (s 2) (s 2) $ 
+--           singletonTree "ROOT" (s 4) (s 5)
+--
diff --git a/Bayes/ImportExport/HuginNet.hs b/Bayes/ImportExport/HuginNet.hs
--- a/Bayes/ImportExport/HuginNet.hs
+++ b/Bayes/ImportExport/HuginNet.hs
@@ -12,6 +12,7 @@
 import qualified Data.Map as Map
 import Bayes.Factor
 import Bayes
+import Bayes.PrivateTypes
 
 --import Debug.Trace 
 
@@ -139,7 +140,7 @@
         conds = tail dvs
         oldOrder = conds ++ [dst]
         dvalues = map read values :: [Double]
-        newvalues = changeVariableOrder oldOrder dvs dvalues
+        newvalues = changeVariableOrder (DVSet oldOrder) (DVSet dvs) dvalues
     cpt dst conds ~~ newvalues
     return ()
 
diff --git a/Bayes/PrivateTypes.hs b/Bayes/PrivateTypes.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/PrivateTypes.hs
@@ -0,0 +1,184 @@
+{- | Private types for Bayes and Factors.
+
+Those type are not exported
+
+-}
+
+module Bayes.PrivateTypes( 
+ -- * Classes
+   BayesianDiscreteVariable(..)
+ , Set(..)
+ -- * Variables
+ , DV(..)
+ , DVSet(..)
+ -- * Instantiations
+ , DVI(..)
+ , setDVValue
+ , (=:)
+ , instantiationValue
+ , instantiationVariable
+ -- * Vertices 
+ , Vertex(..)
+ -- * Misc
+ , getMinBound
+ -- * Indices 
+ , Index(..)
+ , forAllInstantiations 
+ , indicesForDomain
+ , fromIndex
+ , instantiationDetails
+ , allInstantiationsForOneVariable
+ ) where
+
+
+import qualified Data.List as L 
+
+{-
+	Set
+-}
+
+-- | A Set of variables used in a factor. s is the set and a the variable
+class Set s where
+    -- | Empty set
+    emptySet :: s a
+    -- | Union of two sets
+    union :: Eq a => s a -> s a -> s a
+    -- | Intersection of two sets
+    intersection :: Eq a => s a -> s a -> s a
+    -- | Difference of two sets
+    difference :: Eq a => s a -> s a -> s a
+    -- | Check if the set is empty
+    isEmpty :: s a -> Bool
+    -- | Check if an element is member of the set
+    isElem :: Eq a => a -> s a -> Bool
+    -- | Add an element to the set
+    addElem :: Eq a => a -> s a -> s a
+    -- | Number of elements in the set
+    nbElements :: s a -> Int
+
+    -- | Check if a set is subset of another one
+    subset :: Eq a => s a -> s a -> Bool
+
+    -- | Check set equality
+    equal :: Eq a => s a -> s a -> Bool
+    equal sa sb = (sa `subset` sb) && (sb `subset` sa)
+
+instance Set [] where
+    emptySet = []
+    union = L.union
+    intersection = L.intersect
+    difference a b = a L.\\ b
+    isEmpty [] = True 
+    isEmpty _ = False
+    isElem = L.elem 
+    addElem a l = if a `elem` l then l else a:l
+    nbElements = length
+    subset sa sb = all (`elem` sb) sa
+
+{-
+
+Misc
+
+-}
+-- | Vertex type used to identify a vertex in a graph
+newtype Vertex = Vertex {vertexId :: Int} deriving(Eq,Ord)
+
+instance Show Vertex where 
+    show (Vertex v) = "v" ++ show v
+
+-- | A discrete variable has a number of levels which is required to size the factors
+class BayesianDiscreteVariable v where
+    dimension :: v -> Int 
+
+-- | Get the minimum bound for a type
+getMinBound :: Bounded a => a -> a 
+getMinBound _ = minBound
+
+
+
+{-
+
+Variables
+	
+-}
+
+-- | A discrete variable
+data DV = DV !Vertex !Int deriving(Eq,Ord)
+
+-- | A set of discrete variables
+-- The tag is used to check that an index is used with the right set of DV
+newtype DVSet s = DVSet [DV] deriving(Eq)
+
+-- | Remove the type tag when not needed
+fromDVSet :: DVSet s -> [DV]
+fromDVSet (DVSet l) = l
+
+instance Show DV where
+    show (DV v d) = show v ++ "(" ++ show d ++ ")"
+
+instance BayesianDiscreteVariable DV where
+    dimension (DV _ d) = d
+
+{-
+
+Index
+
+-}
+
+-- | An index with meaning only for a given DVSet
+newtype Index s = Index Int deriving(Eq)
+
+-- | Used to forget the type tag
+fromIndex :: Index s -> Int 
+fromIndex (Index i) = i 
+
+-- | Generate all the indices for a set of variables
+indicesForDomain :: DVSet s -> [[Index s]]
+indicesForDomain (DVSet l) = mapM indicesForOneDomain l
+ where 
+ 	indicesForOneDomain (DV _ d) = map Index [0..d-1]
+
+allInstantiationsForOneVariable :: DV -> [DVI Int]
+allInstantiationsForOneVariable v@(DV _ d) = map (setDVValue v) [0..d-1]
+
+-- | Generate all instantiations of variables
+-- The DVInt can be in any order so the tag s is not used
+forAllInstantiations :: DVSet s -> [[DVI Int]]
+forAllInstantiations (DVSet l) = mapM allInstantiationsForOneVariable l
+ 
+
+{- 
+
+Instantiations
+
+-}
+-- | Discrete Variable instantiation. A variable and its value
+data DVI a = DVI DV !a deriving(Eq)
+
+instance Show a => Show (DVI a) where 
+   show (DVI (DV v _) i) = show v ++ "=" ++ show i
+
+   -- | A set of variable instantiations
+type DVISet a = [DVI a]
+
+-- | Create a discrete variable instantiation for a given discrete variable
+setDVValue :: DV -> a -> DVI a
+setDVValue v a = DVI v a
+
+-- | Create a variable instantiation using values from
+-- an enumeration
+(=:) :: (Bounded b, Enum b) => DV -> b -> DVI Int 
+(=:) a b = setDVValue a (fromEnum b - fromEnum (getMinBound b))
+
+instance BayesianDiscreteVariable (DVI a) where
+    dimension (DVI v _) = dimension v
+
+-- | Get the variables and their values with a type constraint
+instantiationDetails :: [DVI Int] -> (DVSet s, [Index s])
+instantiationDetails l = (DVSet $ map instantiationVariable l, map (Index . instantiationValue) l)
+
+-- | Extract value of the instantiation
+instantiationValue (DVI _ v) = v
+
+-- | Discrete variable from the instantiation
+instantiationVariable (DVI dv _) = dv
diff --git a/Bayes/Test.hs b/Bayes/Test.hs
--- a/Bayes/Test.hs
+++ b/Bayes/Test.hs
@@ -13,6 +13,10 @@
 import Bayes.Factor(testProductProject_prop,testScale_prop,testProjectCommut_prop,testScalarProduct_prop,testProjectionToScalar_prop)
 import Bayes.FactorElimination(junctionTreeProperty_prop)
 
+#ifdef LOCAL
+import Bayes.Test.ReferencePatterns(compareAsiaReference,compareCancerReference,comparePokerReference,compareFarmReference)
+#endif 
+
 -- | Run all the tests
 runTests = defaultMain tests
 
@@ -32,6 +36,14 @@
                 testProperty "Test the junction tree property" junctionTreeProperty_prop,
                 testCase "Test variable elimination == factor elimination" compareVariableFactor
             ]
+#ifdef LOCAL
+        , testGroup "Reference patterns" [ 
+                testCase "Asia reference pattern" compareAsiaReference,
+                testCase "Cancer reference pattern" compareCancerReference,
+                testCase "Poker reference pattern" comparePokerReference,
+                testCase "Farm reference pattern" compareFarmReference
+        ]
+#endif
 
     ]
 
diff --git a/Bayes/Test/CompareEliminations.hs b/Bayes/Test/CompareEliminations.hs
--- a/Bayes/Test/CompareEliminations.hs
+++ b/Bayes/Test/CompareEliminations.hs
@@ -32,8 +32,8 @@
         jt = createJunctionTree nodeComparisonForTriangulation exampleG
     compareFactors "PRIOR FOR RAIN" (posterior jt rain) (priorMarginal exampleG [winter,sprinkler,wet,road] [rain])
 
-    let jt1 = updateEvidence [wet =: True] jt 
-        jt2 = updateEvidence [wet =: True, sprinkler =: True] jt1 
+    let jt1 = changeEvidence [wet =: True] jt 
+        jt2 = changeEvidence [wet =: True, sprinkler =: True] jt1 
 
     compareFactors "POSTERIOR RAIN FOR WET" (posterior jt1 rain) 
          (posteriorMarginal exampleG [winter,sprinkler,wet,road] [rain]  [wet =: True])
diff --git a/Bayes/Test/ReferencePatterns.hs b/Bayes/Test/ReferencePatterns.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Test/ReferencePatterns.hs
@@ -0,0 +1,122 @@
+{- | A comparison of factor elimination with reference values generated with another bayesian network software
+
+It is a non regression test. The test patterns are not provided with this package.
+So, those tests are disabled by default in the hackage version.
+
+-}
+module Bayes.Test.ReferencePatterns(
+#ifdef LOCAL
+   compareAsiaReference
+ , compareCancerReference
+ , comparePokerReference
+ , compareFarmReference
+#endif
+ ) where
+
+import Test.HUnit.Base(assertBool)
+import Data.Maybe(fromJust)
+import qualified Data.Map as Map
+import Bayes.Factor
+import Bayes
+import Bayes.FactorElimination
+import Bayes.Examples(anyExample)
+import Bayes.FactorElimination.JTree(root)
+
+value varmap jt s = 
+  let v =  fromJust $ Map.lookup s varmap
+    in 
+    factorToList (fromJust $ posterior jt v) 
+
+testWithRef varmap jt s l = assertBool s $ value varmap jt s ~=~ l
+testWithRefAndPrint varmap jt s l = do
+  let r = value varmap jt s 
+  putStrLn $ "Computed:" ++ show r
+  putStrLn $ "Reference:" ++ show l
+  putStrLn ""
+  assertBool s $ r ~=~ l
+
+-- Check that the float values are equal with an accuracy < 0.01%
+comparePercent :: Double -> Double -> Bool
+comparePercent a b = abs (a-b) < 1e-4
+
+(~=~) a b = and (zipWith comparePercent a b)
+
+#ifdef LOCAL
+compareFarmReference = do 
+  (varmap,g) <- anyExample "studfarm.net"
+  let jt = createJunctionTree nodeComparisonForTriangulation g
+  
+  assertBool "Junction Tree property" $ junctionTreeProperty jt [] (root jt)
+  testWithRef varmap jt "L"  [0.01,0.99]
+  testWithRef varmap jt "Ann"  [0.01,0.99]
+  testWithRef varmap jt "Brian"  [0.01,0.99]
+  testWithRef varmap jt "Cecily"  [0.01,0.99]
+  testWithRef varmap jt "K"  [0.01,0.99]
+  testWithRef varmap jt "Fred"  [0.01,0.99]
+  testWithRef varmap jt "Dorothy"  [0.01,0.99]
+  testWithRef varmap jt "Eric"  [0.01,0.99]
+  testWithRef varmap jt "Gwenn"  [0.01,0.99]
+  testWithRef varmap jt "Henry"  [0.0091,0.9909]
+  testWithRef varmap jt "Irene"  [0.0099,0.9901]
+  testWithRef varmap jt "John"  [0.0004,0.0087,0.9909]
+
+
+comparePokerReference = do 
+  (varmap,g) <- anyExample "poker.net"
+  let jt = createJunctionTree nodeComparisonForTriangulation g
+
+  assertBool "Junction Tree property" $ junctionTreeProperty jt [] (root jt)
+  testWithRef varmap jt "OH0"  [0.1672, 0.0445,0.0635,0.4659,0.1694,0.0494,0.0353,0.0024,0.0024]
+  testWithRef varmap jt "OH1"  [0.0265,0.0170,0.0357,0.4125,0.2633,0.1599,0.0676,0.0098,0.0077]
+  testWithRef varmap jt "OH2"  [0.2472,0.0628,0.2903,0.0258,0.2526,0.0881,0.0212,0.0121]
+  testWithRef varmap jt "SC"  [0.2450,0.7116,0.0435]
+  testWithRef varmap jt "FC"  [0.0895,0.6988,0.0445,0.1672]
+  testWithRef varmap jt "Besthand"  [0.6396,0.3604]
+  testWithRef varmap jt "MH"  [0.1250,0.1250,0.1250,0.1250,0.1250,0.1250,0.1250,0.1250]
+
+
+compareAsiaReference = do 
+  (varmap,g) <- anyExample "asia.net"
+  let jt = createJunctionTree nodeComparisonForTriangulation g
+
+  assertBool "Junction Tree property" $ junctionTreeProperty jt [] (root jt)
+  testWithRef varmap jt "A"  [0.0100, 0.9900]
+  testWithRef varmap jt "S"  [0.5000, 0.5000]
+  testWithRef varmap jt "T"  [0.0104, 0.9896]
+  testWithRef varmap jt "L"  [0.0550, 0.9450]
+  testWithRef varmap jt "B"  [0.4500, 0.5500]
+  testWithRef varmap jt "E"  [0.0648, 0.9352]
+  testWithRef varmap jt "X"  [0.1103, 0.8897]
+  testWithRef varmap jt "D"  [0.4360, 0.5640]
+
+-- | Type defined to set the evidence on the Coma variable
+-- from the cancer network.
+data Coma = Present | Absent deriving(Eq,Enum,Bounded)
+
+compareCancerReference = do 
+  (varmap,g) <- anyExample "cancer.net"
+  let jt = createJunctionTree nodeComparisonForTriangulation g
+
+  assertBool "Junction Tree property" $ junctionTreeProperty jt [] (root jt)
+  testWithRef varmap jt "A"  [0.2000, 0.8000]
+  testWithRef varmap jt "B"  [0.3200, 0.6800]
+  testWithRef varmap jt "C"  [0.0800, 0.9200]
+  testWithRef varmap jt "D"  [0.3200, 0.6800]
+  testWithRef varmap jt "E"  [0.6160, 0.3840]
+
+  let varD = fromJust $ Map.lookup "D" varmap
+  let jt' = changeEvidence [varD =: Present] jt 
+  testWithRef varmap jt' "A"  [0.4250, 0.5750]
+  testWithRef varmap jt' "B"  [0.8000, 0.2000]
+  testWithRef varmap jt' "C"  [0.2000, 0.8000]
+  testWithRef varmap jt' "D"  [1.0000, 0.0000]
+  testWithRef varmap jt' "E"  [0.6400, 0.3600]
+
+  let jt'' = changeEvidence [varD =: Absent] jt'
+  testWithRef varmap jt'' "A"  [0.0941, 0.9059]
+  testWithRef varmap jt'' "B"  [0.0941, 0.9059]
+  testWithRef varmap jt'' "C"  [0.0235, 0.9765]
+  testWithRef varmap jt'' "D"  [0.0000, 1.0000]
+  testWithRef varmap jt'' "E"  [0.6047, 0.3953]
+
+#endif
diff --git a/Bayes/VariableElimination.hs b/Bayes/VariableElimination.hs
--- a/Bayes/VariableElimination.hs
+++ b/Bayes/VariableElimination.hs
@@ -26,12 +26,12 @@
 --debug s a = trace (s  ++ "\n" ++ show a ++ "\n") a
 
 -- | Elimination order
-type EliminationOrder = DVSet
+type EliminationOrder = [DV]
 
 -- | Get all variables from a Bayesian Network
 allVariables :: (Graph g, Factor f) 
              => BayesianNetwork g f 
-             -> DVSet
+             -> [DV]
 allVariables g = 
   let s = allVertexValues g 
       createDV = factorMainVariable 
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.1
+Version:             0.2
 
 -- A short (one-line) description of the package.
 Synopsis:            Inference with Discrete Bayesian Networks
@@ -15,9 +15,7 @@
 -- A longer description of the package.
 Description:  Algorithms for inference with Discrete Bayesian Networks.  
  It is a very preliminary version. It has only been tested on very simple
- examples where it worked. On bigger networks, imported from Hugin files, it was very very very slow.
- So, you can use this software as a toy. Much more work is needed to validate
- and optimize it.     
+ examples where it worked. This 0.2 version is using new faster and cleaner algorithms.
 
 -- URL for the project homepage or repository.
 Homepage:            http://www.alpheccar.org
@@ -48,7 +46,16 @@
 
 data-files: cancer.net
 
+Flag local {
+  Description: Enable local tests by the author. They can only be run with some test patterns not distributed with this package.
+  Default:     False
+}
 
+Flag profile {
+  Description: Build profiling version of the library too.
+  Default:     False
+}
+
 Library
   -- Modules exported by the library.
   Exposed-modules:
@@ -61,11 +68,19 @@
     Bayes.Test.CompareEliminations
     Bayes.Examples
     Bayes.Examples.Tutorial
+    Bayes.Test.ReferencePatterns
   other-modules:
     Paths_hbayes
     Bayes.ImportExport.HuginNet.Splitting
+    Bayes.PrivateTypes
+    Bayes.FactorElimination.JTree
 
   GHC-Options: -O2 -funbox-strict-fields
+  Extensions: CPP
+  if flag(local)
+    cpp-options: -DLOCAL
+  if flag(profile)
+    ghc-options: -auto-all
 
   
   -- Packages needed in order to build this package.