diff --git a/Bayes.hs b/Bayes.hs
--- a/Bayes.hs
+++ b/Bayes.hs
@@ -2,6 +2,7 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 {- | Discrete Bayesian Network Library.
 
 It is a very preliminary version. It has only been tested on very simple
@@ -21,6 +22,7 @@
   , DirectedGraph(..)
   , FoldableWithVertex(..)
   , NamedGraph(..)
+  , Distribution(..)
   -- ** Graph Monad
   , GraphMonad
   , GMState(..)
@@ -49,13 +51,26 @@
   , runBN 
   , evalBN
   , execBN
+  -- ** Variable creation
   , variable
   , unamedVariable
   , variableWithSize
+  , tdv
+  , t
+  -- ** Creation of conditional probability tables
   , cpt
   , proba
-  , t
   , (~~)
+  , softEvidence
+  , se
+  -- ** Creation of truth tables
+  , logical 
+  , (.==.)
+  , (.!.)
+  , (.|.)
+  , (.&.)
+  -- ** Noisy OR
+  , noisyOR
   -- * Testing
   , testEdgeRemoval_prop
   , testVertexRemoval_prop
@@ -73,11 +88,13 @@
 import qualified Data.Traversable as T 
 import Control.Applicative 
 import qualified Data.Set as Set
+import qualified Data.List as L(find)
 
-import Test.QuickCheck
+import Test.QuickCheck hiding ((.&.),Testable)
 import Test.QuickCheck.Arbitrary
 import Data.List(sort,intercalate,nub)
 import Bayes.PrivateTypes hiding(isEmpty)
+import GHC.Float(float2Double)
 
 --import Debug.Trace
 --debug a = trace (show a) a
@@ -324,6 +341,7 @@
 
 -- | Create an edge description
 edge :: Vertex -> Vertex -> Edge 
+{-# INLINE edge #-}
 edge a b = Edge a b
 
 -- | Endpoints of an edge
@@ -411,6 +429,9 @@
 -- | Used to prevent adding duplicates to a graph
 noRedundancy new old = old
 
+instance FactorContainer (SimpleGraph local edge) where 
+   changeFactor = changeFactorInFunctor 
+
 instance Functor (SimpleGraph local edge) where 
   fmap f (SP em vm nm) = SP em (IM.map (\(l,d) -> (l, f d)) vm) nm
 
@@ -530,6 +551,8 @@
 
 _allEdgeValues (SP em _ _) = M.elems em
 
+_isLinkedWithAnEdge :: SimpleGraph n e v -> Vertex -> Vertex -> Bool 
+{-# INLINE _isLinkedWithAnEdge #-}
 _isLinkedWithAnEdge (SP em _ _) va vb = M.member (edge va vb) em || M.member (edge vb va) em
 
 _someVertex (SP _ vm _) = 
@@ -548,13 +571,13 @@
 
 _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)
+_removeVertex v@(Vertex vertex) g@(SP _ vm _)  = maybe g removeVertexWithValue $! (IM.lookup vertex vm)
   where
     removeVertexWithValue (n,_) = let g' = foldr _removeEdge g (ingoingNeighbors n)
                                       SP em vm' nm' = foldr _removeEdge g' (outgoingNeighbors n)
                                   in 
                                   SP em (IM.delete vertex vm') nm'
-_vertexValue g@(SP _ vm _) (Vertex i) = maybe Nothing (Just . extractValue) (IM.lookup i vm)
+_vertexValue g@(SP _ vm _) (Vertex i) = maybe Nothing (Just . extractValue) $! (IM.lookup i vm)
   where
     extractValue (_,d) = d
 
@@ -588,16 +611,32 @@
       else
         Nothing
 
+_edgeValue :: SimpleGraph n e v -> Edge -> Maybe e 
+{-# INLINE _edgeValue #-}
 _edgeValue (SP em _ _) e = do
      v <- M.lookup e em
      return v
 
+addEdgeReference :: NeighborhoodStructure local 
+                 => Edge 
+                 -> IM.IntMap (local, vertexdata) 
+                 -> Vertex 
+                 -> Vertex 
+                 -> IM.IntMap (local, vertexdata)
+{-# INLINE addEdgeReference #-}
+addEdgeReference newEdge vm (Vertex vsi) (Vertex vei) = id $! IM.adjust addi vei $! (IM.adjust addo vsi vm)
+ where
+   addi (n,v) = (addIngoingEdge newEdge n,v)  
+   addo (n,v) = (addOutgoingEdge newEdge n,v)  
+
+_addEdge :: (NeighborhoodStructure n,Graph (SimpleGraph n)) => Edge -> e -> SimpleGraph n e v -> SimpleGraph n e v 
+{-# INLINE _addEdge #-}
 _addEdge newEdge@(Edge vs ve) value g@(SP em vm nm)   = 
   if testEdgeExistence g em vs ve 
     then 
       g
     else
-      SP (M.insert newEdge value em) (addEdgeReference vm vs ve) nm
+      SP (M.insert newEdge value em) (addEdgeReference newEdge vm vs ve) nm
   where
     testEdgeExistence g em va vb = 
       if (oriented g)
@@ -605,9 +644,7 @@
           M.member (Edge va vb) em
         else 
           M.member (Edge va vb) em || M.member (Edge vb va) em 
-    addEdgeReference vm (Vertex vsi) (Vertex vei) = IM.adjust addi vei (IM.adjust addo vsi vm)
-    addi (n,v) = (addIngoingEdge newEdge n,v)  
-    addo (n,v) = (addOutgoingEdge newEdge n,v)  
+    
 
 instance UndirectedGraph UndirectedSG where
   edges g@(SP _ vm _) v@(Vertex vi) =
@@ -778,26 +815,42 @@
   whenJust newGraph $ \nvm -> do
      put $! (aux, nvm)
 
+-- | A distribution which can be used to create a factor
+class Distribution d where
+  -- | Create a factor from variables and a distributions for those variables
+  createFactor :: Factor f => [DV] -> d -> Maybe f
+
+instance Real a => Distribution [a] where 
+  createFactor dvs l = factorWithVariables dvs (map realToFrac l)
+
+setCpt :: (DirectedGraph g, Distribution d, Factor f) 
+       => g () (MaybeBNode f )
+       -> d 
+       -> Vertex 
+       -> Maybe DV 
+       -> MaybeBNode f 
+       -> BNMonad g f () 
+setCpt g _ _ _ (InitializedBNode _ _ _) = return ()
+setCpt g l v current (UninitializedBNode s dim) = do 
+  let vertices = map (fromJust . startVertex g) . fromJust . ingoing g $ v
+  fv <- mapM factorVariable vertices
+  let cpt = createFactor (map fromJust (current:fv)) l
+      newValue r = InitializedBNode s dim r
+  maybe (return ()) (setBayesianNode v . newValue) cpt
+
 -- | Initialize the values of a factor
-(~~) :: (DirectedGraph g, Factor f) 
-     => BNMonad g f DV -- ^ Discrete variable in the graph
-     -> [Double] -- ^ List of values
+(~~) :: (DirectedGraph g, Factor f, Distribution d, BayesianDiscreteVariable v) 
+     => BNMonad g f v -- ^ Discrete variable in the graph
+     -> d -- ^ List of values
      -> BNMonad g f ()
 (~~) mv l = do 
-  (DV v _) <- mv -- This is updating the state and so the graph
+  (DV v _) <- mv >>= return . dv -- This is updating the state and so the graph
   g <- gets snd
   current <- factorVariable v
   mvalue <- getBayesianNode v
-  maybe (return ()) (setCpt g v current) mvalue
- where
-  setCpt g _ _ (InitializedBNode _ _ _) = return ()
-  setCpt g v current (UninitializedBNode s dim) = do 
-    let vertices = map (fromJust . startVertex g) . fromJust . ingoing g $ v
-    fv <- mapM factorVariable vertices
-    let cpt = factorWithVariables (map fromJust (current:fv)) l
-        newValue r = InitializedBNode s dim r
-    maybe (return ()) (setBayesianNode v . newValue) cpt
+  maybe (return ()) (setCpt g l v current) mvalue
 
+
     
 minBoundForEnum :: Bounded a => a -> a
 minBoundForEnum _ = minBound
@@ -839,7 +892,7 @@
 graphNode :: NamedGraph g => String -> f -> GraphMonad g e f Vertex 
 graphNode vertexName initValue = do
   ((namemap,_),_) <- get
-  maybe (getNewEmptyVariable (Just vertexName) initValue) returnVertex (M.lookup vertexName namemap)
+  maybe (getNewEmptyVariable (Just vertexName) initValue) returnVertex $! (M.lookup vertexName namemap)
    where
     returnVertex i = return (Vertex i)
 
@@ -857,19 +910,19 @@
 _initializeNewVariable :: (Enum a, Bounded a, NamedGraph g)
                        => Vertex 
                        -> a 
-                       -> BNMonad g f DV
+                       -> BNMonad g f (TDV a)
 _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) 
+     UninitializedBNode s d -> return (tdv $ DV va d)
+     InitializedBNode _ d _ -> return (tdv $ DV va d) 
 
 -- | Create a new unamed variable
 unamedVariable :: (Enum a, Bounded a, NamedGraph g)
                => a -- ^ Variable bounds 
-               -> BNMonad g f DV 
+               -> BNMonad g f (TDV a)
 unamedVariable e = do 
   va <- getNewEmptyVariable Nothing (UninitializedBNode "unamed" 0)
   _initializeNewVariable va e
@@ -878,7 +931,7 @@
 variable :: (Enum a, Bounded a, NamedGraph g) 
         => String -- ^ Variable name
         -> a -- ^ Variable bounds
-        -> BNMonad g f DV
+        -> BNMonad g f (TDV a)
 variable name e = do
   va <- addVariableIfNotFound name
   _initializeNewVariable va e
@@ -890,6 +943,22 @@
         -> BNMonad g f DV
 variableWithSize name e = do
   va <- addVariableIfNotFound name
+  _initializeNewVariableWithSize va e
+
+-- | Define a Bayesian variable (name and bounds)
+unNamedVariableWithSize :: NamedGraph g
+                        => Int -- ^ Variable size
+                        -> BNMonad g f DV
+unNamedVariableWithSize e = do
+  va <- getNewEmptyVariable Nothing (UninitializedBNode "unamed" 0)
+  _initializeNewVariableWithSize va e
+
+-- | Initialize a new variable with size
+_initializeNewVariableWithSize :: NamedGraph g
+                               => Vertex -- ^ Variable name
+                               -> Int -- ^ Variable size
+                               -> BNMonad g f DV
+_initializeNewVariableWithSize va e = do
   setVariableBoundWithSize va 0 (e-1)
   maybeValue <- getBayesianNode va 
   setBayesianNode va (fromJust maybeValue)
@@ -901,17 +970,140 @@
 -- Variables are ordered like
 -- FFF FFT FTF FTT TFF TFT TTF TTT
 -- and same for other enumeration keeping enumeration order
-cpt :: DirectedGraph g => DV -> [DV] -> BNMonad g f DV
+cpt :: (DirectedGraph g , BayesianDiscreteVariable v,BayesianDiscreteVariable vb) => v -> [vb] -> BNMonad g f v
 cpt node conditions = do
-  mapM_ (node <--) (reverse conditions)
+  mapM_ ((dv node) <--) (reverse (map dv conditions))
   return node
 
 -- | Define proba for a variable
 -- Values are ordered like
 -- FFF FFT FTF FTT TFF TFT TTF TTT
 -- and same for other enumeration keeping enumeration order
-proba :: DirectedGraph g => DV -> BNMonad g f DV
-proba node = cpt node []
+proba :: (DirectedGraph g, BayesianDiscreteVariable v) => v -> BNMonad g f v
+proba node = cpt node ([] :: [DV])
+
+-- | Create an auxiliairy node to force soft evidence
+softEvidence :: (NamedGraph g, DirectedGraph g, Factor f) 
+             => TDV Bool -- ^ Variable on which we want to define Soft evidence
+             -> BNMonad g f (TDV Bool) -- ^ Return a soft evidence node (for the factor encoding the soft evidence values)
+             -- and an hard evidence node to activate the soft evidence observation
+softEvidence d = do 
+  se <- unNamedVariableWithSize (dimension d) 
+  --seEnabled <- unNamedVariableWithSize (dimension d) 
+
+  cpt se [dv d] ~~ [1.0,0.0,1.0,0.0]
+  --cpt seEnabled [dv se] ~~ [1.0,0.0,0.0,1.0] -- No info about the observation of the soft evidence node
+  return (tdv se) 
+
+-- | Soft evidence factor
+se :: Factor f 
+   => TDV s -- ^ Soft evidence node
+   -> TDV s -- ^ Node on which the soft evidence is imposed
+   -> Double -- ^ Soft evidence (probability of right detection)
+   -> Maybe f
+se s orgNode p = factorWithVariables [dv s,dv orgNode] [p,1-p,1-p,p]
+
+{-
+
+Helper functions to create logical distributions  
+
+-}
+
+data LE = LETest DVI
+        | LEAnd LE LE 
+        | LEOr LE LE 
+        | LENot LE 
+        deriving(Eq)
+
+-- | Generate the variables used in the expression
+varsFromLE :: LE -> [DV]
+varsFromLE le = nub $ _getVars le 
+ where 
+  _getVars  (LETest dvi) = [dv dvi] 
+  _getVars (LEAnd a b) = _getVars a ++ _getVars b
+  _getVars (LEOr a b) = _getVars a ++ _getVars b
+  _getVars (LENot a) = _getVars a
+
+boolValue :: Maybe Bool -> Bool 
+boolValue (Just True) = True 
+boolValue _ = False
+
+-- | Generate values for the LE
+functionFromLE :: LE -> ([DVI] -> Bool)
+functionFromLE (LETest dvi) = \i -> boolValue $ do 
+  var <- L.find (== dvi) i
+  return (instantiationValue dvi == instantiationValue var)
+functionFromLE (LENot l) = \i -> not (functionFromLE l i)
+functionFromLE (LEAnd la lb) = \i -> (functionFromLE la i) && (functionFromLE lb i)
+functionFromLE (LEOr la lb) = \i -> (functionFromLE la i) || (functionFromLE lb i)
+
+class Testable d v where 
+  -- | Create a variable instantiation using values from
+  -- an enumeration
+  (.==.) :: d -> v -> LE 
+
+instance Instantiable d v => Testable d v where 
+  (.==.) a b = LETest (a =: b)
+
+infixl 8 .==.
+infixl 6 .&.
+infixl 5 .|.
+
+(.|.) :: LE -> LE -> LE
+(.|.)  = LEOr 
+
+(.&.) :: LE -> LE -> LE
+(.&.) = LEAnd
+
+(.!.) :: LE -> LE
+(.!.) = LENot
+
+logical :: (Factor f, DirectedGraph g) => TDV Bool -> LE -> BNMonad g f () 
+logical dv l = 
+  let theVars = varsFromLE l
+      logicalF = functionFromLE l 
+      probaVal True = 1.0 :: Double
+      probaVal False = 0.0 :: Double
+      valuesF = [probaVal (logicalF i == False) | i <-forAllInstantiations (DVSet theVars)]
+      valuesT = [probaVal (logicalF i == True) | i <-forAllInstantiations (DVSet theVars)]
+
+  in 
+  cpt dv theVars ~~ (valuesF ++ valuesT)
+
+{-
+
+Noisy OR
+
+-}
+
+-- | Noisy AND. Variable A is passed with probability 1-p
+noisyAND :: (DirectedGraph g, Factor f, NamedGraph g) => TDV Bool -> Double -> BNMonad g f (TDV Bool) 
+noisyAND a p = do 
+    na <- unamedVariable (t::Bool)
+    cpt na [dv a] ~~ [1-p,p,p,1-p]
+    return na 
+
+-- | OR Gate
+orG :: (DirectedGraph g, Factor f, NamedGraph g) => TDV Bool -> TDV Bool -> BNMonad g f (TDV Bool)
+orG a b = do 
+    no <- unamedVariable (t::Bool)
+    logical no ((a .==. True) .|. (b .==. True))
+    return no 
+
+-- | Noisy OR. The Noisy-OR with leak can be implemented by using the
+-- standard Noisy-OR and a leak variable.
+noisyOR :: (DirectedGraph g, Factor f, NamedGraph g) 
+        => [(TDV Bool,Double)] -- ^ Variables and probability of no influence
+        -> BNMonad g f (TDV Bool) 
+noisyOR l = do 
+    a <- mapM (\(a,p) -> noisyAND a p) l
+    foldM orG (head a) (tail a)
+
+{-
+ 
+Graph creation from the Monad.  
+
+-}
 
 
 runGraph :: Graph g => GraphMonad g e f a -> (a,g e f)
diff --git a/Bayes/Examples.hs b/Bayes/Examples.hs
--- a/Bayes/Examples.hs
+++ b/Bayes/Examples.hs
@@ -29,7 +29,7 @@
 construction because those variables are used for making inferences.
 
 @
-example :: ('DVSet','SBN' 'CPT')
+example :: (['TDV' Bool],'SBN' 'CPT')
 example = 'runBN' $ do 
     winter <- 'variable' \"winter\" (t :: Bool)
     sprinkler <- 'variable' \"sprinkler\" (t :: Bool) 
@@ -45,6 +45,66 @@
     return [winter,sprinkler,rain,wet,road]
 @
 
+By default, all variables are typed ('TDV' Bool). 'TDV' means Typed Discrete Variable.
+
+In case you are mixing several types, you'll need to remove the type
+to build the 'cpt' since the list can't be heterogeneous. Just use 'dv' for this. It will convert the variable into the 
+type 'DV' of untyped discrete variable.
+
+/Creating truth tables/
+
+In practise, it is easy to compute the posterior of a variable because it is always possible
+to find a cluster containing the variable in the junction tree. But, it is more difficult
+to compute the posterior of a logical assertion or just a conjunction of assertions.
+
+If a query is likely to be done often, then it may be a good idea to add a new node 
+to the Bayesian network to represent this query. So, some functions to create truth tables are provided.
+
+@
+exampleLogical :: (['TDV' Bool], 'SBN' 'CPT')
+exampleLogical = 'runBN' $ do 
+    a <- 'variable' \"a\" (t :: Bool)
+    b <- 'variable' \"b\" (t :: Bool)
+    notV <- 'variable' \"notV\" (t :: Bool)
+    andV <- 'variable' \"andV\" (t :: Bool)
+    orV <- 'variable' \"orV\" (t :: Bool)
+    let ta = a '.==.' True 
+        tb = b '.==.' True
+    'logical' notV (('.!.') ta)
+    'logical' andV (ta '.&.' tb)
+    'logical' orV (ta '.|.' tb)
+    return $ [a,b,notV,andV,orV]
+@
+
+In the previous example, we force a type on the discrete variables 'DV' to avoid futur errors
+in the instantiations. It is done through the 'tdv' function.
+
+But, it is also possible to use the untyped variables and write:
+
+@
+    'logical' andV ((a '.==.' True) '.&.' (b '.==.' True))
+@
+
+The goal of a Bayesian network is to factorize a big probability table because otherwise the algorithms
+can't process it. So, of course it is not a good idea to represent a complex logical assertion with a huge
+probability table. So, the 'logical' keyword should only be used to build small tables.
+
+If you need to encode a complex logical assertion, use 'logical' several times to build a network representing
+the assertion instead of building just one node to represent it.
+
+/Noisy OR/
+
+The Noisy OR is a combination of logical tables (OR) and conditional probability tables which is often used
+during modeling to avoid generating big conditional probability tables.
+
+It is easy to use:
+
+@
+    no <- 'noisyOR' [(a,0.1),(b,0.2),(c,0.3)] 
+@
+
+Each probability is the probability that a given variable has no effect (so is inhibited in the OR).
+
 /Importing a network from a Hugin file/
 
 The 'exampleImport' function can be used to import a file in Hugin format.
@@ -69,6 +129,8 @@
 module Bayes.Examples(
    example
  , exampleJunction
+ , exampleWithFactorChange
+ , exampleSoftEvidence
 #ifndef LOCAL
  , exampleImport
 #endif
@@ -77,6 +139,7 @@
  , examplePoker
  , exampleFarm
  , examplePerso
+ , exampleLogical
  , testJunction
  , anyExample
  ) where 
@@ -88,7 +151,9 @@
 import qualified Data.Map as Map
 import System.Directory(getHomeDirectory)
 import System.FilePath((</>))
+import Bayes.Factor.CPT
 
+
 #ifndef LOCAL
 import Paths_hbayes
 
@@ -140,21 +205,68 @@
     genericExample $ h </> "Dropbox/bayes_examples/mytest.net"
 
 
+-- | Example of soft evidence use
+exampleSoftEvidence :: ((TDV Bool,TDV Bool),SBN CPT)
+exampleSoftEvidence = runBN $ do
+  a <- variable "a" (t :: Bool)
+  proba a ~~ [0.5,0.5]
+  se <- softEvidence a 
+  return (a,se) 
+
 -- | Standard example found in many books about Bayesian Networks.
-example :: ([DV],SBN CPT)
+example :: ([TDV Bool],SBN CPT)
 example = runBN $ do 
     winter <- variable "winter" (t :: Bool)
     sprinkler <- variable "sprinkler" (t :: Bool) 
     wet <- variable "wet grass" (t :: Bool) 
     rain <- variable "rain" (t :: Bool) 
     road <- variable "slippery road" (t :: Bool) 
+    roadandrain <- variable "rain and slippery road" (t :: Bool)
 
     proba winter ~~ [0.4,0.6]
     cpt sprinkler [winter] ~~ [0.25,0.8,0.75,0.2]
     cpt rain [winter] ~~ [0.9,0.2,0.1,0.8]
     cpt wet [sprinkler,rain] ~~ [1,0.2,0.1,0.05,0,0.8,0.9,0.95]
     cpt road [rain] ~~ [1,0.3,0,0.7]
-    return [winter,sprinkler,rain,wet,road]
+
+    logical roadandrain ((rain .==. True) .&. (road .==. True))
+    return $ [winter,sprinkler,rain,wet,road,roadandrain]
+
+-- | Standard example but with a wrong factor that is changed
+-- in the tests using factor replacement functions
+exampleWithFactorChange :: ([TDV Bool],SBN CPT)
+exampleWithFactorChange = runBN $ do 
+    winter <- variable "winter" (t :: Bool)
+    sprinkler <- variable "sprinkler" (t :: Bool) 
+    wet <- variable "wet grass" (t :: Bool) 
+    rain <- variable "rain" (t :: Bool) 
+    road <- variable "slippery road" (t :: Bool) 
+    roadandrain <- variable "rain and slippery road" (t :: Bool)
+
+    proba winter ~~ [0.4,0.6]
+    cpt sprinkler [winter] ~~ [0.25,0.8,0.75,0.2]
+    cpt rain [winter] ~~ [0.9,0.2,0.1,0.8]
+    cpt wet [sprinkler,rain] ~~ [1,1,1,1,1,1,1,1]
+    cpt road [rain] ~~ [1,0.3,0,0.7]
+
+    logical roadandrain ((rain .==. True) .&. (road .==. True))
+    return $ [winter,sprinkler,rain,wet,road,roadandrain]
+
+
+exampleLogical :: ([TDV Bool], SBN CPT)
+exampleLogical = runBN $ do 
+    a <- variable "a" (t :: Bool)
+    b <- variable "b" (t :: Bool)
+    notV <- variable "notV" (t :: Bool)
+    andV <- variable "andV" (t :: Bool)
+    orV <- variable "orV" (t :: Bool)
+    let ta = a .==. True 
+        tb = b .==. True
+    logical notV ((.!.) ta)
+    logical andV (ta .&. tb)
+    logical orV (ta .|. tb)
+    return $ [a,b,notV,andV,orV]
+
 
 testJunction  :: DirectedSG () Vertex
 testJunction = execGraph $ do
diff --git a/Bayes/Examples/Tutorial.hs b/Bayes/Examples/Tutorial.hs
--- a/Bayes/Examples/Tutorial.hs
+++ b/Bayes/Examples/Tutorial.hs
@@ -12,7 +12,7 @@
 First, the 'example' is loaded to make its variables and its bayesian network available:
 
 @
-    let ([winter,sprinkler,rain,wet,road],exampleG) = 'example'
+    let ([winter,sprinkler,rain,wet,road],exampleG) = example
 @
 
 Then, we compute a prior marginal. Prior means that no evidence is used. A bayesian
@@ -76,10 +76,37 @@
 To set evidence, you need to update the junction tree with new evidence:
 
 @
-    let jt' = 'updateEvidence' [wet '=:'' True] jt 
+    let jt' = 'updateEvidence' [wet '=:' True] jt 
     'posterior' jt' rain
 @
 
+If you want to compute the posterior for a combination of variables, you have two possibilities : either going back to the
+variable elimination methods. Or, introduce new nodes in the network to represent the query.
+
+It is easily done through the new 'logical' function when building the Bayesian graph.
+
+Once you have a node to represent a complex query, you can use it to compute a posterior. For instance, in the rain example,
+there is a new variable:
+
+@
+    roadandrain <- 'variable' \"rain and slippery road\" (t :: Bool)
+    'logical' roadandrain ((rain '.==.' True) '.&.' (road '.==.' True))
+@
+
+This variable is representing the assertion : rain True AND slippery road True. This variable can be used
+ to answer different queries, like for instance:
+
+@
+    let jt4 = 'changeEvidence' [wet '=:' True] jt 
+    print \"Posterior Marginal : probability of rain and road slippery if grass wet\"
+    let m = 'posterior' jt4 roadandrain
+    print m
+--
+    let jt5 = 'changeEvidence' [wet '=:' True, sprinkler '=:' False] jt 
+    print \"Posterior Marginal : probability of rain and road slippery if grass wet and srinkler not used\"
+    let m = 'posterior' jt5 roadandrain
+@
+
 /Inferences with an imported network/
 
 There is a slight additional difficulty with imported networks : you need
@@ -95,53 +122,126 @@
     data Coma = Present | Absent deriving(Eq,Enum,Bounded)
 @
 
-and check that 'Present' is corresponding to the level 0 in the importd network.
+and check that 'Present' is corresponding to the level 0 in the imported network.
 
 Once this datatype is created, you can easily use the cancer network. First we load
 the network and import the discrete variables of type 'DV' from the names of the nodes in the
-network (not the label of the nodes)
+network (not the label of the nodes).
 
 @
     print \"CANCER NETWORK\"
     (varmap,cancer) <- 'exampleImport'
     print cancer
-    let [varA,varB,varC,varD,varE] = fromJust $ mapM (flip Map.lookup varmap) ["A","B","C","D","E"]
+    let [varA,varB,varC,varE] = fromJust $ mapM (flip Map.lookup varmap) [\"A\",\"B\",\"C\",\"E\"]
 @
 
+To avoid any errors with the future queries, some imported variables can be transformed into typed variables:
+
+@
+    varD = 'tdv' (fromJust $ Map.lookup \"D\" varmap) :: 'TDV' Coma
+@
+
 Once the variables are available, you can create the junction tree and start making inferences:
 
 @
     let jtcancer = 'createJunctionTree' 'nodeComparisonForTriangulation' cancer
 --
-    mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer $ x)) [varA,varB,varC,varD,varE]
+    mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer $ x)) [varA,varB,varC,varE]
 --
     print \"UPDATED EVIDENCE\"
     let jtcancer' = 'updateEvidence' [varD '=:' Present] jtcancer 
-    mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer' $ x)) [varA,varB,varC,varD,varE]
+    mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer' $ x)) [varA,varB,varC,varE]
 @
 
+The '=:' operator will check that the assignment is type compatible because varD is a typed discrete variable of type 'TDV' Coma.
+
+/MPE inferences/
+
+It is possible to compute the Most Probable Explanation for a set of observation. The syntax is very similar to the
+posterior computation with variable elimination:
+
+@
+    let m = 'mpe' exampleG [wet,road] [winter,sprinkler,rain,roadandrain] [wet '=:' True, road '=:' True]
+@
+
+The first list of variables (which should containg the evidence variables) is summed out.
+The second list of variables is used to maximize the probability.
+Both lists should contain all variables of the Bayesian network and are defining an elimination order.
+
+The result of the mpe functions is a list of instantiations. The result is easier to read when the type information is
+reintroduced. It can be done with the 'tdvi' function:
+
+@
+    let typedResult = map (map 'tdvi') m :: [[('TDV' Bool,Bool)]]
+@
+
+In this example, all variables are boolean ones.
+
+/ Soft Evidence /
+
+Soft evidence is more complex to handle since new node have to be added to the graph.
+And the node factor has to be changed when the node evidence is changed.
+
+Here is how you could do it. First you load an example graph containg a soft evidence node created with 'softEvidence'.
+
+@
+inferencesWithSoftEvidence = do 
+    let ((a,seNode),exampleG) = 'exampleSoftEvidence' 
+@
+
+Then, you create the junction tree as usual and force an hard evidence on the soft evidence node.
+
+@
+        jt = 'createJunctionTree' 'nodeComparisonForTriangulation' exampleG
+        jt' = 'changeEvidence' [seNode '=:' True] jt
+@
+
+This junction tree cannot be used because the soft evidence node created in 'exampleSoftEvidence' has a 
+probability table which is meaningless. You need to update the probability table for a given soft evidence.
+You create a new factor for this:
+
+@
+        theNewFactor x = fromJust $ 'se' seNode a x -- x % success for the sensor
+@
+
+This new factor, can then be used to do inference with different soft evidences.
+
+@
+    print \"Sensor 90%\"
+    print $ posterior ('changeFactor' (theNewFactor 0.9) jt') a
+--
+    print \"Sensor 50%\"
+    print $ posterior ('changeFactor' (theNewFactor 0.5) jt') a
+--
+    print \"Sensor 10%\"
+    print $ posterior ('changeFactor' (theNewFactor 0.1) jt') a
+@
+
 -}
 module Bayes.Examples.Tutorial(
     -- * Tests with the standard network 
       inferencesOnStandardNetwork
+    , mpeStandardNetwork
+    , inferencesWithSoftEvidence
 #ifndef LOCAL
     -- * Tests with the cancer network
     , inferencesOnCancerNetwork
+    , Coma(..)
 #endif
 #ifdef LOCAL
     , miscDiabete
 #endif
-    , Coma(..)
     , miscTest
+    , logicalTest
 	) where 
 
 import Bayes.Factor
 import Bayes
 import Bayes.VariableElimination
 #ifndef LOCAL
-import Bayes.Examples(example, exampleJunction,exampleImport,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)
+import Bayes.Examples(example, exampleLogical,exampleSoftEvidence,exampleJunction,exampleImport,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)
 #else 
-import Bayes.Examples(example, exampleJunction,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)
+import Bayes.Examples(example, exampleLogical,exampleSoftEvidence,exampleJunction,exampleDiabete, exampleAsia, examplePoker, exampleFarm,examplePerso,anyExample)
 #endif
 import Bayes.FactorElimination
 import Data.Function(on)
@@ -178,33 +278,59 @@
   mapM_ prior (zip l names)
 
 
+#ifndef LOCAL
+
 -- | Type defined to set the evidence on the Coma variable
 -- from the cancer network.
 data Coma = Present | Absent deriving(Eq,Enum,Bounded)
 
-#ifndef LOCAL
 -- | Inferences with the cancer network
 inferencesOnCancerNetwork = do 
     print "CANCER NETWORK"
     (varmap,cancer) <- exampleImport
     print cancer
-    let [varA,varB,varC,varD,varE] = fromJust $ mapM (flip Map.lookup varmap) ["A","B","C","D","E"]
-    let jtcancer = createJunctionTree nodeComparisonForTriangulation cancer
-
-    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer $ x)) [varA,varB,varC,varD,varE]
+    let [varA,varB,varC,varE] = fromJust $ mapM (flip Map.lookup varmap) ["A","B","C","E"]
+        jtcancer = createJunctionTree nodeComparisonForTriangulation cancer
+        varD = tdv (fromJust $ Map.lookup "D" varmap) :: TDV Coma
+    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer $ x)) [varA,varB,varC,varE]
 
     print "UPDATED EVIDENCE : Coma present"
     let jtcancer' = changeEvidence [varD =: Present] jtcancer 
-    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varD,varE]
+    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varE]
 
     print "UPDATED EVIDENCE : Coma absent"
     let jtcancer' = changeEvidence [varD =: Absent] jtcancer 
-    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varD,varE]
+    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varE]
 #endif 
 
+-- | Display of factors generated by the logical keyword
+logicalTest = do 
+    let ([a,b,notV,andV,orV],g) = exampleLogical
+        fnot = fromJust $ vertexValue g (vertex notV)
+        fand = fromJust $ vertexValue g (vertex andV)
+        for = fromJust $ vertexValue g (vertex orV)
+    print fnot 
+    print fand 
+    print for
+
+-- | Inferences with soft evidence
+inferencesWithSoftEvidence = do 
+    let ((a,seNode),exampleG) = exampleSoftEvidence 
+        jt = createJunctionTree nodeComparisonForTriangulation exampleG
+        theNewFactor x = fromJust $ se seNode a x -- x % success for the sensor
+        jt' = changeEvidence [seNode =: True] jt
+    print "Sensor 90%"
+    print $ posterior (changeFactor (theNewFactor 0.9) jt') a
+
+    print "Sensor 50%"
+    print $ posterior (changeFactor (theNewFactor 0.5) jt') a
+
+    print "Sensor 10%"
+    print $ posterior (changeFactor (theNewFactor 0.1) jt') a
+
 -- | Inferences with the standard network
 inferencesOnStandardNetwork = do
-    let ([winter,sprinkler,rain,wet,road],exampleG) = example
+    let ([winter,sprinkler,rain,wet,road,roadandrain],exampleG) = example
 
     print exampleG
     putStrLn ""
@@ -250,10 +376,37 @@
     putStrLn ""
 
     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
     print m
     putStrLn ""
 
+    let jt4 = changeEvidence [wet =: True] jt 
+    print "Posterior Marginal : probability of rain and road slippery if grass wet"
+    let m = posterior jt4 roadandrain
+    print m
+    putStrLn ""
+
+    let jt5 = changeEvidence [wet =: True, sprinkler =: False] jt 
+    print "Posterior Marginal : probability of rain and road slippery if grass wet and srinkler not used"
+    let m = posterior jt5 roadandrain
+    print m
+    putStrLn ""
+
     return ()
+
+
+-- | Most likely explanation on standard network
+mpeStandardNetwork = do
+    let ([winter,sprinkler,rain,wet,road,roadandrain],exampleG) = example
+
+    print exampleG
+    print "Most likely explanation if grass wet and road slippery"
+    let m = mpe exampleG [wet,road] [winter,sprinkler,rain,roadandrain] [wet =: True, road =: True]
+        typedResult = map (map tdvi) m :: [[(TDV Bool,Bool)]]
+    print typedResult
+    putStrLn ""
+    let m = mpe exampleG [wet,road,roadandrain,winter] [sprinkler,rain] [wet =: True, road =: True]
+        typedResult = map (map tdvi) m :: [[(TDV Bool,Bool)]]
+    print typedResult
+    putStrLn ""
diff --git a/Bayes/Factor.hs b/Bayes/Factor.hs
--- a/Bayes/Factor.hs
+++ b/Bayes/Factor.hs
@@ -1,8 +1,6 @@
 {-# LANGUAGE TypeSynonymInstances #-}
 {-# LANGUAGE FlexibleInstances #-}
-{- | Conditional probability table
-
-Conditional Probability Tables and Probability tables
+{- | Factors
 
 -}
 module Bayes.Factor(
@@ -10,15 +8,21 @@
    Factor(..)
  , isomorphicFactor
  , normedFactor
+ , displayFactorBody
+ , changeFactorInFunctor 
+ , FactorContainer(..)
  -- * Set of variables 
  , Set(..)
  , BayesianDiscreteVariable(..)
  -- * Implementation
  , Vertex(..)
  -- ** Discrete variables and instantiations
- , DV
+ , DV(..)
+ , TDV
  --, DVSet(..)
  , DVI
+ , DVISet
+ , tdvi
  , setDVValue
  , instantiationValue
  , instantiationVariable
@@ -26,36 +30,40 @@
  , (=:)
  , forAllInstantiations
  , factorFromInstantiation
- , changeVariableOrder
- -- ** Factor
- , CPT
- -- * Tests
- , testProductProject_prop
- , testAssocProduct_prop
- , testScale_prop
- , testProjectCommut_prop
- , testScalarProduct_prop
- , testProjectionToScalar_prop
  ) where
 
-import qualified Data.Vector.Unboxed as V
-import Data.Vector.Unboxed((!))
-import Data.Maybe(fromJust,mapMaybe,isJust)
-import qualified Data.List as L
-import Text.PrettyPrint.Boxes hiding((//))
-import Test.QuickCheck
-import Test.QuickCheck.Arbitrary
-import qualified Data.IntMap as IM
+import Data.Maybe(fromJust)
 import Control.Monad
-import System.Random(Random)
-import Data.List(partition)
 import Bayes.PrivateTypes
+import Bayes.Tools
+import qualified Data.Vector.Unboxed as V
+import Text.PrettyPrint.Boxes hiding((//))
 
 --import Debug.Trace
 
 --debug a = trace ("\nDEBUG\n" ++ show a ++ "\n") a
 
+-- | Change factor in a functor (only factor values should have been changed)
+-- It assumes that the variables of a factor are enough to identify it.
+-- If the functor is containing several factors with same set of variables then it
+-- won't give a meaningful result.
+-- So it should be used only on functor derived from a Bayesian Network.
+changeFactorInFunctor :: (Factor f, Functor m) => f -> m f -> m f
+changeFactorInFunctor f g = 
+  let replaceFactor cf | cf `isUsingSameVariablesAs` f = f 
+                       | otherwise = cf
+  in
+  fmap replaceFactor g 
 
+-- | Structure containing factors which can be replaced.
+-- It is making sense when the factors are related to the nodes of a Bayesian
+-- network. 
+class FactorContainer m where 
+   changeFactor :: Factor f => f -> m f -> m f 
+
+instance FactorContainer [] where 
+    changeFactor = changeFactorInFunctor
+
 -- | A vertex associated to another value (variable dimension, variable value ...)
 class LabeledVertex l where
     variableVertex :: l -> Vertex
@@ -63,7 +71,7 @@
 
 -- | Convert a variable instantation to a factor
 -- Useful to create evidence factors
-factorFromInstantiation :: Factor f => DVI Int -> f
+factorFromInstantiation :: Factor f => DVI -> f
 factorFromInstantiation (DVI dv a) = 
     let setValue i = if i == a then 1.0 else 0.0 
     in
@@ -72,7 +80,7 @@
 
 
 
-instance LabeledVertex (DVI a) where
+instance LabeledVertex DVI where
     variableVertex (DVI v _) = variableVertex v
 
 instance LabeledVertex DV where
@@ -111,7 +119,11 @@
     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 -> [DVI Int] -> Double
+    factorValue :: f -> [DVI] -> Double
+
+    -- | String representation of a factor value 
+    factorStringValue :: f -> [DVI] -> String
+
     -- | 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
@@ -130,7 +142,13 @@
 
     -- | Create an evidence factor from an instantiation.
     -- If the instantiation is empty then we get nothing
-    evidenceFrom :: [DVI Int] -> Maybe f
+    evidenceFrom :: [DVI] -> Maybe f
+
+    -- | Test if two factors are coding for the same probability dependence.
+    -- It does not test if the factors are equal (same probabilities) but just
+    -- if they involve the same variables so are linked to the same
+    -- node in the Bayesian network
+    isUsingSameVariablesAs ::  f -> f -> Bool
     
 
     -- | Divide all the factor values
@@ -153,151 +171,15 @@
         in 
         factorProjectOut s' f
 
--- | Change the layout of values in the
--- factor to correspond to a new variable order
--- Used to import external files
-changeVariableOrder :: DVSet s -- ^ Old order
-                    -> DVSet s' -- ^ New order 
-                    -> [Double] -- ^ Old values
-                    -> [Double] -- ^ New values
-changeVariableOrder (DVSet oldOrder) newOrder oldValues =
-    let oldFactor = fromJust $ factorWithVariables oldOrder oldValues :: CPT
-    in
-    [factorValue oldFactor i | i <- forAllInstantiations newOrder]
 
 
--- | Mainly used for conditional probability table like p(A B | C D E) but the normalization to 1
--- is not imposed. And the conditionned variables are not different from the conditionning ones.
--- The dimensions for each variables are listed.
--- The variables on the left or right of the condition bar are not tracked. What's matter is that
--- it is encoding a function of several variables.
--- Marginalization of variables will be computed from the bayesian graph where
--- the knowledge of the dependencies is.
--- So, this same structure is used for a probability too (conditioned on nothing)
-data CPT = CPT {
-           dimensions :: ![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
 
-debugCPT (Scalar d) = do 
-   putStrLn "SCALAR CPT"
-   print d
-   putStrLn ""
 
-debugCPT (CPT d m v) = do 
-    putStrLn "CPT"
-    print d 
-    putStrLn ""
-    print m 
-    putStrLn ""
-    print v
-    putStrLn ""
-{-
 
-CPT can't have same same vertex values but with different sizes.
-But, arbitrary CPT generation will general several vertex with same vertex id
-and different vertex size.
 
-So, we introduce a function mapping a vertex ID to a vertex size. So, vertex size are hard coded
 
--}
 
-quickCheckVertexSize :: Int -> Int
-quickCheckVertexSize 0 = 2
-quickCheckVertexSize 1 = 2
-quickCheckVertexSize 2 = 2
-quickCheckVertexSize _ = 2
 
--- | Generate a random value until this value is not already present in the list
-whileIn :: (Arbitrary a, Eq a) => [a] -> Gen a -> Gen a
-whileIn l m = do 
-    newVal <- m 
-    if newVal `elem` l 
-        then
-            whileIn l m 
-        else 
-            return newVal
-
--- | Generate a random vector of n elements without replacement (no duplicate)
--- May loop if the range is too small !
-generateWithoutReplacement :: (Random a, Arbitrary a, Eq a)  
-                           => Int -- ^ Vector size
-                           -> (a,a) -- ^ Bounds
-                           -> Gen [a]
-generateWithoutReplacement n b | n == 1 = generateSingle b 
-                               | n > 1 = generateMultiple n b 
-                               | otherwise = return []
- where
-   generateSingle b = do 
-       r <- choose b
-       return [r]
-   generateMultiple n b = do 
-       l <- generateWithoutReplacement (n-1) b
-       newelem <- whileIn l $ choose b
-       return (newelem:l)
-
-
-
-instance Arbitrary CPT where
-    arbitrary = do 
-        nbVertex <- choose (1,4) :: Gen Int
-        vertexNumbers <- generateWithoutReplacement nbVertex (0,50)
-        let dimensions = map (\i -> (DV (Vertex i)  (quickCheckVertexSize i))) vertexNumbers
-        let valuelen = product (map dimension dimensions)
-        rndValues <- vectorOf valuelen (choose (0.0,1.0) :: Gen Double)
-        return . fromJust . factorWithVariables dimensions $ rndValues
-
--- | Test product followed by a projection when the factors have no
--- common variables
-
--- | Floating point number comparisons which should take into account
--- all the subtleties of that kind of comparison
-nearlyEqual :: Double -> Double -> Bool
-nearlyEqual a b = 
-    let absA = abs a 
-        absB = abs b 
-        diff = abs (a-b)
-        epsilon = 2e-5
-    in
-    case (a,b) of 
-        (x,y) | x == y -> True -- handle infinities
-              | x*y == 0 -> diff < (epsilon * epsilon)
-              | otherwise -> diff / (absA + absB) < epsilon
-
-testScale_prop :: Double -> CPT -> Bool
-testScale_prop s f = (factorNorm (s `factorScale` f)) `nearlyEqual` (s * (factorNorm f))
-
-testProductProject_prop :: CPT -> CPT -> Property
-testProductProject_prop fa fb = isEmpty ((factorVariables fa) `intersection` (factorVariables fb))  ==> 
-    let r = factorProjectOut (factorVariables fb) (factorProduct [fa,fb])
-        fa' = r `factorDivide` (factorNorm fb)
-    in
-    fa' `isomorphicFactor` fa
-
-testScalarProduct_prop :: Double -> CPT -> Bool 
-testScalarProduct_prop v f = (factorProduct [(Scalar v),f]) `isomorphicFactor` (v `factorScale` f)
-
-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 = 
-    let allVars = factorVariables f 
-    in
-    (factorProjectOut allVars f) `isomorphicFactor` (factorFromScalar (factorNorm f))
-
-testProjectCommut_prop:: CPT -> Property 
-testProjectCommut_prop f = length (factorVariables f) >= 3 ==>
-    let a = take 1 (factorVariables f)
-        b = take 1 . drop 1 $ factorVariables f 
-        commuta = factorProjectOut a (factorProjectOut b f)
-        commutb = factorProjectOut b (factorProjectOut a f)
-    in
-    commuta `isomorphicFactor` commutb
-
 -- | Test equality of two factors taking into account the fact
 -- that the variables may be in a different order.
 -- In case there is a distinction between conditionned variable and
@@ -322,17 +204,17 @@
 
 -}
 -- | Display a variable name and its size
-vname :: Int -> DVI Int -> Box
+vname :: Int -> DVI -> Box
 vname vc i = text $ "v" ++ show vc ++ "=" ++ show (instantiationValue i)
 
-dispFactor :: Factor f => f -> DV -> [DVI Int] -> [DV] -> Box
+dispFactor :: Factor f => f -> DV -> [DVI] -> [DV] -> Box
 dispFactor cpt h c [] = 
     let dstIndexes = allInstantiationsForOneVariable h
         dependentIndexes =  reverse c
         factorValueAtPosition p = 
-            let v = factorValue cpt p
+            let v = factorStringValue cpt p
             in
-            text . show  $ v
+            text v
     in
     vsep 0 center1 . map (factorValueAtPosition . (:dependentIndexes)) $ dstIndexes
 
@@ -341,12 +223,10 @@
     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
-    show c@(CPT [] _ v) = "\nEmpty CPT:\n"
-
-    show c@(CPT d _ v) = 
-        let h@(DV (Vertex vc) _) = head d
+displayFactorBody :: Factor f => f -> String 
+displayFactorBody c = 
+        let d = factorVariables c
+            h@(DV (Vertex vc) _) = head d
             table = dispFactor c h [] (tail d)
             dstIndexes = map head (forAllInstantiations . DVSet $ [h])
             -- In P(A | B ...), the dst column is containing the possible values for the
@@ -355,134 +235,3 @@
             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
-    factorFromScalar v = Scalar v
-    factorDimension f@(CPT _ _ _) = product . map dimension . factorVariables$ f
-    factorDimension _ = 1
-    containsVariable (CPT _ m _) (DV (Vertex i) _)   = IM.member i m
-    containsVariable (Scalar _) _ = False
-    factorWithVariables = createCPTWithDims
-    factorVariables (CPT v _ _) = v
-    factorVariables (Scalar _) = []
-    factorNorm f@(CPT 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@(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@(CPT d _ v) i = 
-        let vars = DVSet d
-            (dv,pos) = instantiationDetails i
-            strides = indexStridesFor vars dv
-        in 
-        v!(indexPosition strides pos)
-    evidenceFrom [] = Nothing 
-    evidenceFrom l = 
-        let (variables,index) = instantiationDetails l
-            DVSet nakedVars = variables
-            setValueForIndex i = if i == index then 1.0 else 0.0 
-        in
-        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
-
-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)
-    in 
-    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)
-
-
--- | Create a CPT given some dimensions and a list of Doubles.
--- Returns nothing is the length are not coherents.
-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
-        p = product (map dimension dims) 
-    in
-    if length values == p
-        then
-            Just $ CPT dims m (V.fromList values)
-        else 
-            Nothing
-
diff --git a/Bayes/Factor/CPT.hs b/Bayes/Factor/CPT.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Factor/CPT.hs
@@ -0,0 +1,141 @@
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{- | Conditional probability table
+
+Conditional Probability Tables and Probability tables
+
+-}
+module Bayes.Factor.CPT(
+   -- * CPT Factor
+     CPT
+   , changeVariableOrder
+   -- * Tests
+   , testProductProject_prop
+   , testAssocProduct_prop
+   , testScale_prop
+   , testProjectCommut_prop
+   , testScalarProduct_prop
+   , testProjectionToScalar_prop
+   , debugCPT
+   ) where
+
+
+import Bayes.Factor 
+import Bayes.Tools
+import Test.QuickCheck
+import Test.QuickCheck.Arbitrary
+import System.Random(Random)
+import Data.Maybe(fromJust,mapMaybe,isJust)
+import Bayes.Factor.PrivateCPT
+import Bayes.PrivateTypes
+
+-- | Soft evidence factor can be used to initialize a factor
+--instance Distribution CPT where 
+--  createFactor dvs f = factorWithVariables dvs (factorToList f)
+
+--- | Change the layout of values in the
+--- factor to correspond to a new variable order
+--- Used to import external files
+changeVariableOrder :: DVSet s -- ^ Old order
+                    -> DVSet s' -- ^ New order 
+                    -> [Double] -- ^ Old values
+                    -> [Double] -- ^ New values
+changeVariableOrder (DVSet oldOrder) newOrder oldValues =
+    let oldFactor = fromJust $ factorWithVariables oldOrder oldValues :: CPT
+    in
+    [factorValue oldFactor i | i <- forAllInstantiations newOrder]
+
+instance Arbitrary CPT where
+    arbitrary = do 
+        theTVars <- arbitrary :: Gen (DVSet s)
+        let theVars = fromDVSet theTVars
+        --nbVertex <- choose (1,4) :: Gen Int
+        --vertexNumbers <- generateWithoutReplacement nbVertex (0,50)
+        --let dimensions = map (\i -> (DV (Vertex i)  (quickCheckVertexSize i))) vertexNumbers
+        let valuelen = product (map dimension theVars)
+        rndValues <- vectorOf valuelen (choose (0.0,1.0) :: Gen Double)
+        return . fromJust . factorWithVariables theVars $ rndValues
+
+-- | Test product followed by a projection when the factors have no
+-- common variables
+
+
+testScale_prop :: Double -> CPT -> Bool
+testScale_prop s f = (factorNorm (s `factorScale` f)) `nearlyEqual` (s * (factorNorm f))
+
+testProductProject_prop :: CPT -> CPT -> Property
+testProductProject_prop fa fb = isEmpty ((factorVariables fa) `intersection` (factorVariables fb))  ==> 
+    let r = factorProjectOut (factorVariables fb) (factorProduct [fa,fb])
+        fa' = r `factorDivide` (factorNorm fb)
+    in
+    fa' `isomorphicFactor` fa
+
+testScalarProduct_prop :: Double -> CPT -> Bool 
+testScalarProduct_prop v f = (factorProduct [factorFromScalar 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 = 
+    let allVars = factorVariables f 
+    in
+    (factorProjectOut allVars f) `isomorphicFactor` (factorFromScalar (factorNorm f))
+
+testProjectCommut_prop:: CPT -> Property 
+testProjectCommut_prop f = length (factorVariables f) >= 3 ==>
+    let a = take 1 (factorVariables f)
+        b = take 1 . drop 1 $ factorVariables f 
+        commuta = factorProjectOut a (factorProjectOut b f)
+        commutb = factorProjectOut b (factorProjectOut a f)
+    in
+    commuta `isomorphicFactor` commutb
+
+
+
+instance Show CPT where
+    show (Scalar v) = "\nScalar Factor:\n" ++ show v
+    show c@(Table [] _ v) = "\nEmpty CPT:\n"
+
+    show c = displayFactorBody c    
+
+
+instance FactorElement Double where 
+    doubleValue = id 
+    mkValue = id
+    scale = (*)
+    multiply = (*)
+    elementUnit = 1.0
+
+
+instance Factor CPT where
+    emptyFactor = _emptyFactor
+    factorVariables = _factorVariables
+    isScalarFactor = _isScalarFactor
+    containsVariable = _containsVariable
+    factorDimension = _factorDimension
+    variablePosition = _variablePosition
+    isUsingSameVariablesAs = _isUsingSameVariablesAs
+
+    factorFromScalar = _factorFromScalar
+    factorWithVariables = _factorWithVariables
+    factorToList = _factorToList
+    factorNorm = _factorNorm
+    factorScale = _factorScale 
+    factorValue = _factorValue
+    factorStringValue f d = show (_factorValue f d)
+    
+    evidenceFrom = _evidenceFrom
+
+    factorProduct = _factorProduct
+    
+    factorProjectOut _ f@(Scalar v) = f
+    factorProjectOut s f= cptFactorProjectOutWith (sum . map fst) s f
+
+
+
+
+
+
diff --git a/Bayes/Factor/MaxCPT.hs b/Bayes/Factor/MaxCPT.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Factor/MaxCPT.hs
@@ -0,0 +1,80 @@
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
+{- | Implementation of Max-product factors for MAP queries
+
+-}
+module Bayes.Factor.MaxCPT(
+	-- * Type
+      MAXCPT 
+    , mpeInstantiations
+	) where 
+
+import Bayes.Factor 
+import Bayes.Factor.PrivateCPT
+import Bayes.PrivateTypes
+import qualified Data.Vector as V
+import qualified Data.IntMap as IM
+import Data.List(maximumBy)
+import Data.Function(on)
+
+--import Debug.Trace 
+
+--debug a = trace (show a ++ "\n") a
+
+instance FactorElement (Double,PossibleInstantiations) where 
+    doubleValue = fst 
+    mkValue x = (x,[])
+    scale s (d,l) = (s*d,l)
+    multiply (da,[]) (db,lb) = (da*db,lb)
+    multiply (da,la) (db,[]) = (da*db,la)
+    multiply (da,la) (db,lb) = (da*db,[xa ++ xb | xa <- la, xb <- lb])
+    elementUnit = (1.0,[])
+
+instantiationProba :: ((Double,PossibleInstantiations), DVISet) -> Double
+instantiationProba (a,b) = fst a 
+
+maximization :: [((Double,PossibleInstantiations), DVISet)] -> (Double,PossibleInstantiations)
+maximization [] = elementUnit 
+maximization l = 
+    let ((d,possible),newInst) = maximumBy (compare `on` instantiationProba) l
+        addTo [] l = l
+        addTo i [] = [i]
+        addTo i l = map (i ++) l
+    in 
+    (d,addTo newInst possible)
+
+instance Factor (MAXCPT) where
+    emptyFactor = _emptyFactor
+    factorVariables = _factorVariables
+    isScalarFactor = _isScalarFactor
+    containsVariable = _containsVariable
+    factorDimension = _factorDimension
+    variablePosition = _variablePosition
+    isUsingSameVariablesAs = _isUsingSameVariablesAs
+
+    factorFromScalar = _factorFromScalar
+    factorWithVariables = _factorWithVariables
+    factorToList = _factorToList
+    factorNorm = _factorNorm
+    factorScale = _factorScale
+    factorValue = _factorValue
+    factorStringValue f d = show (privateFactorValue f d)
+
+    evidenceFrom = _evidenceFrom
+
+    factorProduct = _factorProduct
+    
+    factorProjectOut _ f@(Scalar v) = f
+    factorProjectOut s f = cptFactorProjectOutWith maximization s f
+
+mpeInstantiations :: MAXCPT -> [DVISet]
+mpeInstantiations (Scalar (_,i)) = i 
+mpeInstantiations _ = error "The final MAXCPT factor should be a scalar one"
+
+instance Show MAXCPT where
+    show (Scalar v) = "\nScalar Factor:\n" ++ show v
+    show c@(Table [] _ v) = "\nEmpty CPT:\n"
+
+    show c = displayFactorBody c  
+
+
diff --git a/Bayes/Factor/PrivateCPT.hs b/Bayes/Factor/PrivateCPT.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Factor/PrivateCPT.hs
@@ -0,0 +1,295 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
+{- | Private functions used to implement CPT
+
+-}
+module Bayes.Factor.PrivateCPT(
+   -- * Functions
+     cptFactorProjectOutWith
+   -- ** Private functions
+   , FactorElement(..)
+   , _emptyFactor
+   , _factorVariables
+   , _isScalarFactor
+   , _containsVariable
+   , _factorDimension
+   , _variablePosition
+   , _factorFromScalar
+   , _factorWithVariables
+   , _factorToList
+   , _factorNorm
+   , _factorScale 
+   , _factorValue
+   , _evidenceFrom
+   , _factorProduct
+   , _isUsingSameVariablesAs
+   , CPT
+   , MAXCPT
+   , PrivateCPT(..)
+   , PossibleInstantiations
+   , convertToMaxFactor
+   , debugCPT
+   , privateFactorValue
+  ) where
+
+import qualified Data.Vector.Generic as GV
+
+import qualified Data.Vector.Unboxed as V
+import Bayes.PrivateTypes
+import qualified Data.IntMap as IM
+import Bayes.Factor
+import Data.Vector.Generic((!))
+import Data.Maybe(fromJust,mapMaybe,isJust)
+import qualified Data.List as L
+import Data.List(partition)
+import qualified Data.Vector as NV
+
+
+{-
+ 
+CPT used for normal CPT and MAXCPT
+
+-}
+-- | Mainly used for conditional probability table like p(A B | C D E) but the normalization to 1
+-- is not imposed. And the conditionned variables are not different from the conditionning ones.
+-- The dimensions for each variables are listed.
+-- The variables on the left or right of the condition bar are not tracked. What's matter is that
+-- it is encoding a function of several variables.
+-- Marginalization of variables will be computed from the bayesian graph where
+-- the knowledge of the dependencies is.
+-- So, this same structure is used for a probability too (conditioned on nothing)
+data PrivateCPT v a = Table {
+                      dimensions :: ![DV] -- ^ Dimensions for all variables
+                    , mapping :: !(IM.IntMap Int) -- ^ Mapping from vertex number to position in dimensions
+                    , values :: !(v a) -- ^ Table of values
+                    }
+                    | Scalar !a 
+
+debugCPT (Scalar d) = do 
+   putStrLn "SCALAR CPT"
+   print d
+   putStrLn ""
+
+debugCPT (Table d m v) = do 
+    putStrLn "CPT"
+    print d 
+    putStrLn ""
+    print m 
+    putStrLn ""
+    print v
+    putStrLn ""
+
+-- | Convert a factor to MAXCPT to change the behavior of inference algorithms and use max
+convertToMaxFactor :: CPT -> MAXCPT
+convertToMaxFactor (Scalar v) = Scalar (v,[])
+convertToMaxFactor (Table d m v) =
+  let newValues = NV.fromList . map (\x -> (x,[])) . V.toList $ v
+  in
+  Table d m newValues
+
+type CPT = PrivateCPT (V.Vector) Double
+
+
+
+type PossibleInstantiations = [DVISet]
+
+type MAXCPT = PrivateCPT (NV.Vector) (Double,PossibleInstantiations)
+
+class FactorElement a where 
+    doubleValue :: a -> Double 
+    mkValue :: Double -> a
+    scale :: Double -> a -> a
+    multiply :: a -> a -> a
+    elementUnit :: a
+
+_isUsingSameVariablesAs :: PrivateCPT v a -> PrivateCPT v a -> Bool
+_isUsingSameVariablesAs (Table va _ _) (Table vb _ _) = va == vb
+
+_emptyFactor :: GV.Vector v a => PrivateCPT v a
+_emptyFactor = Table [] IM.empty GV.empty
+
+_factorVariables :: PrivateCPT v a -> [DV]
+_factorVariables (Table v _ _) = v
+_factorVariables (Scalar _) = []
+
+_factorToList :: (FactorElement a, GV.Vector v a) => PrivateCPT v a -> [Double]
+_factorToList (Scalar v) = [doubleValue v]
+_factorToList (Table _ _ v) = map doubleValue $ GV.toList v
+
+_isScalarFactor :: PrivateCPT v a -> Bool 
+_isScalarFactor (Scalar _) = True
+_isScalarFactor _ = False
+
+_factorFromScalar :: FactorElement a => Double -> PrivateCPT v a
+_factorFromScalar v = Scalar (mkValue v) 
+
+_factorDimension :: Factor (PrivateCPT v a) => PrivateCPT v a -> Int
+_factorDimension f@(Table _ _ _) = product . map dimension . factorVariables$ f
+_factorDimension _ = 1
+
+_containsVariable :: PrivateCPT v a -> DV  -> Bool
+_containsVariable (Table _ m _) (DV (Vertex i) _)   = IM.member i m
+_containsVariable (Scalar _) _ = False
+
+_factorWithVariables :: (FactorElement a, GV.Vector v a) => [DV] -> [Double] -> Maybe (PrivateCPT v a)
+_factorWithVariables dv v = createCPTWithDims dv (GV.fromList (map mkValue v))
+
+_factorNorm :: (FactorElement a, GV.Vector v a) => PrivateCPT v a -> Double 
+_factorNorm f@(Table d _ vals) = 
+    let vars = DVSet d
+        strides = indexStrides vars
+    in
+    sum [ doubleValue(vals!(indexPosition strides x)) | x <- indicesForDomain vars]
+_factorNorm (Scalar v) = (doubleValue v)
+
+
+_variablePosition :: PrivateCPT v a -> DV -> Maybe Int
+_variablePosition (Table _ m _) (DV (Vertex i) _) = let r = IM.lookup i m in r `seq` r
+_variablePosition (Scalar _) _ = Nothing
+
+_factorScale :: (FactorElement a, GV.Vector v a, Factor (PrivateCPT v a)) => Double -> PrivateCPT v a -> PrivateCPT v a
+_factorScale s (Scalar v) = Scalar (scale s v)
+_factorScale s f@(Table d _ vals) = 
+    let vars = DVSet d
+        strides = indexStrides vars
+        newValues = GV.force . GV.fromList $ map (scale s) [ vals!(indexPosition strides x) | x <- indicesForDomain vars]
+    in 
+    fromJust $ createCPTWithDims (factorVariables f) newValues
+
+privateFactorValue :: (FactorElement a, GV.Vector v a) => PrivateCPT v a  -> [DVI] -> a
+{-# INLINE privateFactorValue #-}
+privateFactorValue (Scalar v) _ = v 
+privateFactorValue f@(Table d _ v) i = 
+    let vars = DVSet d
+        (dv,pos) = instantiationDetails i
+        strides = indexStridesFor vars dv
+    in 
+    v!(indexPosition strides pos)
+
+_evidenceFrom :: (FactorElement a, GV.Vector v a) => [DVI] -> Maybe (PrivateCPT v a)
+_evidenceFrom [] = Nothing 
+_evidenceFrom l = 
+    let (variables,index) = instantiationDetails l
+        DVSet nakedVars = variables
+        setValueForIndex i = if i == index then mkValue 1.0 else mkValue 0.0 
+        values = GV.force . GV.fromList $ map setValueForIndex $ indicesForDomain variables
+    in
+    createCPTWithDims nakedVars values
+
+_factorValue :: (FactorElement a, GV.Vector v a) => PrivateCPT v a  -> [DVI] -> Double
+{-# INLINE _factorValue #-}
+_factorValue f d = doubleValue (privateFactorValue f d)
+
+_factorProduct :: (FactorElement a, GV.Vector v a, Factor (PrivateCPT v a)) => [PrivateCPT v a] -> PrivateCPT v a
+_factorProduct [] = factorFromScalar 1.0
+_factorProduct l = 
+    let nakedVars = L.foldr1 union . map factorVariables $ l
+        allVars = DVSet nakedVars
+        (scalars,cpts) = partition isScalarFactor l
+        stridesFromCPT (Table d _ _) = indexStridesFor (DVSet d) allVars
+        elementProduct [] = elementUnit
+        elementProduct l = foldr1 multiply l
+        ps = elementProduct . map (flip privateFactorValue undefined) $ scalars
+        cptsStrides = map stridesFromCPT cpts
+    in 
+    if L.null cpts 
+        then 
+            Scalar ps
+        else
+            let getFactorValueAtIndex i (strides,factor@(Table _ _ vals)) = vals!(indexPosition strides i)
+                instantiationProduct instantiation = elementProduct . map (getFactorValueAtIndex instantiation) $ (zip cptsStrides cpts)
+                values = GV.force $ GV.fromList [multiply ps (instantiationProduct x) | x <- indicesForDomain allVars]
+            in 
+            fromJust $ createCPTWithDims nakedVars values
+
+
+-- | Genereic projection function which can use any function to project out
+-- some variables in a CPT. Used to implement the MAP queries.
+cptFactorProjectOutWith :: (FactorElement a, GV.Vector v a)
+                        => ([(a,DVISet)] -> a) -- ^ Summing or maximizing function
+                        -> [DV] -- ^ Variables to sum or maximize
+                        -> PrivateCPT v a 
+                        -> PrivateCPT v a
+cptFactorProjectOutWith sumF s f@(Table d _ v) = 
+            let variablesToSum = s
+                variablesToKeep = d `difference` s 
+                keepSet = DVSet variablesToKeep
+                sumSet = DVSet variablesToSum 
+                
+                values = 
+                    let strides = indexStridesFor (DVSet d) (DVSet $ variablesToKeep ++ variablesToSum)
+                    in
+                    GV.force . GV.fromList $ do 
+                      keepIndex <- indicesForDomain keepSet 
+                      let l = do
+                            sumIndex <- indicesForDomain sumSet 
+                            let result = v!(indexPosition strides $ combineIndex strides keepIndex sumIndex) 
+                                sumInstantiation = instantiation sumSet sumIndex
+                            return (result,sumInstantiation)
+                      return $ (sumF l)
+            in
+            fromJust $ createCPTWithDims variablesToKeep values
+
+-- | Create a CPT given some dimensions and a list of Doubles.
+-- Returns nothing is the length are not coherents.
+createCPTWithDims :: GV.Vector v a => [DV] -> v a -> Maybe (PrivateCPT v a)
+createCPTWithDims [] v = Just (Scalar (GV.head v))
+createCPTWithDims dims values = 
+    let createDVIndex i (DV (Vertex v) _)  = (v,i)
+        m = IM.fromList . zipWith createDVIndex ([0,1..]::[Int]) $ dims
+        p = product (map dimension dims) 
+    in
+    if GV.length values == p
+        then
+            Just $! m `seq` Table dims m (GV.force values)
+        else 
+            Nothing
+
+
+
+{-
+
+Stride management
+
+-}
+
+
+-- | Used to combined the keep and sum indices in the factor projection
+combineIndex :: Strides s'' -> MultiIndex s -> MultiIndex s' -> MultiIndex s''
+{-# INLINE combineIndex #-}
+combineIndex _ (MultiIndex la) (MultiIndex lb) = MultiIndex (la V.++ lb)
+
+
+
+newtype Strides s = Strides (V.Vector 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 (V.toList originStrides) 
+        getNewStrides dv = maybe 0 id (lookup dv reference)
+    in 
+    Strides $ (V.force . V.fromList $ 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)
+    in 
+    Strides (V.force . V.fromList $ pos')
+
+-- | Convertion of a multiindex to its
+-- position inside of the data vector of a 'CPT'
+indexPosition :: Strides s -> MultiIndex s -> Int
+{-# INLINE indexPosition #-}
+indexPositions _ []  = 0
+indexPosition (Strides pos') (MultiIndex pos) = V.sum . V.zipWith (\x y -> x * y) pos' $ pos
+
+
diff --git a/Bayes/FactorElimination.hs b/Bayes/FactorElimination.hs
--- a/Bayes/FactorElimination.hs
+++ b/Bayes/FactorElimination.hs
@@ -1,8 +1,8 @@
-{-# LANGUAGE TypeSynonymInstances #-}
-{-# LANGUAGE FlexibleInstances #-}
 {- | Algorithms for factor elimination
 
 -}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
 module Bayes.FactorElimination(
     -- * Moral graph
       moralGraph
@@ -54,6 +54,8 @@
 import Test.QuickCheck hiding ((.||.), collect)
 import Test.QuickCheck.Arbitrary
 
+import Bayes.Factor.CPT -- This import is only used for quickcheck tests
+
 --import Debug.Trace
 --debug s a = trace (s ++ "\n" ++ show a ++ "\n") a
 
@@ -104,7 +106,7 @@
                                -> Vertex 
                                -> Vertex 
                                -> Ordering 
-nodeComparisonForTriangulation g = (compare `on` (numberOfAddedEdges g)) .||. (compare `on` (weight g))
+nodeComparisonForTriangulation g = (compare `on` (weightedEdges g))
 
 {-
 
@@ -225,7 +227,7 @@
 possibilities g currentT remaining leavesClusters = do 
   rv <- remaining
   lv <- leavesClusters
-  let NodeValue lvVertex lvCluster _ = nodeValue currentT lv
+  let NodeValue lvVertex _ _ = nodeValue currentT lv
   guard (isLinkedWithAnEdge g rv lvVertex)
   let ev = fromJust $ edgeValue g (edge rv lvVertex)
   return $ (rv,lv,ev)
@@ -305,8 +307,10 @@
 
 -- | 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 = 
+posterior :: (BayesianDiscreteVariable dv, Factor f) => JunctionTree f -> dv -> Maybe f
+posterior t someDv = 
+  let v = dv someDv
+  in
   case snd $ traverseTree (findClusterFor v) Nothing t of 
     Nothing -> Nothing
     Just c -> let NodeValue ver f e = nodeValue t c 
diff --git a/Bayes/FactorElimination/JTree.hs b/Bayes/FactorElimination/JTree.hs
--- a/Bayes/FactorElimination/JTree.hs
+++ b/Bayes/FactorElimination/JTree.hs
@@ -102,6 +102,16 @@
                         ,  separatorClusterMap :: !(Map.Map Sep c)
                         } deriving(Eq)
 
+instance FactorContainer (JTree Cluster) where 
+  changeFactor f t = 
+    let changeNodeValue (NodeValue v fa ev) = NodeValue v (changeFactor f fa) ev
+    in
+    distribute .
+    collect $
+    t { nodeValueMap = Map.map changeNodeValue (nodeValueMap t)
+      , separatorValueMap = Map.map (const EmptySeparator) (separatorValueMap t)
+      }
+
 -- | Create a singleton tree with just one root node
 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
@@ -116,38 +126,47 @@
 
 -- | Get the cluster for a separator
 separatorCluster :: JTree c a -> Sep -> c 
-separatorCluster t s = fromJust $ Map.lookup s (separatorClusterMap t)
+{-# INLINE separatorCluster #-}
+separatorCluster t s = fromJust $! Map.lookup s (separatorClusterMap t)
 
 -- | Leaves of the tree
 leaves :: JTree c a -> [c]
+{-# INLINE leaves #-}
 leaves = Set.toList . leavesSet
 
 -- | All nodes of the tree
 treeNodes :: JTree c a -> [c]
+{-# INLINE treeNodes #-}
 treeNodes = Map.keys . nodeValueMap
 
 treeValues :: JTree c f -> [(c,NodeValue f)]
+{-# INLINE treeValues #-}
 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)
+{-# INLINE nodeValue #-}
+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
+{-# INLINE setNodeValue #-}
 setNodeValue c v t = t {nodeValueMap = Map.insert c v (nodeValueMap t)} 
 
 -- | Parent of a node
 nodeParent :: Ord c => JTree c a -> c -> Maybe Sep 
-nodeParent t e = Map.lookup e (parentMap t)
+{-# INLINE nodeParent #-}
+nodeParent t e = let r = Map.lookup e (parentMap t) in r `seq` r
 
 -- | Value of a node
 separatorValue :: Ord c => JTree c a -> Sep -> SeparatorValue a 
-separatorValue t e = fromJust $ Map.lookup e (separatorValueMap t)
+{-# INLINE separatorValue #-}
+separatorValue t e = fromJust $! Map.lookup e (separatorValueMap t)
 
 -- | Parent of a separator
 separatorParent :: Ord c => JTree c a -> Sep -> c 
-separatorParent t e = fromJust $ Map.lookup e (separatorParentMap t)
+{-# INLINE separatorParent #-}
+separatorParent t e = fromJust $! Map.lookup e (separatorParentMap t)
 
 -- | UpMessage for a separator node
 upMessage :: Ord c => JTree c a -> Sep -> a
@@ -164,14 +183,17 @@
 
 -- | Return the separator childrens of a node
 nodeChildren :: Ord c => JTree c a -> c -> [Sep]
-nodeChildren t e = maybe [] id $ Map.lookup e (childrenMap t)
+{-# INLINE nodeChildren #-}
+nodeChildren t e = maybe [] id $! Map.lookup e (childrenMap t)
 
 -- | Return the child of a separator
 separatorChild :: Ord c => JTree c a -> Sep -> c 
-separatorChild t e = fromJust $ Map.lookup e (separatorChildMap t)
+{-# INLINE separatorChild #-}
+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 
+{-# INLINE nodeIsMemberOfTree #-}
 nodeIsMemberOfTree c t = Map.member c (nodeValueMap t)
 
 -- | Add a separator between two nodes.
@@ -298,6 +320,7 @@
     updateDownMessage child generatedMessage t
 
 unique :: Ord c => [c] -> [c]
+{-# INLINE unique #-}
 unique = Set.toList . Set.fromList
 
 -- | Collect message taking into account that the tree depth may be different for different leaves.
@@ -353,7 +376,7 @@
 -- for a cluster
 class IsCluster c where 
   -- | Evidence contained in the cluster
-  overlappingEvidence :: c -> [DVI Int] -> [DVI Int]
+  overlappingEvidence :: c -> [DVI] -> [DVI]
   -- | Cluser variables
   clusterVariables :: c -> [DV]
   -- | Intersection of two clusters
@@ -427,7 +450,7 @@
 
 -- | Change evidence in the network
 changeEvidence :: (IsCluster c, Ord c, Factor f, Message f c, Show c, Show f)
-               => [DVI Int] -- ^ Evidence
+               => [DVI] -- ^ Evidence
                -> JTree c f 
                -> JTree c f 
 changeEvidence e t =  
@@ -457,10 +480,10 @@
 fromCluster (Cluster s) = Set.toList s 
 
 
-instance Factor f => Message f Cluster where 
-  newMessage input (NodeValue _ f e) dv = 
+instance (Factor f) => Message f Cluster where 
+  newMessage input (NodeValue _ f e) c = 
     let allFactors = f ++ e ++ input 
-        variablesToKeep = fromCluster dv 
+        variablesToKeep = fromCluster c 
         variablesToRemove = (nub (concatMap factorVariables allFactors)) \\ variablesToKeep
     in 
     marginal allFactors variablesToRemove variablesToKeep []
diff --git a/Bayes/ImportExport/HuginNet.hs b/Bayes/ImportExport/HuginNet.hs
--- a/Bayes/ImportExport/HuginNet.hs
+++ b/Bayes/ImportExport/HuginNet.hs
@@ -13,6 +13,7 @@
 import Bayes.Factor
 import Bayes
 import Bayes.PrivateTypes
+import Bayes.Factor.CPT(changeVariableOrder)
 
 --import Debug.Trace 
 
diff --git a/Bayes/PrivateTypes.hs b/Bayes/PrivateTypes.hs
--- a/Bayes/PrivateTypes.hs
+++ b/Bayes/PrivateTypes.hs
@@ -3,36 +3,51 @@
 Those type are not exported
 
 -}
-
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
 module Bayes.PrivateTypes( 
- -- * Classes
+ -- * Classes and types
    BayesianDiscreteVariable(..)
  , Set(..)
  -- * Variables
  , DV(..)
  , DVSet(..)
+ , DVISet(..)
+ , TDV
+ , tdv
+ , tdvi
  -- * Instantiations
+ , Instantiable(..)
  , DVI(..)
  , setDVValue
- , (=:)
  , instantiationValue
  , instantiationVariable
+ , fromDVSet
  -- * Vertices 
  , Vertex(..)
  -- * Misc
  , getMinBound
  -- * Indices 
- , Index(..)
+ , MultiIndex(..)
  , forAllInstantiations 
  , indicesForDomain
- , fromIndex
  , instantiationDetails
+ , instantiation
  , allInstantiationsForOneVariable
+ -- * Tests 
+ , instantiationProp
  ) where
 
 
 import qualified Data.List as L 
+import qualified Data.Vector.Unboxed as V
+import Test.QuickCheck
+import Test.QuickCheck.Arbitrary
+import System.Random(Random)
+import qualified Data.IntMap as IM
 
+
+
 {-
 	Set
 -}
@@ -89,6 +104,8 @@
 -- | A discrete variable has a number of levels which is required to size the factors
 class BayesianDiscreteVariable v where
     dimension :: v -> Int 
+    dv :: v -> DV
+    vertex :: v -> Vertex
 
 -- | Get the minimum bound for a type
 getMinBound :: Bounded a => a -> a 
@@ -107,7 +124,7 @@
 
 -- | 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)
+newtype DVSet s = DVSet [DV] deriving(Eq,Show)
 
 -- | Remove the type tag when not needed
 fromDVSet :: DVSet s -> [DV]
@@ -118,67 +135,169 @@
 
 instance BayesianDiscreteVariable DV where
     dimension (DV _ d) = d
+    dv = id
+    vertex (DV v _) = v
 
+-- | A typed discrete variable
+data TDV s = TDV !Vertex !Int deriving(Eq,Ord)
+
+instance Show (TDV s) where
+    show (TDV v d) = show v
+
+instance BayesianDiscreteVariable (TDV s) where
+    dimension (TDV _ d) = d
+    dv (TDV v nb) = DV v nb
+    vertex (TDV v _) = v
+
+-- | Typed discrete variable
+tdv :: DV -> TDV s 
+tdv (DV v nb) = TDV v nb
+
+-- | Typed instantiation 
+tdvi :: Enum s => DVI -> (TDV s,s)
+tdvi (DVI dv value) = (tdv dv, toEnum value)
 {-
 
 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 
+newtype MultiIndex s = MultiIndex (V.Vector Int) deriving(Eq,Show)
 
+-- | Get the instantiations for a given multindex
+instantiation :: DVSet s -> MultiIndex s -> [DVI]
+instantiation (DVSet l) (MultiIndex v) = zipWith setDVValue l (V.toList v)
+
 -- | Generate all the indices for a set of variables
-indicesForDomain :: DVSet s -> [[Index s]]
-indicesForDomain (DVSet l) = mapM indicesForOneDomain l
+indicesForDomain :: DVSet s -> [MultiIndex s]
+{-# INLINE indicesForDomain #-}
+indicesForDomain (DVSet l) = map (MultiIndex . V.fromList) $ (mapM indicesForOneDomain l)
  where 
- 	indicesForOneDomain (DV _ d) = map Index [0..d-1]
+ 	indicesForOneDomain (DV _ d) = [0..d-1]
 
-allInstantiationsForOneVariable :: DV -> [DVI Int]
+allInstantiationsForOneVariable :: DV -> [DVI]
 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 s -> [[DVI]]
 forAllInstantiations (DVSet l) = mapM allInstantiationsForOneVariable l
  
 
+
+
 {- 
 
 Instantiations
 
 -}
 -- | Discrete Variable instantiation. A variable and its value
-data DVI a = DVI DV !a deriving(Eq)
+data DVI = DVI DV !Int deriving(Eq)
 
-instance Show a => Show (DVI a) where 
+instance Show (DVI) where 
    show (DVI (DV v _) i) = show v ++ "=" ++ show i
 
    -- | A set of variable instantiations
-type DVISet a = [DVI a]
+type DVISet = [DVI]
 
+class Instantiable d v where 
+  -- | Create a variable instantiation using values from
+  -- an enumeration
+  (=:) :: d -> v -> DVI 
+
+
+instance (Bounded b, Enum b) => Instantiable DV b where
+  (=:) a b = setDVValue a (fromEnum b - fromEnum (getMinBound b))
+
+instance (Bounded b, Enum b) => Instantiable (TDV b) b where
+  (=:) (TDV v nb) b = setDVValue (DV v nb) (fromEnum b - fromEnum (getMinBound b))
+
+
 -- | Create a discrete variable instantiation for a given discrete variable
-setDVValue :: DV -> a -> DVI a
+setDVValue :: DV -> Int -> DVI
 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
+
+instance BayesianDiscreteVariable DVI where
     dimension (DVI v _) = dimension v
+    dv = instantiationVariable
+    vertex (DVI dv _) = vertex dv
 
 -- | 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)
+instantiationDetails :: [DVI] -> (DVSet s, MultiIndex s)
+instantiationDetails l = (DVSet $ map instantiationVariable l, MultiIndex . V.fromList . map (instantiationValue) $ l)
 
 -- | Extract value of the instantiation
 instantiationValue (DVI _ v) = v
 
 -- | Discrete variable from the instantiation
 instantiationVariable (DVI dv _) = dv
+
+
+{-
+
+QuickCheck
+
+-}
+
+{-
+
+CPT can't have same same vertex values but with different sizes.
+But, arbitrary CPT generation will general several vertex with same vertex id
+and different vertex size.
+
+So, we introduce a function mapping a vertex ID to a vertex size. So, vertex size are hard coded
+
+-}
+
+quickCheckVertexSize :: Int -> Int
+quickCheckVertexSize 0 = 2
+quickCheckVertexSize 1 = 2
+quickCheckVertexSize 2 = 2
+quickCheckVertexSize _ = 2
+
+-- | Generate a random value until this value is not already present in the list
+whileIn :: (Arbitrary a, Eq a) => [a] -> Gen a -> Gen a
+whileIn l m = do 
+    newVal <- m 
+    if newVal `elem` l 
+        then
+            whileIn l m 
+        else 
+            return newVal
+
+-- | Generate a random vector of n elements without replacement (no duplicate)
+-- May loop if the range is too small !
+generateWithoutReplacement :: (Random a, Arbitrary a, Eq a)  
+                           => Int -- ^ Vector size
+                           -> (a,a) -- ^ Bounds
+                           -> Gen [a]
+generateWithoutReplacement n b | n == 1 = generateSingle b 
+                               | n > 1 = generateMultiple n b 
+                               | otherwise = return []
+ where
+   generateSingle b = do 
+       r <- choose b
+       return [r]
+   generateMultiple n b = do 
+       l <- generateWithoutReplacement (n-1) b
+       newelem <- whileIn l $ choose b
+       return (newelem:l)
+
+-- | Check that we can recover an instantiation from a MultiIndex
+instantiationProp :: DVSet s -> Bool 
+instantiationProp dvl = 
+    let dvs = DVSet (fromDVSet dvl) 
+    in 
+    forAllInstantiations dvs == map (instantiation dvs) (indicesForDomain dvs) 
+
+
+instance Arbitrary (DVSet s) where 
+    arbitrary =  do
+        nbVertex <- choose (1,4) :: Gen Int
+        vertexNumbers <- generateWithoutReplacement nbVertex (0,50)
+        let dimensions = map (\i -> (DV (Vertex i)  (quickCheckVertexSize i))) vertexNumbers
+        return (DVSet dimensions)
+
diff --git a/Bayes/Test.hs b/Bayes/Test.hs
--- a/Bayes/Test.hs
+++ b/Bayes/Test.hs
@@ -7,14 +7,15 @@
 import Test.Framework (defaultMain, testGroup)
 import Test.Framework.Providers.QuickCheck2 (testProperty)
 import Test.Framework.Providers.HUnit(testCase)
-import Bayes.Test.CompareEliminations(compareVariableFactor)
+import Bayes.Test.CompareEliminations(compareVariableFactor,compareFactorChange)
 
 import Bayes(testEdgeRemoval_prop,testVertexRemoval_prop)
-import Bayes.Factor(testProductProject_prop,testScale_prop,testProjectCommut_prop,testScalarProduct_prop,testProjectionToScalar_prop,testAssocProduct_prop)
+import Bayes.Factor.CPT(testProductProject_prop,testScale_prop,testProjectCommut_prop,testScalarProduct_prop,testProjectionToScalar_prop,testAssocProduct_prop)
 import Bayes.FactorElimination(junctionTreeProperty_prop,junctionTreeAllClusters_prop)
+import Bayes.PrivateTypes(instantiationProp)
 
 #ifdef LOCAL
-import Bayes.Test.ReferencePatterns(compareAsiaReference,compareCancerReference,comparePokerReference,compareFarmReference)
+import Bayes.Test.ReferencePatterns(compareAsiaReference,compareCancerReference,comparePokerReference,compareFarmReference,compareMpeCancer)
 #endif 
 
 -- | Run all the tests
@@ -31,19 +32,24 @@
                 testProperty "Commutativity of project" testProjectCommut_prop,
                 testProperty "Product with scalar factor" testScalarProduct_prop,
                 testProperty "Test projection to scalar" testProjectionToScalar_prop,
-                testProperty "Test associativity of factor" testAssocProduct_prop
+                testProperty "Test associativity of factor" testAssocProduct_prop,
+                testCase "Test factors can be changed in containers" compareFactorChange
             ]
         , testGroup "Junction Tree" [
                 testProperty "Test the junction tree property" junctionTreeProperty_prop,
                 testCase "Test variable elimination == factor elimination" compareVariableFactor,
                 testProperty "Test all clusters are included in the junction tree" junctionTreeAllClusters_prop
             ]
+        , testGroup "Misc functions" [
+                testProperty "Instantiation from multiindex" instantiationProp
+            ]
 #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
+                testCase "Farm reference pattern" compareFarmReference,
+                testCase "Test MPE and MAP with Cancer network" compareMpeCancer
         ]
 #endif
 
diff --git a/Bayes/Test/CompareEliminations.hs b/Bayes/Test/CompareEliminations.hs
--- a/Bayes/Test/CompareEliminations.hs
+++ b/Bayes/Test/CompareEliminations.hs
@@ -4,16 +4,19 @@
 
 -}
 module Bayes.Test.CompareEliminations(
-    compareVariableFactor
+   compareVariableFactor
+ , compareFactorChange
  ) where
 
 import Test.HUnit.Lang(assertFailure)
 
-import Bayes.Examples(example)
+import Bayes.Examples(example,exampleWithFactorChange)
 import Bayes.Factor
 import Bayes
 import Bayes.VariableElimination
 import Bayes.FactorElimination
+import Bayes.Factor.CPT
+import Data.Maybe(fromJust)
 
 compareFactors :: String -> Maybe CPT -> CPT -> IO ()
 compareFactors s Nothing _ = assertFailure s
@@ -24,24 +27,59 @@
         else 
             assertFailure s
 
+compareFactorChange :: IO () 
+compareFactorChange = do 
+    let ([winter,sprinkler,rain,wet,road,roadandrain],exampleG) = (map dv . fst $ example,snd example)
+        (_,exampleGF) = exampleWithFactorChange 
+        theNewFactor = fromJust $ factorWithVariables [wet,sprinkler,rain] [1,0.2,0.1,0.05,0,0.8,0.9,0.95]
+        exampleGF' = changeFactor theNewFactor exampleGF
+        compareG s a b = compareFactors s (Just $ priorMarginal exampleGF' a b) (priorMarginal exampleG a b)
+    compareG "PRIOR FOR RAIN" [winter,sprinkler,wet,road,roadandrain] [rain]
+    compareG "PRIOR FOR WINTER" [sprinkler,wet,road,rain,roadandrain] [winter]
+    compareG "PRIOR FOR SPRINKLER" [winter,wet,road,rain,roadandrain] [sprinkler]
+    compareG "PRIOR FOR WET" [winter,sprinkler,road,rain,roadandrain] [wet]
+    compareG "PRIOR FOR ROAD" [winter,sprinkler,wet,rain,roadandrain] [road]
+
+    -- Now we test that the junction tree factors can be changed and the result will be right.
+    -- We first create a wrong junction tree using the wrong graph (where factor values are wrong)
+    let jt' = createJunctionTree nodeComparisonForTriangulation exampleGF 
+    -- We then update the values to the right ones by using the right factor
+        jt = changeFactor theNewFactor jt'
+
+    compareFactors "PRIOR FOR RAIN" (posterior jt rain) (priorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain])
+
+    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,roadandrain] [rain]  [wet =: True])
+    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt2 rain) 
+         (posteriorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain]  [wet =: True, sprinkler =: True])
+
+    compareFactors "PRIOR FOR WINTER" (posterior jt winter) (priorMarginal exampleG [sprinkler,wet,road,rain,roadandrain] [winter])
+    compareFactors "PRIOR FOR SPRINKLER" (posterior jt sprinkler) (priorMarginal exampleG [winter,wet,road,rain,roadandrain] [sprinkler])
+    compareFactors "PRIOR FOR WET" (posterior jt wet) (priorMarginal exampleG [winter,sprinkler,road,rain,roadandrain] [wet])
+    compareFactors "PRIOR FOR ROAD" (posterior jt road) (priorMarginal exampleG [winter,sprinkler,wet,rain,roadandrain] [road])
+
+
 -- | Compare that variable elemination and factor elimination are giving
 -- similar results on a simple example
 compareVariableFactor :: IO ()
 compareVariableFactor = do 
-    let ([winter,sprinkler,rain,wet,road],exampleG) = example
+    let ([winter,sprinkler,rain,wet,road,roadandrain],exampleG) = example
         jt = createJunctionTree nodeComparisonForTriangulation exampleG
-    compareFactors "PRIOR FOR RAIN" (posterior jt rain) (priorMarginal exampleG [winter,sprinkler,wet,road] [rain])
+    compareFactors "PRIOR FOR RAIN" (posterior jt rain) (priorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain])
 
     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])
+         (posteriorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain]  [wet =: True])
     compareFactors "POSTERIOR RAIN FOR WET" (posterior jt2 rain) 
-         (posteriorMarginal exampleG [winter,sprinkler,wet,road] [rain]  [wet =: True, sprinkler =: True])
+         (posteriorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain]  [wet =: True, sprinkler =: True])
 
-    compareFactors "PRIOR FOR WINTER" (posterior jt winter) (priorMarginal exampleG [sprinkler,wet,road,rain] [winter])
-    compareFactors "PRIOR FOR SPRINKLER" (posterior jt sprinkler) (priorMarginal exampleG [winter,wet,road,rain] [sprinkler])
-    compareFactors "PRIOR FOR WET" (posterior jt wet) (priorMarginal exampleG [winter,sprinkler,road,rain] [wet])
-    compareFactors "PRIOR FOR ROAD" (posterior jt road) (priorMarginal exampleG [winter,sprinkler,wet,rain] [road])
+    compareFactors "PRIOR FOR WINTER" (posterior jt winter) (priorMarginal exampleG [sprinkler,wet,road,rain,roadandrain] [winter])
+    compareFactors "PRIOR FOR SPRINKLER" (posterior jt sprinkler) (priorMarginal exampleG [winter,wet,road,rain,roadandrain] [sprinkler])
+    compareFactors "PRIOR FOR WET" (posterior jt wet) (priorMarginal exampleG [winter,sprinkler,road,rain,roadandrain] [wet])
+    compareFactors "PRIOR FOR ROAD" (posterior jt road) (priorMarginal exampleG [winter,sprinkler,wet,rain,roadandrain] [road])
 
diff --git a/Bayes/Test/ReferencePatterns.hs b/Bayes/Test/ReferencePatterns.hs
--- a/Bayes/Test/ReferencePatterns.hs
+++ b/Bayes/Test/ReferencePatterns.hs
@@ -10,6 +10,7 @@
  , compareCancerReference
  , comparePokerReference
  , compareFarmReference
+ , compareMpeCancer
 #endif
  ) where
 
@@ -19,6 +20,7 @@
 import Bayes.Factor
 import Bayes
 import Bayes.FactorElimination
+import Bayes.VariableElimination(mpe)
 import Bayes.Examples(anyExample)
 import Bayes.FactorElimination.JTree(root)
 
@@ -42,6 +44,24 @@
 (~=~) a b = and (zipWith comparePercent a b)
 
 #ifdef LOCAL
+
+data Positive = Yes | No deriving(Eq,Enum,Bounded,Show) 
+
+--  | Get the name of vertices in the graph for display
+rename :: NamedGraph g => g a b -> (TDV s,s) -> (String,s)
+rename g = \(a,s) -> (fromJust . vertexLabel g . vertex $ a, s)
+
+-- | Test that a MAP is not always the projection of a MPE
+compareMpeCancer = do 
+  (varmap,g) <- anyExample "asia.net"
+  let [x,b,d,a,s,l,t,e] = map tdv . fromJust $ mapM (flip Map.lookup varmap) ["X","B","D","A","S","L","T","E"] :: [TDV Positive]
+      m = mpe g [x,d] [b,a,s,l,t,e] [x =: Yes, d =: No]
+      typedResult = map (map (rename g . tdvi)) m :: [[(String,Positive)]]
+  assertBool "Test MPE" $ typedResult == [[("E",No),("T",No),("L",No),("S",No),("B",No),("A",No)]]
+  let m = mpe g [x,d,b,l,t,e] [a,s] [x =: Yes, d =: No]
+      typedResult = map (map (rename g . tdvi)) m :: [[(String,Positive)]]
+  assertBool "Test MAP" $ typedResult == [[("S",Yes),("A",No)]]
+
 compareFarmReference = do 
   (varmap,g) <- anyExample "studfarm.net"
   let jt = createJunctionTree nodeComparisonForTriangulation g
diff --git a/Bayes/Tools.hs b/Bayes/Tools.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Tools.hs
@@ -0,0 +1,22 @@
+{- | Miscellaneous tools
+	
+-}
+module Bayes.Tools (
+	nearlyEqual
+	) where 
+
+
+
+-- | Floating point number comparisons which should take into account
+-- all the subtleties of that kind of comparison
+nearlyEqual :: Double -> Double -> Bool
+nearlyEqual a b = 
+    let absA = abs a 
+        absB = abs b 
+        diff = abs (a-b)
+        epsilon = 2e-5
+    in
+    case (a,b) of 
+        (x,y) | x == y -> True -- handle infinities
+              | x*y == 0 -> diff < (epsilon * epsilon)
+              | otherwise -> diff / (absA + absB) < epsilon
diff --git a/Bayes/VariableElimination.hs b/Bayes/VariableElimination.hs
--- a/Bayes/VariableElimination.hs
+++ b/Bayes/VariableElimination.hs
@@ -12,6 +12,8 @@
  , minFillOrder
  , allVariables
  , marginal
+ , mpemarginal
+ , mpe
  , EliminationOrder
  ) where
 
@@ -21,13 +23,17 @@
 import Data.Maybe(fromJust)
 import Data.Function(on)
 import qualified Data.Map as M
+import Bayes.Factor.PrivateCPT(convertToMaxFactor,CPT,MAXCPT)
+import Bayes.Factor.CPT 
+import Bayes.Factor.MaxCPT
+import Bayes.PrivateTypes(DVISet)
 
 --import Debug.Trace 
 
 --debug s a = trace (s  ++ "\n" ++ show a ++ "\n") a
 
 -- | Elimination order
-type EliminationOrder = [DV]
+type EliminationOrder dv = [dv]
 
 -- | Get all variables from a Bayesian Network
 allVariables :: (Graph g, Factor f) 
@@ -40,12 +46,22 @@
   map createDV s
 
 -- | Used for bucket elimination. Factor are organized by their first DV
-data Buckets f = Buckets !EliminationOrder !(M.Map DV [f])
+data Buckets f = Buckets !(EliminationOrder DV) !(M.Map DV [f])
 
+instance Show f => Show (Buckets f) where 
+  show (Buckets v m) = "BUCKET\n" ++ show v ++ "\n" ++ concatMap disp (M.toList m)
+   where
+    disp (v,f) = "Bucket for " ++ show v ++ "\n" ++ concatMap dispElem f ++ "\n----\n"
+    dispElem f = show f ++ "\n"
+
+convertToMaxCPT :: Buckets CPT -> Buckets MAXCPT 
+convertToMaxCPT (Buckets e m) = Buckets e (M.map (map convertToMaxFactor) m) 
+
+
 createBuckets ::  (Factor f) 
               => [f] -- ^ Factor to use for computing the marginal one
-              -> EliminationOrder -- ^ Variables to eliminate
-              -> EliminationOrder -- ^ Remaining variables
+              -> EliminationOrder DV-- ^ Variables to eliminate
+              -> EliminationOrder DV -- ^ Remaining variables
               -> Buckets f 
 createBuckets s e r = 
   let -- We put the selected variables for elimination in the right order at the beginning
@@ -58,7 +74,7 @@
         (remaining, M.insert dv fk m)
       (_,b) = foldl' addDVToBucket (s,M.empty) theOrder
   in
-  Buckets (tail theOrder) b
+  Buckets theOrder b
 
 -- | Get the factors for a bucket
 getBucket :: DV 
@@ -75,12 +91,14 @@
 updateBucket dv f b@(Buckets e m) = 
   if isScalarFactor f 
     then 
-      Buckets (tail e) (M.insert dv [f] m)
+      Buckets (remainingVarsToProcess e) (M.insert dv [f] m)
     else
       let b' = removeFromBucket dv b
-          Buckets e' m' = addBucket b' f 
-      in 
-      Buckets (tail e') m'
+      in
+      addBucket b' f 
+ where 
+  remainingVarsToProcess [] = []
+  remainingVarsToProcess l = tail l
 
 -- | Add a factor to the right bucket
 addBucket :: Factor f => Buckets f -> f -> Buckets f
@@ -93,15 +111,24 @@
 
 -- | Remove a variable from the bucket
 removeFromBucket :: DV -> Buckets f -> Buckets f 
-removeFromBucket dv (Buckets e m) = Buckets e (M.delete dv m) 
+removeFromBucket dv (Buckets [] m) = Buckets [] (M.delete dv m) 
+removeFromBucket dv (Buckets e m) = Buckets (tail e) (M.delete dv m) 
 
+marginalizeOneVariable :: Factor f => Buckets f -> DV -> Buckets f
+marginalizeOneVariable currentBucket dv   = 
+  let fk = getBucket dv currentBucket
+      p = factorProduct fk
+      f' = factorProjectOut [dv] p
+  in
+  updateBucket dv f' currentBucket
+
 -- | Compute the prior marginal. All the variables in the
 -- elimination order are conditionning variables ( p( . | conditionning variables) )
-marginal :: (Factor f) 
+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
+         -> EliminationOrder DV -- ^ Ordering of variables to marginalize
+         -> EliminationOrder DV -- ^ Ordering of remaining variables
+         -> [DVI] -- ^ Assignment for some factors in variables to marginalize
          -> f
 marginal lf p r assignment = 
   -- The elimintation order are the variables to eliminate.
@@ -115,22 +142,58 @@
   in
   -- We get P(Q , e)
   resultFactor
- where 
-  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) 
+-- | Compute the prior marginal. All the variables in the
+-- elimination order are conditionning variables ( p( . | conditionning variables) )
+-- First we sum, then we maximize for the remaining variables
+mpemarginal :: [CPT] -- ^ Bayesian Network
+            -> EliminationOrder DV -- ^ Ordering of variables to marginalize
+            -> EliminationOrder DV -- ^ Ordering of remaining variables
+            -> [DVI] -- ^ Assignment for some factors in variables to marginalize
+            -> MAXCPT
+mpemarginal lf p r assignment = 
+  -- The elimintation order are the variables to eliminate.
+  -- But the algorithm also needs the remaining variables
+  let bucket = createBuckets lf p r
+      assignmentFactors = map factorFromInstantiation assignment
+      bucket' = foldl' addBucket bucket assignmentFactors
+      bucket'' = foldl' marginalizeOneVariable bucket' p
+      bucketMax = convertToMaxCPT bucket''
+      Buckets _ resultBucket = foldl' marginalizeOneVariable bucketMax r
+      resultFactor = factorProduct  . concat . M.elems $ resultBucket
+      -- The norm is P(e) and result factor is P(Q,e)
+  in
+  -- We get P(Q , e)
+  resultFactor
+ 
+
+-- | Most Probable Explanation (or Maximum A Posteriori estimator)
+-- when restricted to a subest of variables in output
+mpe :: (Graph g, BayesianDiscreteVariable dva, BayesianDiscreteVariable dvb) 
+    => BayesianNetwork g CPT -- ^ Bayesian network defining the factors
+    -> EliminationOrder dva -- ^ Ordering of variables to sum out (should contain evidence variables)
+    -> EliminationOrder dvb -- ^ Ordering of remaining variables (to maximize)
+    -> [DVI] -- ^ Assignment
+    -> [DVISet] -- ^ MPE or MAP instantiation
+mpe g someP someR assignment = 
+    let p = map dv someP
+        r = map dv someR
+        s = allVertexValues g 
+        resultFactor = mpemarginal s p r assignment
+    in 
+    mpeInstantiations (resultFactor)
+
+posteriorMarginal :: (Graph g, Factor f, Show f, BayesianDiscreteVariable dva, BayesianDiscreteVariable dvb) 
                   => 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
+                  -> EliminationOrder dva -- ^ Ordering of variables to marginzalie
+                  -> EliminationOrder dvb-- ^ Ordering of remaining variables
+                  -> [DVI] -- ^ Assignment for some factors in variables to marginalize
                   -> f
-posteriorMarginal g p r assignment = 
-  let s = allVertexValues g 
+posteriorMarginal g someP someR assignment = 
+  let p = map dv someP 
+      r = map dv someR
+      s = allVertexValues g 
       resultFactor = marginal s p r assignment
       norm = factorNorm resultFactor
   in
@@ -139,13 +202,15 @@
 
 -- | Compute the prior marginal. All the variables in the
 -- elimination order are conditionning variables ( p( . | conditionning variables) )
-priorMarginal :: (Graph g, Factor f, Show f) 
+priorMarginal :: (Graph g, Factor f, Show f, BayesianDiscreteVariable dva, BayesianDiscreteVariable dvb) 
               => BayesianNetwork g f -- ^ Bayesian Network
-              -> EliminationOrder -- ^ Ordering of variables to marginalize
-              -> EliminationOrder -- ^ Ordering of remaining to keep in result
+              -> EliminationOrder dva-- ^ Ordering of variables to marginalize
+              -> EliminationOrder dvb-- ^ Ordering of remaining to keep in result
               -> f
-priorMarginal g ea eb = 
-  let s = allVertexValues g 
+priorMarginal g someEA someEB = 
+  let ea = map dv someEA 
+      eb = map dv someEB
+      s = allVertexValues g 
       resultFactor = marginal s ea eb []
       norm = factorNorm resultFactor
   in
@@ -191,7 +256,7 @@
 -- | Compute the degree order of an elimination order
 degreeOrder :: (FoldableWithVertex g, Factor f, Graph g)
             => BayesianNetwork g f
-            -> EliminationOrder 
+            -> EliminationOrder DV 
             -> Int 
 degreeOrder g p =
   let  ig = interactionGraph g :: UndirectedSG () DV
@@ -216,7 +281,7 @@
 eliminationOrderForMetric :: (Graph g, Factor f, FoldableWithVertex g, UndirectedGraph g')
                           => (g' () DV -> DV -> Int)
                           -> BayesianNetwork g f 
-                          -> EliminationOrder 
+                          -> EliminationOrder DV  
 eliminationOrderForMetric metric g = 
   let ig = interactionGraph g
       s = allVertexValues ig
@@ -232,11 +297,11 @@
 -- | Elimination order minimizing the degree
 minDegreeOrder :: (Graph g, Factor f, FoldableWithVertex g)
                => BayesianNetwork g f 
-               -> EliminationOrder 
+               -> EliminationOrder DV  
 minDegreeOrder = eliminationOrderForMetric nbNeighbors
 
 -- | Elimination order minimizing the filling
 minFillOrder :: (Graph g, Factor f, FoldableWithVertex g)
                => BayesianNetwork g f 
-               -> EliminationOrder 
+               -> EliminationOrder DV  
 minFillOrder = eliminationOrderForMetric nbMissingLinks
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.1
+Version:             0.3
 
 -- A short (one-line) description of the package.
 Synopsis:            Inference with Discrete Bayesian Networks
@@ -15,7 +15,8 @@
 -- 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.1 version is using new faster and cleaner algorithms and correcting a problem with graph triangulation.
+ examples where it worked. It should be considered as experimental and not used
+ in any production work.
 
 -- URL for the project homepage or repository.
 Homepage:            http://www.alpheccar.org
@@ -34,8 +35,10 @@
 Maintainer:          misc@alpheccar.org
 
 -- A copyright notice.
-Copyright: Copyright (c) 2012, alpheccar       
+Copyright: Copyright (c) 2012, alpheccar 
 
+Stability: experimental      
+
 Category:            Math
 
 Build-type:          Simple
@@ -61,6 +64,8 @@
   Exposed-modules:
     Bayes
     Bayes.Factor
+    Bayes.Factor.CPT
+    Bayes.Factor.MaxCPT
     Bayes.ImportExport.HuginNet
     Bayes.VariableElimination
     Bayes.FactorElimination
@@ -74,8 +79,10 @@
     Bayes.ImportExport.HuginNet.Splitting
     Bayes.PrivateTypes
     Bayes.FactorElimination.JTree
+    Bayes.Tools
+    Bayes.Factor.PrivateCPT
 
-  GHC-Options: -O2 -funbox-strict-fields
+  GHC-Options: -funbox-strict-fields
   Extensions: CPP
   if flag(local)
     cpp-options: -DLOCAL