diff --git a/Bayes.hs b/Bayes.hs
--- a/Bayes.hs
+++ b/Bayes.hs
@@ -1,18 +1,18 @@
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# 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
-examples where it worked. On bigger networks, imported from Hugin files, it was very very very slow.
-So, you can use this software as a toy. Much more work is needed to validate
-and optimize it. 
+examples where it worked. It should be considered as experimental and not used
+in any production work.
 
 Look at the "Bayes.Examples" and "Bayes.Examples.Tutorial" in this package 
 to see how to use the library.
 
+In "Bayes.Examples.Influence" you'll find additional examples about influence diagrams.
+
 -}
 module Bayes(
   -- * Graph
@@ -22,7 +22,6 @@
   , DirectedGraph(..)
   , FoldableWithVertex(..)
   , NamedGraph(..)
-  , Distribution(..)
   -- ** Graph Monad
   , GraphMonad
   , GMState(..)
@@ -30,11 +29,18 @@
   , runGraph
   , execGraph
   , evalGraph
+  , emptyAuxiliaryState
+  , getNewEmptyVariable
+  , isRoot
+  , rootNode
+  , parentNodes
+  , childrenNodes
   -- ** Support functions for Graph constructions
   , Vertex
   , Edge 
   , edge
   , newEdge
+  , getVertex
   , edgeEndPoints
   , connectedGraph
   , dag
@@ -43,34 +49,11 @@
   -- ** The SimpleGraph type
   , DirectedSG
   , UndirectedSG
+  , SBN(..)
+  , varMap
+  , displaySimpleGraph
   -- ** Bayesian network
-  , SBN
   , BayesianNetwork(..)
-  -- * Bayesian Monad used to ease creation of Bayesian Networks
-  , BNMonad
-  , runBN 
-  , evalBN
-  , execBN
-  -- ** Variable creation
-  , variable
-  , unamedVariable
-  , variableWithSize
-  , tdv
-  , t
-  -- ** Creation of conditional probability tables
-  , cpt
-  , proba
-  , (~~)
-  , softEvidence
-  , se
-  -- ** Creation of truth tables
-  , logical 
-  , (.==.)
-  , (.!.)
-  , (.|.)
-  , (.&.)
-  -- ** Noisy OR
-  , noisyOR
   -- * Testing
   , testEdgeRemoval_prop
   , testVertexRemoval_prop
@@ -82,23 +65,27 @@
 import Control.Monad.Writer.Strict
 import Control.Applicative((<$>))
 import Bayes.Factor hiding(isEmpty)
+import Bayes.Factor.CPT(CPT(..))
+import Bayes.Factor.MaxCPT(MAXCPT(..))
 import Data.Maybe
 import qualified Data.Map as Map
 import qualified Data.Foldable as F
 import qualified Data.Traversable as T 
 import Control.Applicative 
 import qualified Data.Set as Set
-import qualified Data.List as L(find)
 
 import Test.QuickCheck hiding ((.&.),Testable)
 import Test.QuickCheck.Arbitrary
-import Data.List(sort,intercalate,nub)
+import Data.List(sort,intercalate,nub,foldl')
 import Bayes.PrivateTypes hiding(isEmpty)
 import GHC.Float(float2Double)
 
 --import Debug.Trace
 --debug a = trace (show a) a
 
+-- | An implementation of the BayesianNetwork using the simple graph and no value for the edges
+type SBN f = DirectedSG () f
+
 -- | Bayesian network. g must be a directed graph and f a factor
 type BayesianNetwork g f = g () f
 
@@ -292,6 +279,27 @@
     ingoing :: g a b -> Vertex -> Maybe [Edge]
     outgoing :: g a b -> Vertex -> Maybe [Edge]
 
+-- | Return the parents of a node
+parentNodes :: DirectedGraph g => g a b -> Vertex -> [Vertex]
+parentNodes g v = maybe [] id $ do 
+  ie <- ingoing g v
+  p <- mapM (startVertex g) ie
+  return p
+
+-- | Return the children of a node
+childrenNodes :: DirectedGraph g => g a b -> Vertex -> [Vertex]
+childrenNodes g v = maybe [] id $ do 
+  ie <- outgoing g v
+  p <- mapM (endVertex g) ie
+  return p
+
+isRoot :: DirectedGraph g => g a b -> Vertex -> Bool
+{-# INLINE isRoot #-}
+isRoot g v =
+  case ingoing g v of 
+    Just [] -> True 
+    _ -> False
+
 -- | Get the root node for the graph
 rootNode :: DirectedGraph g => g a b -> Maybe Vertex
 rootNode g = 
@@ -300,11 +308,6 @@
   case someRoots of 
     (h:l) -> Just h 
     _ -> Nothing
-  where 
-    isRoot g v =
-      case ingoing g v of 
-        Just [] -> True 
-        _ -> False
 
 -- | Check if the graph is a directed Acyclic graph
 dag :: DirectedGraph g => g a b -> Bool 
@@ -333,12 +336,6 @@
 
 
 
-                          
-
-
--- | Edge type used to identify and edge in a graph
-data Edge = Edge !Vertex !Vertex deriving(Eq,Ord,Show)
-
 -- | Create an edge description
 edge :: Vertex -> Vertex -> Edge 
 {-# INLINE edge #-}
@@ -349,15 +346,8 @@
 edgeEndPoints (Edge va vb) = (va,vb)
 
 
--- | Synonym for undefined because it is clearer to use t to set the Enum bounds of a variable
-t = undefined
 
--- | Neighborhood structure for directed or undirected edges
--- | Directed edges
-data DE = DE ![Edge] ![Edge] deriving(Eq,Show)
 
--- | Undirected edges
-data UE = UE ![Edge] deriving(Eq,Show)
 
 -- | Class used to share as much code as possible between
 -- directed and undirected graphs without
@@ -402,24 +392,18 @@
   addOutgoingEdge e (UE l) = UE (e:l)
   addIngoingEdge e (UE l) = UE (e:l)
 
--- | Implementtaion of a SimpleGraph
-data SimpleGraph local edgedata vertexdata = SP {
- -- | Mapping of edge to edge data
-    edgeMap :: !(M.Map Edge edgedata) 
- -- ^ Mapping of vertex number to vertex neighborhood and vertex data
- ,  vertexMap :: !(IM.IntMap (local, vertexdata))
- -- ^ Vertex names. Used only to generate the graphviz representtaion. Names are useless for the algorithms
- -- and I don't want them to appear in the vetex values which should only be factor. Otherwise, the algorithms
- -- are less elegant since I have to extract the factors from the values
- , nameMap :: !(IM.IntMap String)
- } 
 
+
 -- | Directed simple graph
 type DirectedSG = SimpleGraph DE
 
 -- | Undirected simple graph
 type UndirectedSG = SimpleGraph UE
 
+-- | Get the variable name mapping
+varMap :: SimpleGraph n e v -> M.Map String Vertex 
+varMap (SP _ _ n) = M.fromList . map (\(i,s) -> (s, Vertex i)) . IM.toList $ n
+
 instance (Eq a, Eq b) => Eq (SimpleGraph DE a b) where
     (==) (SP a b _) (SP a' b' _) = a == a' && b == b'
 
@@ -668,6 +652,40 @@
 Following code is used to display a graph in a form adapted to humans.
 
 -}
+
+bracketS :: String -> String 
+bracketS [] = []
+bracketS s = " [" ++ s ++ "];"
+
+createNodeStyle :: (MonadWriter String m) 
+                => (Vertex -> n -> Maybe String)
+                -> (Vertex -> n -> Maybe String)
+                -> Maybe String 
+                -> Vertex 
+                -> n 
+                -> m ()
+createNodeStyle nodeShape nodeColor maybeLabel v n = 
+  let apply f = f v n
+      label _ _ = case maybeLabel of 
+         Nothing -> Nothing 
+         Just s -> Just $ "label=\"" ++ s ++ "\""
+  in 
+  tell $ bracketS . intercalate "," . mapMaybe apply $ [nodeShape,nodeColor, label]
+  
+
+createEdgeStyle :: (MonadWriter String m) 
+                => (Edge -> e -> Maybe String)
+                -> (Edge -> e -> Maybe String)
+                -> Edge
+                -> e 
+                -> m ()
+createEdgeStyle edgeShape edgeColor e n = 
+  let apply f = f e n
+  in 
+  tell $ bracketS . intercalate "," . mapMaybe apply $ [edgeShape,edgeColor]
+
+
+
 printNode nm (Vertex k,v) = do 
   tell "\n"
   let r = IM.lookup k nm
@@ -676,59 +694,66 @@
   tell "\n"
   tell $ show v
   tell "\n"
-addVertexToGraphviz nm (k,(_,v)) = do
+
+addVertexToGraphviz nodeShape nodeColor nm (k,(_,v)) = do
   tell $ show k
   let r = IM.lookup k $ nm 
-  when (isJust r) $ do
-    tell " [label=\""
-    tell $ fromJust r
-    tell "\"] ;" 
+  createNodeStyle nodeShape nodeColor r (Vertex k) v
   tell "\n"
 
+addVertexToUndirectedGraphviz nm (k,(_,v)) = do
+  tell $ show k
+  tell "\n"
+
 -- | Print the values of the graph vertices
 printGraphValues :: (Graph (SimpleGraph n), Show b) => SimpleGraph n e b -> IO () 
 printGraphValues g@(SP _ _ nm) = putStrLn . execWriter $ mapM_ (printNode nm) (allNodes g)
 
-instance (Show b, Show e) => Show (DirectedSG e b)where
-  show g@(SP em vm nm) = execWriter $ do
+displaySimpleGraph :: (Vertex -> n -> Maybe String)
+                   -> (Vertex -> n -> Maybe String)
+                   -> (Edge -> e -> Maybe String)
+                   -> (Edge -> e -> Maybe String)
+                   -> SimpleGraph local e n 
+                   -> String 
+displaySimpleGraph  nodeShape nodeColor edgeShape edgeColor g@(SP em vm nm) = execWriter $ do 
   tell "digraph dot {\n"
-  mapM_ (addVertexToGraphviz nm) $ IM.toList vm
+  mapM_ (addVertexToGraphviz nodeShape nodeColor  nm) $ IM.toList vm
   tell "\n"
-  mapM_ addEdgeToGraphviz $ M.toList em
+  mapM_ (addEdgeToGraphviz edgeShape edgeColor ) $ M.toList em
   tell "}\n"
    where
-     addEdgeToGraphviz (Edge (Vertex vs) (Vertex ve),l) = do
+     addEdgeToGraphviz es ec (e@(Edge (Vertex vs) (Vertex ve)),l) = do
        tell $ show vs 
        tell " -> "
        tell $ show ve
-       tell " [label=\""
-       tell $ show l
-       tell "\"]"
-       tell ";\n"
+       createEdgeStyle es ec e l
+       tell "\n"  
 
+noNodeStyle _ _ = Nothing 
+noEdgeStyle _ _ = Nothing
+
+instance Show (DirectedSG () CPT) where
+  show g = displaySimpleGraph noNodeStyle noNodeStyle noEdgeStyle noEdgeStyle g
+
+instance Show (DirectedSG () MAXCPT) where
+  show g = displaySimpleGraph noNodeStyle noNodeStyle noEdgeStyle noEdgeStyle g
+
+instance Show (DirectedSG String String) where
+  show g = displaySimpleGraph noNodeStyle noNodeStyle noEdgeStyle noEdgeStyle g
+
 instance (Show b, Show e) => Show (UndirectedSG e b)where
   show g@(SP em vm nm) = execWriter $ do
   tell "graph dot {\n"
-  mapM_ (addVertexToGraphviz nm) $ IM.toList vm
+  mapM_ (addVertexToUndirectedGraphviz nm) $ IM.toList vm
   tell "\n"
-  mapM_ addEdgeToGraphviz $ M.toList em
+  mapM_ (addEdgeToGraphviz) $ M.toList em
   tell "}\n"
    where
-     addEdgeToGraphviz (Edge (Vertex vs) (Vertex ve),l) = do
+     addEdgeToGraphviz (e@(Edge (Vertex vs) (Vertex ve)),l) = do
        tell $ show vs 
        tell " -- "
        tell $ show ve
-       tell " [label=\""
-       tell $ show l
-       tell "\"]"
-       tell ";\n"
-
-
--- | Bayesian variable : name,dimension, factor
--- When initialized it is using a factor with bayesian variables.
--- But the factor value are not yet set
-data MaybeBNode f = UninitializedBNode String Int
-                  | InitializedBNode String Int f
+       tell "\n"
 
 
 displayFactors :: (NeighborhoodStructure n, Show f, Factor f, Graph (SimpleGraph n)) => SimpleGraph n a f -> String
@@ -741,21 +766,20 @@
   in
   intercalate "\n" $ map displayFactor nodes
 
--- | An implementation of the BayesianNetwork using the simple graph and no value of edges
-type SBN f = DirectedSG () f
+{-
 
+Graph Monad
+
+-}
+
+
+
 -- | State used for the construction of the graph in the monad and containing
 -- auxiliary informations like vertex name to vertex id and vertex count
 type AuxiliaryState = (M.Map String Int, Int)
 
 emptyAuxiliaryState = (M.empty,0)
 
--- | The State for the monad with a mapping from variable name to variable ID.
-type BNState g f = GMState g () (MaybeBNode f)
-
--- | The Bayesian monad
-type BNMonad g f a = GraphMonad g () (MaybeBNode f) a
-
 -- | The state of the graph monad : the graph and auxiliary data
 -- useful during the construction
 type GMState g e f = (AuxiliaryState,g e f)
@@ -765,18 +789,6 @@
 -- g is the graph type. e the edge type. f the node type (generally a 'Factor')
 newtype GraphMonad g e f a = GM {runGraphMonad :: State (GMState g e f) a} deriving(Monad, MonadState (GMState g e f))
 
--- | Get the Bayesian Discrete Variable for a vertex.
--- It works because we keep the variable dimension
-factorVariable :: Graph g => Vertex -> BNMonad g f (Maybe DV)  
-factorVariable v = do 
-  g <- gets snd 
-  let value = vertexValue g v
-  case value of
-    Nothing -> return Nothing
-    Just (UninitializedBNode _ d) -> return $ Just $ DV v d
-    Just (InitializedBNode _ d _) -> return $ Just $ DV v d
-  
-
 -- | Get a named vertex from the graph monad
 getVertex :: Graph g => String -> GraphMonad g e f (Maybe Vertex)
 getVertex a = do
@@ -785,10 +797,6 @@
     i <- M.lookup a namemap
     return (Vertex i)
 
--- | Create an edge between two vertex of the Bayesian network
-(<--) :: Graph g => DV -> DV -> BNMonad g f ()
-DV va _ <-- DV vb _ = newEdge vb va ()
-
 -- | Add a new labeled edge to the graph
 newEdge :: Graph g => Vertex -> Vertex -> e -> GraphMonad g e f ()
 newEdge va vb e = do
@@ -797,97 +805,6 @@
   put $! (aux,g1)
   return ()
 
-whenJust Nothing _ = return ()
-whenJust (Just i) f = f i >> return ()
-
--- | Get the node of a bayesian network under creation
-getBayesianNode :: Graph g => Vertex -> BNMonad g f (Maybe (MaybeBNode f))
-getBayesianNode v = do
-  g <- gets snd
-  return $ vertexValue g v
-
--- | Set the node of a bayesian network under creation
-setBayesianNode :: Graph g => Vertex -> MaybeBNode f -> BNMonad g f ()
-setBayesianNode v newValue = do
-  (aux,oldGraph) <- get
-  let newGraph = changeVertexValue v newValue oldGraph
- 
-  whenJust newGraph $ \nvm -> do
-     put $! (aux, nvm)
-
--- | 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, 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 >>= return . dv -- This is updating the state and so the graph
-  g <- gets snd
-  current <- factorVariable v
-  mvalue <- getBayesianNode v
-  maybe (return ()) (setCpt g l v current) mvalue
-
-
-    
-minBoundForEnum :: Bounded a => a -> a
-minBoundForEnum _ = minBound
-
-maxBoundForEnum :: Bounded a => a -> a
-maxBoundForEnum _ = maxBound
-
-intValue :: Enum a => a -> Int
-intValue = fromEnum
-
-
--- | Set the bound of a bayesian variable (number of levels)
-setVariableBoundWithSize :: Graph g
-                         => Vertex -- ^ Vertex
-                         -> Int -- ^ Inf limit  (0 for instance)
-                         -> Int -- ^ Sup limit (1 for instance for 2 elements)
-                         -> BNMonad g f ()
-setVariableBoundWithSize a bmin bmax = do
-    v <- getBayesianNode a
-    whenJust v $ \(UninitializedBNode s _) -> do
-      setBayesianNode a (UninitializedBNode s (bmax - bmin + 1))
-
-setVariableBound :: (Enum a, Bounded a, Graph g) 
-                 => Vertex -- ^ Vertex
-                 -> a -- ^ Bounded variable (t :: type where t is undefined)
-                 -> BNMonad g f ()
-setVariableBound a e = 
-  let bmin = intValue $ minBoundForEnum e
-      bmax = intValue $ maxBoundForEnum e
-  in 
-  setVariableBoundWithSize a bmin bmax
-
--- | Create a new named Bayesian variable if not found.
--- Otherwise, return the found one.
-addVariableIfNotFound :: NamedGraph g => String -> BNMonad g f Vertex
-addVariableIfNotFound vertexName = graphNode vertexName (UninitializedBNode vertexName 0)
-
 -- | Add a node in the graph using the graph monad
 graphNode :: NamedGraph g => String -> f -> GraphMonad g e f Vertex 
 graphNode vertexName initValue = do
@@ -906,206 +823,6 @@
   put $! ((namemap1,count+1),g1)
   return (Vertex count)
 
--- | Initialize a new variable
-_initializeNewVariable :: (Enum a, Bounded a, NamedGraph g)
-                       => Vertex 
-                       -> a 
-                       -> 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 (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 (TDV a)
-unamedVariable e = do 
-  va <- getNewEmptyVariable Nothing (UninitializedBNode "unamed" 0)
-  _initializeNewVariable va e
-
--- | Define a Bayesian variable (name and bounds)
-variable :: (Enum a, Bounded a, NamedGraph g) 
-        => String -- ^ Variable name
-        -> a -- ^ Variable bounds
-        -> BNMonad g f (TDV a)
-variable name e = do
-  va <- addVariableIfNotFound name
-  _initializeNewVariable va e
-
--- | Define a Bayesian variable (name and bounds)
-variableWithSize :: NamedGraph g
-        => String -- ^ Variable name
-        -> Int -- ^ Variable size
-        -> BNMonad g f DV
-variableWithSize name e = do
-  va <- addVariableIfNotFound name
-  _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)
-  case fromJust maybeValue of 
-     UninitializedBNode s d -> return (DV va d)
-     InitializedBNode _ d _ -> return (DV va d)
-
--- | Define a conditional probability between different variables
--- Variables are ordered like
--- FFF FFT FTF FTT TFF TFT TTF TTT
--- and same for other enumeration keeping enumeration order
-cpt :: (DirectedGraph g , BayesianDiscreteVariable v,BayesianDiscreteVariable vb) => v -> [vb] -> BNMonad g f v
-cpt node conditions = do
-  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, 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)
 runGraph = removeAuxiliaryState . flip runState (emptyAuxiliaryState,emptyGraph) . runGraphMonad 
  where 
@@ -1116,29 +833,3 @@
 
 execGraph :: Graph g => GraphMonad g e f a -> g e f
 execGraph = snd . flip execState (emptyAuxiliaryState,emptyGraph) . runGraphMonad 
-
--- | Create a bayesian network using the simple graph implementation
--- The initialized nodes are replaced by the factor.
--- Returns the monad values and the built graph.
-runBN :: BNMonad DirectedSG f a -> (a,DirectedSG () f)
-runBN x = 
-  let (r,g) = runGraph x
-      convertBNodes (InitializedBNode s d f) = f 
-      convertBNodes (UninitializedBNode s d) = error $ "All variables must be initialized with a factor: " ++ s ++ "(" ++ show d ++ ")"
-  in 
-  (r,fmap convertBNodes g)
-
--- | Create a bayesian network but only returns the monad value.
--- Mainly used for testing.
-evalBN :: BNMonad DirectedSG f a -> a
-evalBN = evalGraph
-
--- | Create a bayesian network but only returns the monad value.
--- Mainly used for testing.
-execBN :: BNMonad DirectedSG f a -> DirectedSG () f
-execBN x = 
-  let g = execGraph x
-      convertBNodes (InitializedBNode s d f) = f 
-      convertBNodes (UninitializedBNode s d) = error $ "All variables must be initialized with a factor: " ++ s ++ "(" ++ show d ++ ")"
-  in 
-  fmap convertBNodes g
diff --git a/Bayes/BayesianNetwork.hs b/Bayes/BayesianNetwork.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/BayesianNetwork.hs
@@ -0,0 +1,224 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{- | Module for building Bayesian Networks
+
+-}
+module Bayes.BayesianNetwork(
+  -- * Bayesian Monad used to ease creation of Bayesian Networks
+    BNMonad
+  , runBN 
+  , evalBN
+  , execBN
+  , Distribution(..)
+  -- ** Variable creation
+  , variable
+  , unamedVariable
+  , variableWithSize
+  , tdv
+  , t
+  -- ** Creation of conditional probability tables
+  , cpt
+  , proba
+  , (~~)
+  , softEvidence
+  , se
+  -- ** Creation of truth tables
+  , logical 
+  , (.==.)
+  , (.!.)
+  , (.|.)
+  , (.&.)
+  -- ** Noisy OR
+  , noisyOR
+  ) where
+
+import Bayes
+import Bayes.PrivateTypes
+import Control.Monad.State.Strict
+import Bayes.Factor
+import Data.Maybe(fromJust)
+import qualified Data.List as L(find)
+import Data.List(sort,intercalate,nub)
+import Bayes.Tools(minBoundForEnum,maxBoundForEnum,intValue)
+import Bayes.Network 
+
+-- | Synonym for undefined because it is clearer to use t to set the Enum bounds of a variable
+t = undefined
+
+-- | The Bayesian monad
+type BNMonad g f a = NetworkMonad g () f a
+
+-- | Initialize the values of a factor
+(~~) :: (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 >>= return . dv -- This is updating the state and so the graph
+  maybeNewValue <- getCpt v l
+  currentValue <- getBayesianNode v
+  case (currentValue, maybeNewValue) of 
+    (Just c, Just n) -> initializeNodeWithValue v c n
+    _ -> return ()
+
+
+-- | Define a conditional probability between different variables
+-- Variables are ordered like
+-- FFF FFT FTF FTT TFF TFT TTF TTT
+-- and same for other enumeration keeping enumeration order
+-- Note that the reverse is important. We add the parents in such a way that 'ingoing'
+-- will give a list of parents in the right order.
+-- This order must correspond to the order of values in the initialization.
+cpt :: (DirectedGraph g , BayesianDiscreteVariable v,BayesianDiscreteVariable vb) => v -> [vb] -> BNMonad g f v
+cpt node conditions = do
+  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, 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.  
+
+-}
+
+-- | Create a  network using the simple graph implementation
+-- The initialized nodes are replaced by the value.
+-- Returns the monad values and the built graph.
+runBN :: BNMonad DirectedSG f a -> (a,SBN f)
+runBN = runNetwork
+
+-- | Create a  network but only returns the monad value.
+-- Mainly used for testing.
+execBN :: BNMonad DirectedSG f a -> SBN f
+execBN = execNetwork
+
+-- | Create a bayesian network but only returns the monad value.
+-- Mainly used for testing.
+evalBN :: BNMonad DirectedSG f a -> a
+evalBN = evalNetwork
+
+
diff --git a/Bayes/Examples.hs b/Bayes/Examples.hs
--- a/Bayes/Examples.hs
+++ b/Bayes/Examples.hs
@@ -152,7 +152,7 @@
 import System.Directory(getHomeDirectory)
 import System.FilePath((</>))
 import Bayes.Factor.CPT
-
+import Bayes.BayesianNetwork
 
 #ifndef LOCAL
 import Paths_hbayes
diff --git a/Bayes/Examples/Influence.hs b/Bayes/Examples/Influence.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Examples/Influence.hs
@@ -0,0 +1,289 @@
+{-# LANGUAGE ViewPatterns #-}
+{- | Examples of influence diagrams
+
+An influence diagram is an extension of a Bayesian network with can be used to solve some decision problems.
+In an influence diagram, there are two new kind of nodes : decision nodes and utility nodes.
+
+Solving an influence diagram means determining the strategies for each decision variable that will maximize the average utility.
+
+There must be an ordering of the decision variables : a path through all the decisions.
+
+A decision variable can depend on other past decisions and probabilistic nodes. In the later case, the variable of 
+the probabilistic node is assumed to be observed before the decision is taken. So, the decision is only trying to 
+maximize the average utility based on what has not been observed (the future and some past probabilistic variables).
+
+A probabilistic node can depend on other probabilistic nodes (like in a Bayesian network) and decision nodes.
+
+An utility is a leaf of the graph.
+
+/Example graph/
+
+Building an influence diagram is done like for a Bayesian network : by using the right monad.
+
+@
+import Bayes.InfluenceDiagram 
+studentSimple = snd . 'runID' $ do
+@
+
+Then, you create the different nodes of the graph:
+
+@
+    e <- 'decisionNode' \"E\" ('t' :: E)
+    uc <- 'utilityNode' \"UC\"
+    ub <- 'utilityNode' \"UB\"
+    i <- 'chance' "I" ('t' :: I)
+    pr <- 'chance' "P" ('t' :: Bool)
+@
+
+The types used above are:
+
+@
+data E = Dont | Do deriving(Eq,Enum,Bounded)
+data I = Low | Average | High deriving(Eq,Enum,Bounded)
+@
+
+Then, you need to define the dependencies and the numerical values. For probabilistic nodes, it is done like
+for Bayesian network:
+
+@
+    cpt pr ['d' e] ~~ [1-0.0000001,1 - 0.001,0.0000001, 0.001]
+    cpt i ['p' pr, 'd' e] ~~ [0.2,0.1,0.01,0.01,0.6,0.5,0.04,0.04,0.2,0.4,0.95,0.95]
+@
+
+The list may contain decision variables of type 'DEV' and probabilistic variables of type 'DV' or 'TDV'. So, the 
+functions 'p' an 'd' are used for the embedding in the heterogenous list.
+
+For decision nodes, the method is similar but with two differences : The first decision may depend on nothing (just on the assumed future).
+And there are no values to define for a decision variable since the goal of the influence diagram is to compute them.
+
+@
+    'decision' e 'noDependencies'
+@
+
+For the utility nodes, it is similar to probabilistic nodes. You define the dependencies and the numerical values:
+
+@
+    'utility' uc [e] ~~ [0,-50000]
+    'utility' ub [i] ~~ [100000,200000,500000]
+@
+  
+Once the influence diagram is defined, you can solve it:
+
+@
+    'solveInfluenceDiagram' studentSimple
+@
+
+The result of this function is the solution : the decision strategies. You may want to display also the original
+graph to see to which node are corresponding the vertex numbers.
+
+/Policy Network/
+
+You can transform a solved influence diagram into a policy network : a Bayesian network where decision variables have been replaced
+with probabilistic variables where the conditional probability table is containing 1 for a choice of variables corresponding
+to the decision and 0 otherwise.
+
+@
+    let l = 'solveInfluenceDiagram' student
+        g = 'policyNetwork' l student
+    print g 
+    'printGraphValues' g
+@ 
+
+-}
+module Bayes.Examples.Influence(
+    -- * Influence diagrams
+      exampleID
+    , student
+    , studentSimple
+    , market 
+    -- * Variables for some networks
+    , studentDecisionVars
+    , studentSimpleDecisionVar
+    -- * Tests for the networks
+    , theTest
+    , policyTest
+    , marketTest
+    ) where 
+
+import Bayes.InfluenceDiagram 
+import Bayes(printGraphValues)
+import Bayes.Factor(forAllInstantiations,dv,instantiationValue,DVSet(..))
+
+-- | Very simple example with one decision node
+exampleID :: InfluenceDiagram
+exampleID = snd . runID $ do 
+    a <- chance "A" (t :: Bool)
+    d1 <- decisionNode "D" (t :: Bool)
+    u <- utilityNode "U"
+
+    proba a ~~ [0.8,0.1]
+    decision d1 [a]
+    utility u [d d1,p a] ~~ [1,10,8,2]
+    return ()
+
+data E = Dont | Do deriving(Eq,Enum,Bounded)
+data I = Low | Average | High deriving(Eq,Enum,Bounded)
+data S = Found | DontFound deriving(Eq,Enum,Bounded)
+
+studentSimpleDecisionVar :: DEV 
+
+-- | Student network as found in the book by Barber
+studentSimple :: InfluenceDiagram
+(studentSimpleDecisionVar,studentSimple) = runID $ do 
+    e <- decisionNode "E" (t :: E)
+
+    uc <- utilityNode "UC"
+    ub <- utilityNode "UB"
+
+    i <- chance "I" (t :: I)
+    pr <- chance "P" (t :: Bool)
+
+    cpt pr [d e] ~~ [1-0.0000001,1 - 0.001,0.0000001, 0.001]
+
+    cpt i [p pr, d e] ~~ [0.2,0.1,0.01,0.01,0.6,0.5,0.04,0.04,0.2,0.4,0.95,0.95]
+    decision e noDependencies
+
+    utility uc [e] ~~ [0,-50000]
+    utility ub [i] ~~ [100000,200000,500000]
+    return e
+
+-- | Solve the influences diagrams for the both student network.
+-- Also displays each network
+theTest = do
+    print studentSimple
+    printGraphValues studentSimple
+    putStrLn "RESULT"
+    print $ solveInfluenceDiagram studentSimple
+    putStrLn "----"
+    print student
+    printGraphValues student
+    putStrLn "RESULT"
+    print $ solveInfluenceDiagram student
+
+
+-- | Solve the influence diagram 'student' and convert it into
+-- a policy network
+policyTest = do 
+    print student 
+    printGraphValues student
+    let l = solveInfluenceDiagram student
+        g = policyNetwork l student
+    print g 
+    printGraphValues g
+
+studentDecisionVars :: (DEV,TDV Bool,DEV)   
+
+-- | Student network as found in the book by Barber
+student :: InfluenceDiagram
+(studentDecisionVars,student) = runID $ do 
+    e <- decisionNode "E" (t :: E)
+    s <- decisionNode "S" (t :: S)
+
+    uc <- utilityNode "UC"
+    ub <- utilityNode "UB"
+    us <- utilityNode "US"
+
+    pr <- chance "P" (t :: Bool)
+    i <- chance "I" (t :: I)
+
+    cpt pr [d e] ~~ [1-0.0000001,1 - 0.001,0.0000001, 0.001]
+
+    cpt i [p pr, d s] ~~ [0.2,0.1,0.05,0.005, 0.6,0.5,0.15,0.005,0.2,0.4,0.8,0.99]
+    decision s [pr]
+    decision e noDependencies
+
+    utility uc [e] ~~ [0,-50000]
+    utility ub [i] ~~ [100000,200000,500000]
+    utility us [s] ~~ [0,-200000]
+    return (e,pr,s)
+
+
+{- 
+
+Test with a market network
+    
+-}
+data F = Forecast | NoForecast deriving(Eq,Enum,Bounded)
+data IN = Choice0 | Choice1 | Choice2 deriving(Eq,Enum,Bounded)
+data EF = Up | Flat | Down deriving(Eq,Enum,Bounded)
+
+genValues :: ([DVI] -> Double) -> [DV] -> [Double]
+genValues f l = [f x | x <- forAllInstantiations (DVSet l)]
+
+e :: Enum a => DVI -> a
+e = toEnum . instantiationValue 
+
+vf :: [DVI] -> (EF,F,EF)
+vf [a,b,c] = (e a, e b, e c)
+vf _ = (toEnum 0, toEnum 0, toEnum 0)
+
+uf :: [DVI] -> (EF,IN,F)
+uf [a,b,c] = (e a, e b, e c)
+uf _ = (toEnum 0, toEnum 0, toEnum 0)
+
+getForecastUtility (uf -> (Up, Choice0, Forecast)) = 1500
+getForecastUtility (uf -> (Up, Choice0, NoForecast)) = 1500
+getForecastUtility (uf -> (Up, Choice1, Forecast)) = 1000
+getForecastUtility (uf -> (Up, Choice1, NoForecast)) = 1000
+getForecastUtility (uf -> (Up, Choice2, Forecast)) = 500
+getForecastUtility (uf -> (Up, Choice2, NoForecast)) = 500
+
+getForecastUtility (uf -> (Flat, Choice0, Forecast)) = 100
+getForecastUtility (uf -> (Flat, Choice0, NoForecast)) = 100
+getForecastUtility (uf -> (Flat, Choice1, Forecast)) = 200
+getForecastUtility (uf -> (Flat, Choice1, NoForecast)) = 200
+getForecastUtility (uf -> (Flat, Choice2, Forecast)) = 500
+getForecastUtility (uf -> (Flat, Choice2, NoForecast)) = 500
+
+getForecastUtility (uf -> (Down, Choice0, Forecast)) = -1000
+getForecastUtility (uf -> (Down, Choice0, NoForecast)) = -1000
+getForecastUtility (uf -> (Down, Choice1, Forecast)) = -100
+getForecastUtility (uf -> (Down, Choice1, NoForecast)) = -100
+getForecastUtility (uf -> (Down, Choice2, Forecast)) = 500
+getForecastUtility (uf -> (Down, Choice2, NoForecast)) = 500
+
+getForecastProba (vf -> (Up,Forecast,Up)) = 0.8
+getForecastProba (vf -> (Up,Forecast,Flat)) = 0.15
+getForecastProba (vf -> (Up,Forecast,Down)) = 0.2
+getForecastProba (vf -> (Up,NoForecast,Up)) = 0.33
+getForecastProba (vf -> (Up,NoForecast,Flat)) = 0.33
+getForecastProba (vf -> (Up,NoForecast,Down)) = 0.33
+
+getForecastProba (vf -> (Flat,Forecast,Up)) = 0.1
+getForecastProba (vf -> (Flat,Forecast,Flat)) = 0.7
+getForecastProba (vf -> (Flat,Forecast,Down)) = 0.2
+getForecastProba (vf -> (Flat,NoForecast,Up)) = 0.33
+getForecastProba (vf -> (Flat,NoForecast,Flat)) = 0.33
+getForecastProba (vf -> (Flat,NoForecast,Down)) = 0.33
+
+getForecastProba (vf -> (Down,Forecast,Up)) = 0.1
+getForecastProba (vf -> (Down,Forecast,Flat)) = 0.15
+getForecastProba (vf -> (Down,Forecast,Down)) = 0.6
+getForecastProba (vf -> (Down,NoForecast,Up)) = 0.33
+getForecastProba (vf -> (Down,NoForecast,Flat)) = 0.33
+getForecastProba (vf -> (Down,NoForecast,Down)) = 0.33
+
+-- | Market diagram
+market :: InfluenceDiagram
+market = snd . runID $ do 
+    o <- decisionNode "Obtain Forecast" (t :: F)
+    i <- decisionNode "Investment" (t :: IN)
+
+    ef <- chance "Economy Forecast" (t :: EF)
+    ma <- chance "Market Activity" (t :: EF)
+
+    u <- utilityNode "Payoff"
+
+    proba ma ~~ [0.5,0.3,0.2]
+    decision o noDependencies 
+    decision i [d o,p ef]
+    cpt ef [d o, p ma] ~~ (genValues getForecastProba [dv ef, dv o, dv ma])
+    utility u [p ma, d i, d o] ~~ (genValues getForecastUtility [dv ma, dv i, dv o])
+    return ()
+
+-- | Solve the 'market' influence diagram
+marketTest = do 
+    print market 
+    printGraphValues market
+    let l = solveInfluenceDiagram market
+    print l 
diff --git a/Bayes/Examples/Tutorial.hs b/Bayes/Examples/Tutorial.hs
--- a/Bayes/Examples/Tutorial.hs
+++ b/Bayes/Examples/Tutorial.hs
@@ -249,6 +249,8 @@
 import Data.Maybe(fromJust,mapMaybe)
 import System.Exit(exitSuccess)
 import qualified Data.List as L((\\))
+import Bayes.BayesianNetwork(se)
+
 
 #ifdef LOCAL
 miscDiabete = do 
diff --git a/Bayes/Factor.hs b/Bayes/Factor.hs
--- a/Bayes/Factor.hs
+++ b/Bayes/Factor.hs
@@ -6,6 +6,8 @@
 module Bayes.Factor(
  -- * Factor
    Factor(..)
+ , Distribution(..)
+ , MultiDimTable(..)
  , isomorphicFactor
  , normedFactor
  , displayFactorBody
@@ -19,10 +21,11 @@
  -- ** Discrete variables and instantiations
  , DV(..)
  , TDV
- --, DVSet(..)
+ , DVSet(..)
  , DVI
  , DVISet
  , tdvi
+ , tdv
  , setDVValue
  , instantiationValue
  , instantiationVariable
@@ -38,11 +41,23 @@
 import Bayes.Tools
 import qualified Data.Vector.Unboxed as V
 import Text.PrettyPrint.Boxes hiding((//))
+import Bayes.VariableElimination.Buckets(IsBucketItem(..))
 
 --import Debug.Trace
 
 --debug a = trace ("\nDEBUG\n" ++ show a ++ "\n") a
 
+
+
+-- | 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)
+
+
 -- | 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
@@ -59,7 +74,7 @@
 -- 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 
+   changeFactor :: (IsBucketItem f,Factor f) => f -> m f -> m f 
 
 instance FactorContainer [] where 
     changeFactor = changeFactorInFunctor
@@ -171,15 +186,6 @@
         in 
         factorProjectOut s' f
 
-
-
-
-
-
-
-
-
-
 -- | 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
@@ -203,16 +209,22 @@
 Following functions are used to typeset the factor when displaying it
 
 -}
+
+-- | Class used to display multidimensional tables
+class MultiDimTable f where 
+    elementStringValue :: f -> [DVI] -> String
+    tableVariables :: f -> [DV]   
+
 -- | Display a variable name and its size
 vname :: Int -> DVI -> Box
 vname vc i = text $ "v" ++ show vc ++ "=" ++ show (instantiationValue i)
 
-dispFactor :: Factor f => f -> DV -> [DVI] -> [DV] -> Box
+dispFactor :: MultiDimTable f => f -> DV -> [DVI] -> [DV] -> Box
 dispFactor cpt h c [] = 
     let dstIndexes = allInstantiationsForOneVariable h
         dependentIndexes =  reverse c
         factorValueAtPosition p = 
-            let v = factorStringValue cpt p
+            let v = elementStringValue cpt p
             in
             text v
     in
@@ -223,9 +235,9 @@
     in
     hsep 1 top . map (\i -> vcat center1 [vname vc i,dispFactor cpt dst (i:c) l])  $ allInst
 
-displayFactorBody :: Factor f => f -> String 
+displayFactorBody :: MultiDimTable f => f -> String 
 displayFactorBody c = 
-        let d = factorVariables c
+        let d = tableVariables c
             h@(DV (Vertex vc) _) = head d
             table = dispFactor c h [] (tail d)
             dstIndexes = map head (forAllInstantiations . DVSet $ [h])
diff --git a/Bayes/Factor/CPT.hs b/Bayes/Factor/CPT.hs
--- a/Bayes/Factor/CPT.hs
+++ b/Bayes/Factor/CPT.hs
@@ -10,6 +10,8 @@
    -- * CPT Factor
      CPT
    , changeVariableOrder
+   , cptDivide
+   , cptSum
    -- * Tests
    , testProductProject_prop
    , testAssocProduct_prop
@@ -29,6 +31,8 @@
 import Data.Maybe(fromJust,mapMaybe,isJust)
 import Bayes.Factor.PrivateCPT
 import Bayes.PrivateTypes
+import Bayes.VariableElimination.Buckets(IsBucketItem(..))
+import qualified Data.Vector.Unboxed as V
 
 -- | Soft evidence factor can be used to initialize a factor
 --instance Distribution CPT where 
@@ -108,6 +112,9 @@
     scale = (*)
     multiply = (*)
     elementUnit = 1.0
+    divide _ 0 = 0
+    divide a b = a / b 
+    elementSum = (+)
 
 
 instance Factor CPT where
@@ -129,13 +136,32 @@
     
     evidenceFrom = _evidenceFrom
 
-    factorProduct = _factorProduct
+    factorProduct = _factorProduct multiply
     
     factorProjectOut _ f@(Scalar v) = f
     factorProjectOut s f= cptFactorProjectOutWith (sum . map fst) s f
 
-
+--  Divie two CPT
+cptDivide :: CPT -> CPT -> CPT
+cptDivide a b = _factorProduct multiply [a,invertCPT b]
 
+invertCPT (Scalar a) = Scalar (1.0 / a)
+invertCPT (Table d m v) = Table d m (V.map inv v)
+ where 
+    inv x = 1.0 / x
+-- Warning, the fold in factorProduct may not do the operations in the order you expect.
+-- Since a/b is different from b/a we have to take this into account. So, the order is changed in the list
 
+cptSum :: [CPT] -> CPT
+cptSum = _factorProduct elementSum
 
+instance IsBucketItem CPT where 
+    scalarItem = isScalarFactor
+    itemProduct = factorProduct
+    itemProjectOut d = factorProjectOut [d]
+    itemContainsVariable = containsVariable
 
+instance MultiDimTable CPT where 
+    elementStringValue = factorStringValue
+    tableVariables = factorVariables
+      
diff --git a/Bayes/Factor/MaxCPT.hs b/Bayes/Factor/MaxCPT.hs
--- a/Bayes/Factor/MaxCPT.hs
+++ b/Bayes/Factor/MaxCPT.hs
@@ -16,6 +16,7 @@
 import qualified Data.IntMap as IM
 import Data.List(maximumBy)
 import Data.Function(on)
+import Bayes.VariableElimination.Buckets(IsBucketItem(..))
 
 --import Debug.Trace 
 
@@ -29,6 +30,8 @@
     multiply (da,la) (db,[]) = (da*db,la)
     multiply (da,la) (db,lb) = (da*db,[xa ++ xb | xa <- la, xb <- lb])
     elementUnit = (1.0,[])
+    divide _ _ = error "It does not make sense to divide MAXCPT elements"
+    elementSum _ _ = error "It does not make sense to sum MAXCPT elements"
 
 instantiationProba :: ((Double,PossibleInstantiations), DVISet) -> Double
 instantiationProba (a,b) = fst a 
@@ -62,7 +65,7 @@
 
     evidenceFrom = _evidenceFrom
 
-    factorProduct = _factorProduct
+    factorProduct = _factorProduct multiply
     
     factorProjectOut _ f@(Scalar v) = f
     factorProjectOut s f = cptFactorProjectOutWith maximization s f
@@ -73,8 +76,17 @@
 
 instance Show MAXCPT where
     show (Scalar v) = "\nScalar Factor:\n" ++ show v
-    show c@(Table [] _ v) = "\nEmpty CPT:\n"
+    show c@(Table [] _ v) = "\nEmpty MAXCPT:\n"
 
     show c = displayFactorBody c  
 
+instance IsBucketItem MAXCPT where 
+    scalarItem = isScalarFactor
+    itemProduct = factorProduct
+    itemProjectOut d = factorProjectOut [d]
+    itemContainsVariable = containsVariable
+
+instance MultiDimTable MAXCPT where 
+    elementStringValue = factorStringValue
+    tableVariables = factorVariables
 
diff --git a/Bayes/Factor/PrivateCPT.hs b/Bayes/Factor/PrivateCPT.hs
--- a/Bayes/Factor/PrivateCPT.hs
+++ b/Bayes/Factor/PrivateCPT.hs
@@ -1,6 +1,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE TypeSynonymInstances #-}
 {-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 {- | Private functions used to implement CPT
 
 -}
@@ -9,6 +10,7 @@
      cptFactorProjectOutWith
    -- ** Private functions
    , FactorElement(..)
+   , DecisionFactor(..)
    , _emptyFactor
    , _factorVariables
    , _isScalarFactor
@@ -29,6 +31,8 @@
    , PrivateCPT(..)
    , PossibleInstantiations
    , convertToMaxFactor
+   , convertToNormalFactor 
+   , decisionFactor
    , debugCPT
    , privateFactorValue
   ) where
@@ -64,7 +68,7 @@
                     , mapping :: !(IM.IntMap Int) -- ^ Mapping from vertex number to position in dimensions
                     , values :: !(v a) -- ^ Table of values
                     }
-                    | Scalar !a 
+                    | Scalar !a  deriving(Eq)
 
 debugCPT (Scalar d) = do 
    putStrLn "SCALAR CPT"
@@ -88,10 +92,33 @@
   in
   Table d m newValues
 
-type CPT = PrivateCPT (V.Vector) Double
+type DecisionFactor = PrivateCPT (NV.Vector) DVI 
 
 
 
+
+-- | Return an array of instantiations. Used to compute decisions
+-- in an influence diagram
+decisionFactor :: MAXCPT -> DecisionFactor
+decisionFactor (Scalar (_,l)) = Scalar (head . head $ l) 
+decisionFactor (Table d m v) =
+  let extractElem = head . head -- We only have ONE instantiated variable in influence diagram. So no need
+      -- for a list of list (list of several variables)
+      newValues = NV.fromList . map extractElem . map snd . NV.toList $ v
+  in
+  Table d m newValues
+
+
+
+convertToNormalFactor :: MAXCPT -> CPT 
+convertToNormalFactor (Scalar a) = Scalar (fst a)
+convertToNormalFactor (Table d m v) =
+  let newValues = V.fromList . map fst . NV.toList $ v
+  in
+  Table d m newValues
+
+type CPT = PrivateCPT (V.Vector) Double
+
 type PossibleInstantiations = [DVISet]
 
 type MAXCPT = PrivateCPT (NV.Vector) (Double,PossibleInstantiations)
@@ -101,6 +128,8 @@
     mkValue :: Double -> a
     scale :: Double -> a -> a
     multiply :: a -> a -> a
+    divide :: a -> a -> a
+    elementSum ::  a -> a -> a
     elementUnit :: a
 
 _isUsingSameVariablesAs :: PrivateCPT v a -> PrivateCPT v a -> Bool
@@ -157,7 +186,7 @@
     in 
     fromJust $ createCPTWithDims (factorVariables f) newValues
 
-privateFactorValue :: (FactorElement a, GV.Vector v a) => PrivateCPT v a  -> [DVI] -> a
+privateFactorValue :: (GV.Vector v a) => PrivateCPT v a  -> [DVI] -> a
 {-# INLINE privateFactorValue #-}
 privateFactorValue (Scalar v) _ = v 
 privateFactorValue f@(Table d _ v) i = 
@@ -181,15 +210,18 @@
 {-# 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 = 
+_factorProduct :: (FactorElement a, GV.Vector v a, Factor (PrivateCPT v a)) 
+               => (a -> a -> a)
+               -> [PrivateCPT v a] 
+               -> PrivateCPT v a
+_factorProduct _ [] = factorFromScalar 1.0
+_factorProduct op 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
+        elementProduct l = foldr1 op l
         ps = elementProduct . map (flip privateFactorValue undefined) $ scalars
         cptsStrides = map stridesFromCPT cpts
     in 
@@ -199,7 +231,7 @@
         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]
+                values = GV.force $ GV.fromList [op ps (instantiationProduct x) | x <- indicesForDomain allVars]
             in 
             fromJust $ createCPTWithDims nakedVars values
 
diff --git a/Bayes/FactorElimination.hs b/Bayes/FactorElimination.hs
--- a/Bayes/FactorElimination.hs
+++ b/Bayes/FactorElimination.hs
@@ -53,7 +53,9 @@
 
 import Test.QuickCheck hiding ((.||.), collect)
 import Test.QuickCheck.Arbitrary
+import Bayes.VariableElimination.Buckets(IsBucketItem(..))
 
+
 import Bayes.Factor.CPT -- This import is only used for quickcheck tests
 
 --import Debug.Trace
@@ -293,7 +295,7 @@
   maximumSpanningTree g''
 
 -- | Create a function tree
-createJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f, Show f)
+createJunctionTree :: (DirectedGraph g, FoldableWithVertex g, NamedGraph g, Factor f, IsBucketItem f, Show f)
                   => (UndirectedSG () f -> Vertex -> Vertex -> Ordering) -- ^ Weight function on the moral graph
                   -> BayesianNetwork g f -- ^ Input directed graph
                   -> JunctionTree f -- ^ Junction tree
@@ -307,7 +309,7 @@
 
 -- | Compute the marginal posterior (if some evidence is set on the junction tree)
   -- otherwise compute just the marginal prior.
-posterior :: (BayesianDiscreteVariable dv, Factor f) => JunctionTree f -> dv -> Maybe f
+posterior :: (BayesianDiscreteVariable dv, Factor f, IsBucketItem f) => JunctionTree f -> dv -> Maybe f
 posterior t someDv = 
   let v = dv someDv
   in
diff --git a/Bayes/FactorElimination/JTree.hs b/Bayes/FactorElimination/JTree.hs
--- a/Bayes/FactorElimination/JTree.hs
+++ b/Bayes/FactorElimination/JTree.hs
@@ -50,9 +50,11 @@
 import Bayes
 import Data.Function(on)
 import Bayes.VariableElimination(marginal)
+import Data.Binary
+import Bayes.VariableElimination.Buckets(IsBucketItem(..))
 
-import Debug.Trace 
-debug s a = trace (s ++ " " ++ show a ++ "\n") a
+--import Debug.Trace 
+--debug s a = trace (s ++ " " ++ show a ++ "\n") a
 
 type UpMessage a = a 
 type DownMessage a = Maybe a
@@ -76,7 +78,7 @@
 instance Show a => Show (NodeValue a) where 
     show (NodeValue v f e) = "f(" ++ show f ++ ") e(" ++ show e ++ ")"
 
-newtype Sep = Sep Int deriving(Eq,Ord,Show,Num)
+newtype Sep = Sep Int deriving(Eq,Ord,Show,Num,Binary)
 
 -- | Junction tree.
 -- 'c' is the node / separator identifier (for instance a set of 'DV')
@@ -480,7 +482,7 @@
 fromCluster (Cluster s) = Set.toList s 
 
 
-instance (Factor f) => Message f Cluster where 
+instance (Factor f,IsBucketItem f) => Message f Cluster where 
   newMessage input (NodeValue _ f e) c = 
     let allFactors = f ++ e ++ input 
         variablesToKeep = fromCluster c 
diff --git a/Bayes/ImportExport.hs b/Bayes/ImportExport.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/ImportExport.hs
@@ -0,0 +1,175 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{- | Import / export Bayesian networks and junction tress
+
+-}
+module Bayes.ImportExport (
+	-- * Networks
+	  writeNetworkToFile
+	, readNetworkFromFile
+	-- * Junction Tree 
+	, writeVariableMapAndJunctionTreeToFile
+	, readVariableMapAndJunctionTreeToFile
+	) where 
+
+import Data.Binary
+import Bayes
+import Bayes.PrivateTypes(Vertex(..),Edge(..), SimpleGraph(..),DE(..), UE(..),DV(..))
+import Bayes.Factor.PrivateCPT(CPT(..), MAXCPT(..),PrivateCPT(..)) 
+import System.FilePath
+import qualified Data.Vector as NV
+import qualified Data.Vector.Unboxed as V
+import Bayes.FactorElimination.JTree(SeparatorValue(..),NodeValue(..),JTree(..),Cluster(..),JunctionTree(..))
+import qualified Data.Map as M 
+
+-- | Write a bayesian network to file
+writeNetworkToFile :: FilePath -- ^ File path
+                   -> SBN CPT -- ^ Bayesian network
+                   -> IO () 
+writeNetworkToFile f n = encodeFile f n 
+
+-- | Read bayesian network from file 
+readNetworkFromFile :: FilePath 
+                    -> IO (SBN CPT)
+readNetworkFromFile = decodeFile
+
+-- | Write a junction tree and the variable map to a file 
+writeVariableMapAndJunctionTreeToFile :: FilePath
+                                      -> (M.Map String Vertex)
+                                      -> JunctionTree CPT 
+                                      -> IO ()
+writeVariableMapAndJunctionTreeToFile f vm jt = encodeFile f (vm,jt)
+
+-- | Read variable map and junction tree from file
+readVariableMapAndJunctionTreeToFile :: FilePath 
+                                     -> IO (M.Map String Vertex, JunctionTree CPT)
+readVariableMapAndJunctionTreeToFile f = decodeFile f
+
+instance Binary Cluster where 
+	put (Cluster s) = put s 
+	get = get >>= return . Cluster
+
+instance (Ord c, Binary c, Binary f) => Binary (JTree c f) where 
+	put (JTree r ls cm pm spm scm nvm svm sck sclm) = do 
+		put r 
+		put ls 
+		put cm 
+		put pm 
+		put spm 
+		put scm 
+		put nvm 
+		put svm 
+		put sck 
+		put sclm 
+	get = do 
+		r <- get
+		ls <- get 
+		cm <- get 
+		pm <- get 
+		spm <- get 
+		scm <- get 
+		nvm <- get 
+		svm <- get 
+		sck <- get 
+		sclm <- get
+		return $ JTree r ls cm pm spm scm nvm svm sck sclm
+
+instance Binary a => Binary (NodeValue a) where 
+	put (NodeValue v f e) = do 
+		put v 
+		put f 
+		put e 
+	get = do 
+		v <- get 
+		f <- get 
+		e <- get
+		return $ NodeValue v f e
+
+instance Binary a => Binary (SeparatorValue a) where
+	put (EmptySeparator) = do 
+		putWord8 0 
+	put (SeparatorValue a b) = do 
+		putWord8 1 
+		put a 
+		put b 
+	get = do 
+		tag <- getWord8 
+		case tag of 
+			0 -> return EmptySeparator
+			_ -> do 
+				a <- get 
+				b <- get 
+				return $ SeparatorValue a b 
+
+instance Binary (V.Vector Double) where 
+	put = put . V.toList 
+	get = get >>= return . V.fromList 
+
+instance Binary (NV.Vector Double) where 
+	put = put . NV.toList 
+	get = get >>= return . NV.fromList 
+
+instance Binary DV where 
+	put (DV v i) = do 
+		put v 
+		put i 
+	get = do 
+		v <- get 
+		i <- get 
+		return $ DV v i 
+
+instance Binary (v Double) => Binary (PrivateCPT v Double) where 
+    put (Table d m v) = do 
+    	putWord8 0 
+    	put d 
+    	put m 
+    	put v 
+    put (Scalar v) = do 
+    	putWord8 1 
+    	put v 
+    get = do 
+    	tag <- getWord8 
+    	case tag of 
+    		0 -> do 
+    			d <- get 
+    			m <- get 
+    			v <- get 
+    			return $ Table d m v 
+    		_ -> get >>= return . Scalar 
+
+instance Binary Vertex where 
+  put (Vertex v) = put v 
+  get = get >>= return . Vertex 
+
+instance Binary Edge where 
+	put (Edge va vb) = do 
+		put va 
+		put vb 
+	get = do 
+		va <- get
+		vb <- get 
+		return $ Edge va vb
+
+instance (Binary l, Binary e, Binary v) => Binary (SimpleGraph l e v) where 
+	put (SP e v n) = do 
+		put e 
+		put v 
+		put n 
+	get = do 
+		e <- get 
+		v <- get 
+		n <- get 
+		return $ SP e v n
+
+instance Binary DE where 
+    put (DE a b) = do 
+    	put a 
+    	put b 
+    get = do 
+    	a <- get 
+    	b <- get 
+    	return $ DE a b 
+
+instance Binary UE where 
+	put (UE a) = put a 
+	get = get >>= return . UE
diff --git a/Bayes/ImportExport/HuginNet.hs b/Bayes/ImportExport/HuginNet.hs
--- a/Bayes/ImportExport/HuginNet.hs
+++ b/Bayes/ImportExport/HuginNet.hs
@@ -14,6 +14,7 @@
 import Bayes
 import Bayes.PrivateTypes
 import Bayes.Factor.CPT(changeVariableOrder)
+import Bayes.BayesianNetwork
 
 --import Debug.Trace 
 
diff --git a/Bayes/InfluenceDiagram.hs b/Bayes/InfluenceDiagram.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/InfluenceDiagram.hs
@@ -0,0 +1,570 @@
+{-# LANGUAGE ViewPatterns #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{- | Tools to build influence diagrams
+
+-}
+module Bayes.InfluenceDiagram(
+  -- * Type
+    InfluenceDiagram
+  , DecisionFactor
+  , Instantiable(..)
+  , DEV 
+  , UV
+  , DV
+  , TDV
+  -- * Building 
+  , t 
+  , (~~)
+  , chance 
+  , decisionNode 
+  , utilityNode 
+  , proba 
+  , decision 
+  , utility
+  , cpt
+  , d 
+  , p
+  , noDependencies
+  -- * Solving
+  , decisionsOrder
+  , solveInfluenceDiagram
+  , runID
+  , policyNetwork
+  , decisionToInstantiation
+  -- * Testing 
+  , DVISet 
+  , DVI
+	) where 
+
+import Bayes
+import Bayes.PrivateTypes
+import Bayes.Network 
+import Data.Monoid
+import Bayes.Factor
+import Bayes.Factor.PrivateCPT
+import Bayes.Factor.CPT
+import Bayes.Factor.MaxCPT
+import Data.Maybe(fromJust,mapMaybe)
+import Control.Applicative((<$>))
+import Bayes.VariableElimination.Buckets
+import Bayes.Factor.PrivateCPT(DecisionFactor(..),decisionFactor,convertToMaxFactor,convertToNormalFactor,privateFactorValue,_factorVariables)
+import Data.List(foldl1',foldl')
+import Control.Monad.State.Strict(gets)
+import qualified Data.Vector as NV
+import qualified Data.Map as M 
+import qualified Data.IntMap as IM 
+
+--import Debug.Trace 
+
+--debug s a = trace (s ++ show a) a
+
+replaceDecisionNodeWithPolicy :: InfluenceDiagram -> CPT -> InfluenceDiagram 
+replaceDecisionNodeWithPolicy g f = 
+  let dv = factorMainVariable f 
+      parentVariables = tail (factorVariables f)
+      factorV = vertex dv
+      g' = fromJust $ changeVertexValue factorV (DecisionNode f) g
+      oldParentEdges = fromJust $ ingoing g' factorV 
+      g'' = foldr removeEdge g' oldParentEdges
+      addNewFactorEdge pdv currentG = addEdge (edge (vertex pdv) (vertex dv)) NormalLink currentG
+  in 
+  -- When we add the edges, we don't care about the order because
+  -- the CPT is already created and we don't need the order of the parents in the graph
+  -- to deduce the order of the variables in the CPT.
+  foldr addNewFactorEdge g'' parentVariables
+
+-- | Create a policy factor from a decision factor
+policyFactor :: DecisionFactor -> CPT
+policyFactor (Scalar v) = 
+  let decisionVar = v 
+      originalV = vertex decisionVar 
+      nakedVars = [dv decisionVar]
+      allVars = DVSet nakedVars
+      values = do 
+        x <- forAllInstantiations allVars 
+        if (instantiationValue v == instantiationValue (head x))
+          then 
+            return 1.0 
+          else 
+            return 0.0
+  in 
+  fromJust . factorWithVariables nakedVars $ values 
+policyFactor f@(Table d m v) = 
+  let decisionVar = NV.head v 
+      originalV = vertex decisionVar 
+      nakedVars = dv decisionVar : d
+      allVars = DVSet nakedVars
+      values = do 
+        x <- forAllInstantiations allVars 
+        let v = privateFactorValue f (tail x)
+        if (instantiationValue v == instantiationValue (head x))
+          then 
+            return 1.0 
+          else 
+            return 0.0
+  in
+  fromJust . factorWithVariables nakedVars $ values
+  
+
+-- | Convert a decision policy to a set of possible instantiations
+-- It is the only way to access to the content of a decision factor.
+decisionToInstantiation :: DecisionFactor -> [DVISet]
+decisionToInstantiation f@(Scalar v) = [[v]]
+decisionToInstantiation f@(Table d m v) = 
+  let allVars = DVSet d
+      values = do 
+        x <- forAllInstantiations allVars 
+        let v = privateFactorValue f x 
+        return (v:x)
+  in 
+  values
+
+-- | Create a policy network from an influence diagram and its solution.
+-- A policy network is a Bayesian network where the decision nodes have been replaced
+-- with probability nodes where the probability is 1 when the configuration is corresponding
+-- to the decision and 0 otherwise.
+policyNetwork :: [DecisionFactor] -> InfluenceDiagram -> SBN CPT 
+policyNetwork l idg = 
+  let idg1 = foldl' replaceDecisionNodeWithPolicy idg (map policyFactor l) 
+      utilities = filter (isUtilityNode idg) . allVertices $ idg
+      SP e v b = foldr removeVertex idg1 utilities
+      toBayesNode (l,ChanceNode f) = (l,f) 
+      toBayesNode (l,DecisionNode f) = (l,f)
+      toBayesNode (l,UtilityNode _ _) = error "No utilities nodes should remain to create the policy network"
+      e' = M.map (const ()) e 
+      v' = IM.map toBayesNode v 
+  in 
+  SP e' v' b
+
+
+instance Show DecisionFactor where
+    show (Scalar v) = "\nScalar Factor:\n" ++ show v
+    show c@(Table [] _ v) = "\nEmpty DecisionFactor:\n"
+
+    show c = displayFactorBody c 
+
+instance MultiDimTable DecisionFactor where 
+    elementStringValue f d = show (privateFactorValue f d)
+    tableVariables = _factorVariables
+
+
+data JoinSum = JS !CPT !CPT deriving(Eq)
+
+instance Show JoinSum where 
+  show (JS p u) = "CPT\n" ++ show p ++ "\nUTILITY\n" ++ show u ++ "\n"
+
+chanceFactor f = JS f (factorFromScalar 0.0)
+utilityFactor f = JS (factorFromScalar 1.0) f
+
+jsProduct ::  JoinSum -> JoinSum -> JoinSum
+jsProduct (JS pa ua) (JS pb ub) = JS (itemProduct [pa,pb]) (cptSum [ua,ub])
+
+
+-- | Max out a variable
+maximalize :: DV -> [JoinSum] -> (JoinSum,DecisionFactor) 
+maximalize dv l = 
+  let JS pa ua = itemProduct l
+      maxa = convertToNormalFactor . itemProjectOut dv . convertToMaxFactor $ pa
+      maxu' = itemProjectOut dv . convertToMaxFactor  . itemProduct $ [pa,ua]
+      maxu = convertToNormalFactor maxu'
+      instF = decisionFactor maxu'
+  in 
+  (JS maxa (cptDivide maxu maxa),instF)
+
+
+
+instance IsBucketItem JoinSum where
+    scalarItem (JS a b) = isScalarFactor a && isScalarFactor b
+    itemProduct l = foldl1' jsProduct l
+    itemProjectOut dv (JS pa ua) =  
+      let suma = itemProjectOut dv pa
+          sumu = itemProjectOut dv (itemProduct [pa,ua])
+      in 
+      JS suma (cptDivide sumu suma)
+    itemContainsVariable (JS a b) dv = containsVariable a dv || containsVariable b dv
+
+
+-- | Synonym for undefined because it is clearer to use t to set the Enum bounds of a variable
+t = undefined
+
+-- | Edge kind
+data EdgeKind = NormalLink 
+              deriving(Eq,Show)
+
+isInformationLink :: InfluenceDiagram -> Edge -> Bool
+isInformationLink g (Edge va vb) = 
+  (isChanceNode g va || isDecisionNode g va) && (isDecisionNode g vb)
+
+isRevealedChanceNode :: InfluenceDiagram -> Vertex -> Bool 
+isRevealedChanceNode g v = isChanceNode g v && any (isDecisionNode g) (childrenNodes g v)
+
+edgeShape :: InfluenceDiagram -> Edge -> EdgeKind -> Maybe String
+edgeShape g e NormalLink | isInformationLink g e = Just "style=dashed"
+                         | otherwise = Nothing 
+edgeColor :: InfluenceDiagram -> Edge -> EdgeKind -> Maybe String
+edgeColor _ _ _ = Nothing
+
+nodeShape :: InfluenceDiagram -> Vertex -> IDValue -> Maybe String
+nodeShape _ _ (ChanceNode _) = Just "shape=ellipse" 
+nodeShape _ _ (UtilityNode _ _) = Just "shape=diamond"
+nodeShape _ _ (DecisionNode _) = Just "shape=box"
+
+nodeColor :: InfluenceDiagram -> Vertex -> IDValue -> Maybe String
+nodeColor g v (ChanceNode _) | isRevealedChanceNode g v = Just "style=filled,fillcolor=gray"
+                             | otherwise = Nothing
+nodeColor _ _ _ = Nothing
+
+instance Show InfluenceDiagram where
+  show g = displaySimpleGraph (nodeShape g) (nodeColor g) (edgeShape g) (edgeColor g) g
+
+
+instance Monoid EdgeKind where 
+   mempty = NormalLink 
+   NormalLink `mappend` NormalLink = NormalLink
+
+-- | Influence diagram
+type InfluenceDiagram = DirectedSG EdgeKind IDValue
+
+type IDMonad g a = NetworkMonad g EdgeKind IDValue a
+
+
+-- Most factors are coding for f(abc) where a is the main factor variable
+-- and where the a vertex is the original vertex in the graph.
+-- For an utility, we have U(abcd) where the variables a is NOT the main factor variable.
+-- Indeed, the main factor variable if it was used in the factor would have dimension 1.
+-- So, it is useless in the factor.
+-- So, we need another DV field to track the original vertex
+data IDValue   = ChanceNode !CPT
+               | UtilityNode !DV !CPT 
+               | DecisionNode !CPT
+               deriving(Eq)
+
+dvFromIDValue (ChanceNode f) = factorMainVariable f
+dvFromIDValue (UtilityNode dv f) = dv
+dvFromIDValue (DecisionNode f) = factorMainVariable f
+
+factorVariablesFromIDValue (ChanceNode f) = factorVariables f
+factorVariablesFromIDValue (UtilityNode _ f) = factorVariables f
+factorVariablesFromIDValue (DecisionNode f) = factorVariables f
+
+jsFromIDValue (ChanceNode f) = chanceFactor f
+jsFromIDValue (UtilityNode _ f) = utilityFactor f
+jsFromIDValue (DecisionNode _) = error "You don't need to get the factor for a decision node"
+
+instance Show IDValue where 
+   show (ChanceNode f) = "CHANCE:\n" ++ show f
+   show (UtilityNode _ f) = "UTILITY:\n" ++ show f
+   show (DecisionNode f) = ""
+
+
+-- | Utility variable
+data UV = UV !Vertex !Int deriving(Eq)
+
+-- | Decision variable
+data DEV = DEV !Vertex !Int deriving(Eq,Ord)
+
+instance Show DEV where
+    show (DEV v d) = "D" ++ show v ++ "(" ++ show d ++ ")"
+
+instance BayesianDiscreteVariable DEV where 
+  dimension (DEV _ d) = d 
+  dv (DEV v d) = DV v d 
+  vertex (DEV v _) = v
+
+instance Instantiable DEV Int where 
+  (=:) d@(DEV v dim) value = DVI (dv d) value
+
+data PorD = P DV | D DEV deriving(Eq)
+
+class ChanceVariable m where 
+  toDV :: m -> DV
+
+instance ChanceVariable DV where 
+  toDV = dv 
+
+instance ChanceVariable (TDV s) where 
+  toDV = dv
+
+instance BayesianDiscreteVariable PorD where
+    dimension (D d) = dimension d
+    dimension (P p) = dimension p
+    dv (D x) = dv x
+    dv (P x) = dv x
+    vertex (D d) = vertex d
+    vertex (P p) = vertex p
+
+-- | Used to mix decision and chance variables and a same list
+p :: ChanceVariable c => c -> PorD
+p = P . toDV
+
+-- | Used to mix decision and chance variables and a same list
+d :: DEV -> PorD 
+d = D
+
+
+
+-- | Create a chance node
+chance :: (Bounded a, Enum a, NamedGraph g)
+       => String 
+       -> a 
+       -> IDMonad g (TDV a)
+chance = variable
+
+-- | Create an utility node
+utilityNode :: (NamedGraph g)
+            => String 
+            -> IDMonad g UV
+utilityNode s = do
+  DV v i <- variableWithSize s 1
+  return (UV v i)
+
+-- | Create a decision node
+decisionNode :: (Bounded a, Enum a, NamedGraph g)
+             => String 
+             -> a
+             -> IDMonad g DEV
+decisionNode s a =  do
+  DV v i <- variable s a >>= return . dv
+  return (DEV v i)
+
+utilityCpt :: (DirectedGraph g, Distribution d, Factor f) 
+           => Vertex -- ^ Vertex containing the factor
+           -> d -- ^ Distribution to initialize the factor
+           -> NetworkMonad g e a (Maybe f) 
+utilityCpt v l  = do 
+  g <- gets snd
+  let vertices = map (fromJust . startVertex g) . fromJust . ingoing g $ v
+  fv <- mapM factorVariable vertices
+  let cpt = createFactor (map fromJust fv) l
+  return cpt
+
+class Initializable v where 
+  (~~) :: (DirectedGraph g, Distribution d) 
+     => IDMonad g v -- ^ Discrete variable in the graph
+     -> d -- ^ List of values
+     -> IDMonad g ()
+
+instance Initializable DV where
+  (~~) mv l = do 
+     (DV v _) <- mv >>= return . dv -- This is updating the state and so the graph
+     maybeNewValue <- getCpt v l
+     currentValue <- getBayesianNode v
+     case (currentValue, maybeNewValue) of 
+       (Just c, Just n) -> initializeNodeWithValue v c (ChanceNode n)
+       _ -> return ()
+
+instance Initializable (TDV s) where
+  (~~) mv l = do 
+     (DV v _) <- mv >>= return . dv -- This is updating the state and so the graph
+     maybeNewValue <- getCpt v l
+     currentValue <- getBayesianNode v
+     case (currentValue, maybeNewValue) of 
+       (Just c, Just n) -> initializeNodeWithValue v c (ChanceNode n)
+       _ -> return ()
+
+instance Initializable UV where
+  (~~) mv l = do 
+     (UV v dim) <- mv  -- This is updating the state and so the graph
+     maybeNewValue <- utilityCpt v l
+     currentValue <- getBayesianNode v
+     case (currentValue, maybeNewValue) of 
+       (Just c, Just n) -> initializeNodeWithValue v c (UtilityNode (DV v dim) n)
+       _ -> return ()
+
+instance Initializable DEV where
+  (~~) mv l = do 
+     (DV v _) <- mv >>= return . dv -- This is updating the state and so the graph
+     maybeNewValue <- getCpt v l
+     currentValue <- getBayesianNode v
+     case (currentValue, maybeNewValue) of 
+       (Just c, Just n) -> initializeNodeWithValue v c (DecisionNode n)
+       _ -> return ()
+
+_cpt :: (DirectedGraph g , BayesianDiscreteVariable v,BayesianDiscreteVariable vb) => v -> [vb] -> IDMonad g v
+_cpt node conditions = do
+  mapM_ ((dv node) <--) (reverse (map dv conditions))
+  return node
+
+-- | Define that a chance node is a conditional probability and define the parent variables
+cpt :: (DirectedGraph g ,BayesianDiscreteVariable vb, ChanceVariable c) => c -> [vb] -> IDMonad g c
+cpt node conditions = do
+  mapM_ ((toDV node) <--) (reverse (map dv conditions))
+  return node
+
+-- | Define that a chance node is a probability (not conditional)
+-- Values are ordered like
+-- FFF FFT FTF FTT TFF TFT TTF TTT
+-- and same for other enumeration keeping enumeration order
+proba :: (ChanceVariable c, DirectedGraph g) => c -> IDMonad g c
+proba node = cpt node ([] :: [DV])
+
+-- | Define a utility dependence
+utility :: (DirectedGraph g , BayesianDiscreteVariable dv) => UV -> [dv] -> IDMonad g UV
+utility (UV v d) l = do 
+  DV v' d' <- _cpt (DV v d) l
+  return (UV v' d')
+
+-- | Used to define a root decision which is not dependent on any past node
+noDependencies :: [DV]
+noDependencies = []
+
+-- | Define a decision dependence
+decision :: (DirectedGraph g, BayesianDiscreteVariable dv) => DEV -> [dv] -> IDMonad g DEV
+decision d l = do 
+  let dim = product . map dimension $ dv d:map dv l
+  _cpt d l ~~ (replicate dim 1.0)
+  return d
+
+
+-- | Run an influence monad
+runID :: IDMonad DirectedSG a -> (a,InfluenceDiagram)
+runID = runNetwork 
+
+{-
+ 
+Generation of temporal order
+
+-}
+
+maybeOnlyResult :: [a] -> Maybe a 
+maybeOnlyResult [a] = Just a 
+maybeOnlyResult _ = Nothing
+
+isDecisionNode :: InfluenceDiagram -> Vertex -> Bool 
+{-# INLINE isDecisionNode #-}
+isDecisionNode g v = maybe False (const True) $ do
+  DecisionNode f <- vertexValue g v
+  return f
+
+isUtilityNode :: InfluenceDiagram -> Vertex -> Bool 
+{-# INLINE isUtilityNode #-}
+isUtilityNode g v = maybe False (const True) $ do
+  UtilityNode _ f <- vertexValue g v
+  return f
+
+isChanceNode :: InfluenceDiagram -> Vertex -> Bool 
+{-# INLINE isChanceNode #-}
+isChanceNode g v = maybe False (const True) $ do
+  ChanceNode f <- vertexValue g v
+  return f
+
+-- | Check the node is a decision node and none of its parents are decision nodes
+isRootDecision :: InfluenceDiagram -> Vertex -> Bool
+{-# INLINE isRootDecision #-}
+isRootDecision g v | isDecisionNode g v = 
+  case ingoing g v of 
+    Just [] -> True 
+    _ -> False
+                   | otherwise = False
+
+-- | Get the chance parents of a decision node and a new graph with those parents
+-- removed and the decision node removed
+chanceParents :: DEV -> InfluenceDiagram -> (InfluenceDiagram,[DV])
+chanceParents dev currentG = 
+  let p = filter (isChanceNode currentG) . parentNodes currentG $ (vertex dev) 
+      theParents = map (vertexToDV currentG) p
+      newG = foldr removeVertex currentG (vertex dev : p)
+  in 
+  (newG,theParents)
+
+-- | Return the remaining chance nodes of the graph
+remainingChanceNodes :: InfluenceDiagram -> [DV]
+remainingChanceNodes = chanceNodes 
+
+-- | Return the utility nodes
+utilityNodes :: InfluenceDiagram -> [UV]
+utilityNodes g = map (vertexToUV g) . filter (isUtilityNode g) . allVertices $ g
+
+-- | Return the chance nodes
+chanceNodes :: InfluenceDiagram -> [DV]
+chanceNodes g = map (vertexToDV g) . filter (isChanceNode g) . allVertices $ g
+
+-- | Return all chances factor and utility factors
+chanceAndUtilityFactors :: InfluenceDiagram -> [JoinSum]
+chanceAndUtilityFactors g = map (jsFromIDValue . fromJust . vertexValue g) . filter (not . isDecisionNode g) . allVertices $ g
+
+vertexToDV :: InfluenceDiagram -> Vertex -> DV 
+vertexToDV g v = dvFromIDValue . fromJust . vertexValue g $ v
+
+vertexToDEV :: InfluenceDiagram -> Vertex -> DEV 
+vertexToDEV g v = 
+  let DV v1 d = vertexToDV g v 
+  in 
+  DEV v1 d
+
+vertexToUV :: InfluenceDiagram -> Vertex -> UV 
+vertexToUV g v = 
+  let DV v1 d = vertexToDV g v 
+  in 
+  UV v1 d
+
+-- | Return a root decision
+rootDecision :: InfluenceDiagram -> Maybe Vertex 
+rootDecision g = do 
+    r <- rootNode g 
+    if isDecisionNode g r 
+      then
+        return r 
+      else 
+        rootDecision (removeVertex r g)
+
+-- | Used to encode the temporal order
+data ChancesOrDecision = C ![DV] | DEC !DEV deriving(Eq,Ord,Show)
+
+dvOrder :: [ChancesOrDecision] -> [DV] 
+dvOrder [] = []
+dvOrder (C l:r) = l ++ dvOrder r 
+dvOrder (DEC d:r) = dv d: dvOrder r 
+
+-- | Remove all decisions node and record their chance parents
+removeAndRecordRootDecision :: [ChancesOrDecision] -> InfluenceDiagram -> [ChancesOrDecision]
+removeAndRecordRootDecision currentL currentG = 
+  case vertexToDEV currentG <$> (rootDecision currentG) of 
+    Nothing -> (C (remainingChanceNodes currentG)):currentL 
+    Just newD -> 
+      let (currentG', p) = chanceParents newD currentG
+      in
+      removeAndRecordRootDecision ((DEC newD):(C p):currentL) currentG'
+
+-- | List of decision vertices in reverse temporal order (corresponding to elimination order)
+decisionsOrder :: InfluenceDiagram -> [ChancesOrDecision] 
+decisionsOrder g = removeAndRecordRootDecision [] $ g 
+
+-- | Maximalize the decision variable to find the decision strategy at the current set
+maximalizeOneVariable :: Buckets JoinSum -> DV -> (Buckets JoinSum,DecisionFactor)
+maximalizeOneVariable currentBucket dv   = 
+  let fk = getBucket dv currentBucket
+      (newF, instF) = maximalize dv fk
+  in
+  (updateBucket dv newF currentBucket, instF)
+
+--debugs True f = \a b -> debug ("SUM OUT " ++ show b ++ "\n") (f a b)
+--debugs False f = \a b -> debug ("MAX OUT " ++ show b ++ "\n") (f a b)
+
+marginalizeID :: [ChancesOrDecision]  -> Buckets JoinSum -> [DecisionFactor] -> (Buckets JoinSum,[DecisionFactor])
+marginalizeID [] b r = (b,r)
+marginalizeID (C d:r) currentB currentR =  
+  let bucket' = foldl' marginalizeOneVariable currentB d 
+  in 
+  marginalizeID r bucket' currentR 
+marginalizeID (DEC de:r) currentB currentR = 
+  let (bucket',instF) = maximalizeOneVariable currentB (dv de) 
+  in 
+  marginalizeID r bucket' (instF:currentR) 
+
+-- | Solve an influence diagram. A DecisionFactor is generated for each decision variable.
+-- A decision factor is containing a variable instantiation instead of a double.
+-- This instantiation is giving the decision to take for each value of the parents.
+solveInfluenceDiagram :: InfluenceDiagram -> [DecisionFactor]
+solveInfluenceDiagram g = 
+  let decOrder = decisionsOrder g
+      theFactors = chanceAndUtilityFactors g
+      p = dvOrder decOrder 
+      bucket = createBuckets theFactors p []
+      (_, result) = marginalizeID decOrder bucket []
+  in
+  result 
+
+
diff --git a/Bayes/Network.hs b/Bayes/Network.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Network.hs
@@ -0,0 +1,214 @@
+{-# LANGUAGE ViewPatterns #-}
+{- | Common functions and types for building networks
+
+-}
+module Bayes.Network(
+	-- * Types 
+	  MaybeNode(..)
+	, NetworkMonad(..)
+	-- * Functions
+	, factorVariable
+	, (<--)
+	, getBayesianNode 
+	, setBayesianNode
+	, initializeNodeWithValue
+	, setVariableBoundWithSize
+	, setVariableBound
+	, addVariableIfNotFound
+	, unamedVariable
+	, variable
+	, variableWithSize
+	, unNamedVariableWithSize
+	, runNetwork
+	, execNetwork
+	, evalNetwork
+	, getCpt
+	) where 
+
+import Bayes.PrivateTypes
+import Bayes 
+import Control.Monad.State.Strict
+import Bayes.Tools 
+import Data.Maybe(fromJust)
+import Bayes.Factor
+import Data.Monoid
+
+-- | Bayesian variable : name,dimension, factor
+-- When initialized it is using a factor with bayesian variables.
+-- But the factor value are not yet set
+data MaybeNode f = UninitializedNode String Int
+                 | InitializedNode String Int f
+
+
+-- | The Network monad
+type NetworkMonad g e f a = GraphMonad g e (MaybeNode f) a
+
+
+-- | Get the Bayesian Discrete Variable for a vertex.
+-- It works because we keep the variable dimension during creating of the graph
+factorVariable :: Graph g => Vertex -> NetworkMonad g e f (Maybe DV)  
+factorVariable v = do 
+  g <- gets snd 
+  let value = vertexValue g v
+  case value of
+    Nothing -> return Nothing
+    Just (UninitializedNode _ d) -> return $ Just $ DV v d
+    Just (InitializedNode _ d _) -> return $ Just $ DV v d
+
+
+-- | Create an edge between two vertex of the Bayesian network
+(<--) :: (Graph g, BayesianDiscreteVariable dv, Monoid e) => dv -> dv -> NetworkMonad g e f ()
+(dv -> DV va _) <-- (dv -> DV vb _) = newEdge vb va mempty
+
+whenJust Nothing _ = return ()
+whenJust (Just i) f = f i >> return ()
+
+getCpt :: (DirectedGraph g, Distribution d, Factor f) 
+       => Vertex -- ^ Vertex containing the factor
+       -> d -- ^ Distribution to initialize the factor
+       -> NetworkMonad g e a (Maybe f) 
+getCpt v l  = do 
+  g <- gets snd
+  currentVar <- factorVariable v
+  let vertices = map (fromJust . startVertex g) . fromJust . ingoing g $ v
+  fv <- mapM factorVariable vertices
+  let cpt = createFactor (map fromJust (currentVar:fv)) l
+  return cpt
+
+-- | Get the node of a bayesian network under creation
+getBayesianNode :: Graph g => Vertex -> NetworkMonad g e f (Maybe (MaybeNode f))
+getBayesianNode v = do
+  g <- gets snd
+  return $ vertexValue g v
+
+-- | Set the node of a bayesian network under creation
+setBayesianNode :: Graph g => Vertex -> MaybeNode f -> NetworkMonad g e f ()
+setBayesianNode v newValue = do
+  (aux,oldGraph) <- get
+  let newGraph = changeVertexValue v newValue oldGraph
+ 
+  whenJust newGraph $ \nvm -> do
+     put $! (aux, nvm)
+
+-- | Set the value of uninitialized nodes. Initialized nodes are not changed.
+initializeNodeWithValue :: Graph g 
+                        => Vertex -- ^ Vertex
+                        -> MaybeNode a -- ^ Current uninitialized node
+                        -> a -- ^ Value to set
+                        -> NetworkMonad g e a () 
+initializeNodeWithValue _ (InitializedNode _ _ _) _ = return ()
+initializeNodeWithValue v (UninitializedNode s dim) newValue = do 
+  g <- gets snd
+  setBayesianNode v (InitializedNode s dim newValue)
+
+-- | Set the bound of a bayesian variable (number of levels)
+setVariableBoundWithSize :: Graph g
+                         => Vertex -- ^ Vertex
+                         -> Int -- ^ Inf limit  (0 for instance)
+                         -> Int -- ^ Sup limit (1 for instance for 2 elements)
+                         -> NetworkMonad g e f ()
+setVariableBoundWithSize a bmin bmax = do
+    v <- getBayesianNode a
+    whenJust v $ \(UninitializedNode s _) -> do
+      setBayesianNode a (UninitializedNode s (bmax - bmin + 1))
+
+setVariableBound :: (Enum a, Bounded a, Graph g) 
+                 => Vertex -- ^ Vertex
+                 -> a -- ^ Bounded variable (t :: type where t is undefined)
+                 -> NetworkMonad g e f ()
+setVariableBound a e = 
+  let bmin = intValue $ minBoundForEnum e
+      bmax = intValue $ maxBoundForEnum e
+  in 
+  setVariableBoundWithSize a bmin bmax
+
+-- | Create a new named Bayesian variable if not found.
+-- Otherwise, return the found one.
+addVariableIfNotFound :: NamedGraph g => String -> NetworkMonad g e f Vertex
+addVariableIfNotFound vertexName = graphNode vertexName (UninitializedNode vertexName 0)
+
+-- | Initialize a new variable
+_initializeVariableBounds :: (Enum a, Bounded a, NamedGraph g)
+                          => Vertex 
+                          -> a 
+                          -> NetworkMonad g e f (TDV a)
+_initializeVariableBounds va e = do 
+  setVariableBound va e
+  maybeValue <- getBayesianNode va 
+  case fromJust maybeValue of 
+     UninitializedNode s d -> return (tdv $ DV va d)
+     InitializedNode _ d _ -> return (tdv $ DV va d) 
+
+-- | Initialize a new variable with size
+_initializeVariableBoundsWithSize :: NamedGraph g
+                                  => Vertex -- ^ Variable name
+                                  -> Int -- ^ Variable size
+                                  -> NetworkMonad g e f DV
+_initializeVariableBoundsWithSize va e = do
+  setVariableBoundWithSize va 0 (e-1)
+  maybeValue <- getBayesianNode va 
+  setBayesianNode va (fromJust maybeValue)
+  case fromJust maybeValue of 
+     UninitializedNode s d -> return (DV va d)
+     InitializedNode _ d _ -> return (DV va d)
+
+-- | Create a new unamed variable
+unamedVariable :: (Enum a, Bounded a, NamedGraph g)
+               => a -- ^ Variable bounds 
+               -> NetworkMonad g e f (TDV a)
+unamedVariable e = do 
+  va <- getNewEmptyVariable Nothing (UninitializedNode "unamed" 0)
+  _initializeVariableBounds va e
+
+-- | Define a Bayesian variable (name and bounds)
+variable :: (Enum a, Bounded a, NamedGraph g) 
+        => String -- ^ Variable name
+        -> a -- ^ Variable bounds
+        -> NetworkMonad g e f (TDV a)
+variable name e = do
+  va <- addVariableIfNotFound name
+  _initializeVariableBounds va e
+
+-- | Define a Bayesian variable (name and bounds)
+variableWithSize :: NamedGraph g
+        => String -- ^ Variable name
+        -> Int -- ^ Variable size
+        -> NetworkMonad g e f DV
+variableWithSize name e = do
+  va <- addVariableIfNotFound name
+  _initializeVariableBoundsWithSize va e
+
+-- | Define a Bayesian variable (name and bounds)
+unNamedVariableWithSize :: NamedGraph g
+                        => Int -- ^ Variable size
+                        -> NetworkMonad g e f DV
+unNamedVariableWithSize e = do
+  va <- getNewEmptyVariable Nothing (UninitializedNode "unamed" 0)
+  _initializeVariableBoundsWithSize va e
+
+-- | Create a  network using the simple graph implementation
+-- The initialized nodes are replaced by the value.
+-- Returns the monad values and the built graph.
+runNetwork :: NetworkMonad DirectedSG e f a -> (a,DirectedSG e f)
+runNetwork x = 
+  let (r,g) = runGraph x
+      convertNodes (InitializedNode s d f) = f 
+      convertNodes (UninitializedNode s d) = error $ "All variables must be initialized with a factor: " ++ s ++ "(" ++ show d ++ ")"
+  in 
+  (r,fmap convertNodes g)
+
+-- | Create a  network but only returns the monad value.
+-- Mainly used for testing.
+execNetwork :: NetworkMonad DirectedSG e f a -> DirectedSG e f
+execNetwork x = 
+  let g = execGraph x
+      convertNodes (InitializedNode s d f) = f 
+      convertNodes (UninitializedNode s d) = error $ "All variables must be initialized with a factor: " ++ s ++ "(" ++ show d ++ ")"
+  in 
+  fmap convertNodes g
+
+
+-- | Create a bayesian network but only returns the monad value.
+-- Mainly used for testing.
+evalNetwork :: Graph g => NetworkMonad g e f a -> a
+evalNetwork = evalGraph
diff --git a/Bayes/PrivateTypes.hs b/Bayes/PrivateTypes.hs
--- a/Bayes/PrivateTypes.hs
+++ b/Bayes/PrivateTypes.hs
@@ -23,8 +23,12 @@
  , instantiationValue
  , instantiationVariable
  , fromDVSet
- -- * Vertices 
+ -- * Vertices, Graph
  , Vertex(..)
+ , Edge(..)
+ , SimpleGraph(..)
+ , DE(..)
+ , UE(..)
  -- * Misc
  , getMinBound
  -- * Indices 
@@ -45,6 +49,7 @@
 import Test.QuickCheck.Arbitrary
 import System.Random(Random)
 import qualified Data.IntMap as IM
+import qualified Data.Map as M
 
 
 
@@ -97,6 +102,28 @@
 -}
 -- | Vertex type used to identify a vertex in a graph
 newtype Vertex = Vertex {vertexId :: Int} deriving(Eq,Ord)
+
+-- | Edge type used to identify and edge in a graph
+data Edge = Edge !Vertex !Vertex deriving(Eq,Ord,Show)
+
+-- | Implementtaion of a SimpleGraph
+data SimpleGraph local edgedata vertexdata = SP {
+ -- | Mapping of edge to edge data
+   edgeMap :: !(M.Map Edge edgedata) 
+ -- ^ Mapping of vertex number to vertex neighborhood and vertex data
+ , vertexMap :: !(IM.IntMap (local, vertexdata))
+ -- ^ Vertex names. Used only to generate the graphviz representation. Names are useless for the algorithms
+ -- and I don't want them to appear in the vertex values which should only be factor. Otherwise, the algorithms
+ -- are less elegant since I have to extract the factors from the values
+ , nameMap :: !(IM.IntMap String)
+ } 
+
+-- | Neighborhood structure for directed or undirected edges
+-- | Directed edges
+data DE = DE ![Edge] ![Edge] deriving(Eq,Show)
+
+-- | Undirected edges
+data UE = UE ![Edge] deriving(Eq,Show)
 
 instance Show Vertex where 
     show (Vertex v) = "v" ++ show v
diff --git a/Bayes/Test.hs b/Bayes/Test.hs
--- a/Bayes/Test.hs
+++ b/Bayes/Test.hs
@@ -1,3 +1,6 @@
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 {- | Testing of the implementation.
 
 -}
@@ -13,10 +16,12 @@
 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)
-
+import Bayes.Test.InfluencePatterns(testStudentDecisions)
 #ifdef LOCAL
-import Bayes.Test.ReferencePatterns(compareAsiaReference,compareCancerReference,comparePokerReference,compareFarmReference,compareMpeCancer)
-#endif 
+import Bayes.Test.ReferencePatterns(compareAsiaReference,compareCancerReference,comparePokerReference,compareFarmReference,compareMpeAsia,testFileExport)
+#else 
+import Bayes.Test.ReferencePatterns(testFileExport)
+#endif
 
 -- | Run all the tests
 runTests = defaultMain tests
@@ -41,15 +46,19 @@
                 testProperty "Test all clusters are included in the junction tree" junctionTreeAllClusters_prop
             ]
         , testGroup "Misc functions" [
-                testProperty "Instantiation from multiindex" instantiationProp
+                testProperty "Instantiation from multiindex" instantiationProp,
+                testCase "Test import/export of bayesian network and junction tree" testFileExport
             ]
+        , testGroup "Influence Diagrams" [
+                testCase "Test with reference patterns" testStudentDecisions
+            ]
 #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 "Test MPE and MAP with Cancer network" compareMpeCancer
+                testCase "Test MPE and MAP with Asia network" compareMpeAsia
         ]
 #endif
 
diff --git a/Bayes/Test/InfluencePatterns.hs b/Bayes/Test/InfluencePatterns.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Test/InfluencePatterns.hs
@@ -0,0 +1,24 @@
+{- | A comparison of influence diagram solution with references
+
+-}
+module Bayes.Test.InfluencePatterns(
+   testStudentDecisions
+ ) where
+
+import Test.HUnit.Base(assertBool)
+import Data.Maybe(fromJust)
+
+import Bayes.Examples.Influence
+import Bayes.InfluenceDiagram 
+
+testStudentDecisions = do 
+	let result = solveInfluenceDiagram student
+	    (e,pr,s) = studentDecisionVars
+	    l = map decisionToInstantiation result 
+	assertBool "Student Network" $ l == [[[e =: (0::Int)]],[[s =: (0::Int),pr =: False],[s =: (0::Int),pr =: True]]]
+	let result = solveInfluenceDiagram studentSimple
+	    e = studentSimpleDecisionVar
+	    l = map decisionToInstantiation result 
+	assertBool "Simple Student Network" $ l == [[[ e =: (1 :: Int)]]]
+
+
diff --git a/Bayes/Test/ReferencePatterns.hs b/Bayes/Test/ReferencePatterns.hs
--- a/Bayes/Test/ReferencePatterns.hs
+++ b/Bayes/Test/ReferencePatterns.hs
@@ -5,12 +5,13 @@
 
 -}
 module Bayes.Test.ReferencePatterns(
+   testFileExport
 #ifdef LOCAL
-   compareAsiaReference
+ , compareAsiaReference
  , compareCancerReference
  , comparePokerReference
  , compareFarmReference
- , compareMpeCancer
+ , compareMpeAsia
 #endif
  ) where
 
@@ -21,8 +22,11 @@
 import Bayes
 import Bayes.FactorElimination
 import Bayes.VariableElimination(mpe)
-import Bayes.Examples(anyExample)
+import Bayes.Examples(anyExample,example)
 import Bayes.FactorElimination.JTree(root)
+import Bayes.Tools(withTempFile)
+import Bayes.ImportExport 
+import Bayes.BayesianNetwork
 
 value varmap jt s = 
   let v =  fromJust $ Map.lookup s varmap
@@ -37,6 +41,23 @@
   putStrLn ""
   assertBool s $ r ~=~ l
 
+-- | Test that we can import / export the bayesian network, junction tree and variable map
+testFileExport :: IO () 
+testFileExport = do 
+  let (vars,g) = example 
+      vm = varMap g
+      jt = createJunctionTree nodeComparisonForTriangulation g
+  withTempFile $ \f -> do 
+    writeNetworkToFile f g 
+    g' <- readNetworkFromFile f 
+    assertBool "Test graph import/export" $ g == g'
+  withTempFile $ \f -> do 
+    writeVariableMapAndJunctionTreeToFile f vm jt 
+    (vm',jt') <- readVariableMapAndJunctionTreeToFile f 
+    assertBool "Test jt import/export" $ jt == jt'
+    assertBool "Test variable map import/export" $ vm == vm'
+     
+
 -- Check that the float values are equal with an accuracy < 0.01%
 comparePercent :: Double -> Double -> Bool
 comparePercent a b = abs (a-b) < 1e-4
@@ -52,7 +73,7 @@
 rename g = \(a,s) -> (fromJust . vertexLabel g . vertex $ a, s)
 
 -- | Test that a MAP is not always the projection of a MPE
-compareMpeCancer = do 
+compareMpeAsia = 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]
diff --git a/Bayes/Tools.hs b/Bayes/Tools.hs
--- a/Bayes/Tools.hs
+++ b/Bayes/Tools.hs
@@ -2,10 +2,16 @@
 	
 -}
 module Bayes.Tools (
-	nearlyEqual
+	  nearlyEqual
+  , withTempFile
+  , minBoundForEnum 
+  , maxBoundForEnum 
+  , intValue
 	) where 
 
-
+import System.IO(openTempFile,hClose)
+import System.Directory(getTemporaryDirectory,removeFile)
+import System.FilePath
 
 -- | Floating point number comparisons which should take into account
 -- all the subtleties of that kind of comparison
@@ -20,3 +26,26 @@
         (x,y) | x == y -> True -- handle infinities
               | x*y == 0 -> diff < (epsilon * epsilon)
               | otherwise -> diff / (absA + absB) < epsilon
+
+
+-- | Execute an action with a temporary file. The file is deleted after.
+-- The action must close the file.
+-- (would be better to use handle to force the closing but it is used with action which are
+-- using a filepath)
+withTempFile :: (FilePath -> IO a) -> IO a 
+withTempFile action = do 
+  tempDir <- getTemporaryDirectory 
+  (filePath,fileHandle) <- openTempFile tempDir "bayestest" 
+  hClose fileHandle 
+  result <- action filePath 
+  removeFile (tempDir </> filePath)
+  return result 
+
+minBoundForEnum :: Bounded a => a -> a
+minBoundForEnum _ = minBound
+
+maxBoundForEnum :: Bounded a => a -> a
+maxBoundForEnum _ = maxBound
+
+intValue :: Enum a => a -> Int
+intValue = fromEnum
diff --git a/Bayes/VariableElimination.hs b/Bayes/VariableElimination.hs
--- a/Bayes/VariableElimination.hs
+++ b/Bayes/VariableElimination.hs
@@ -19,7 +19,7 @@
 
 import Bayes
 import Bayes.Factor
-import Data.List(partition,minimumBy,(\\),find,foldl')
+import Data.List(minimumBy,(\\),foldl')
 import Data.Maybe(fromJust)
 import Data.Function(on)
 import qualified Data.Map as M
@@ -27,14 +27,14 @@
 import Bayes.Factor.CPT 
 import Bayes.Factor.MaxCPT
 import Bayes.PrivateTypes(DVISet)
+import Bayes.VariableElimination.Buckets
 
 --import Debug.Trace 
 
 --debug s a = trace (s  ++ "\n" ++ show a ++ "\n") a
 
--- | Elimination order
-type EliminationOrder dv = [dv]
 
+
 -- | Get all variables from a Bayesian Network
 allVariables :: (Graph g, Factor f) 
              => BayesianNetwork g f 
@@ -45,86 +45,15 @@
   in 
   map createDV s
 
--- | Used for bucket elimination. Factor are organized by their first DV
-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 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
-      -- Which means the function can work with a partial order which is completed with other
-      -- variables by default.
-      theOrder = e ++ r
-      addDVToBucket (rf, m) dv  =
-        let (fk,remaining) = partition (flip containsVariable dv) rf
-        in 
-        (remaining, M.insert dv fk m)
-      (_,b) = foldl' addDVToBucket (s,M.empty) theOrder
-  in
-  Buckets theOrder b
 
--- | Get the factors for a bucket
-getBucket :: DV 
-          -> Buckets f 
-          -> [f]
-getBucket dv (Buckets _ m) = fromJust $ M.lookup dv m
-
--- | Update bucket
-updateBucket :: Factor f
-             => DV -- ^ Variable that was eliminated
-             -> f -- ^ New factor resulting from this elimination
-             -> Buckets f 
-             -> Buckets f 
-updateBucket dv f b@(Buckets e m) = 
-  if isScalarFactor f 
-    then 
-      Buckets (remainingVarsToProcess e) (M.insert dv [f] m)
-    else
-      let b' = removeFromBucket dv b
-      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
-addBucket (Buckets e b) f = 
-  let inBucket = find (f `containsVariable`) e
-  in 
-  case inBucket of 
-    Nothing -> Buckets e b
-    Just bucket -> Buckets e (M.insertWith' (++) bucket [f] b)
-
--- | Remove a variable from the bucket
-removeFromBucket :: DV -> Buckets f -> Buckets f 
-removeFromBucket dv (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 :: (IsBucketItem f, Factor f)
          => [f] -- ^ Bayesian Network
          -> EliminationOrder DV -- ^ Ordering of variables to marginalize
          -> EliminationOrder DV -- ^ Ordering of remaining variables
@@ -184,7 +113,7 @@
     in 
     mpeInstantiations (resultFactor)
 
-posteriorMarginal :: (Graph g, Factor f, Show f, BayesianDiscreteVariable dva, BayesianDiscreteVariable dvb) 
+posteriorMarginal :: (Graph g, IsBucketItem f, Factor f,Show f, BayesianDiscreteVariable dva, BayesianDiscreteVariable dvb) 
                   => BayesianNetwork g f -- ^ Bayesian Network
                   -> EliminationOrder dva -- ^ Ordering of variables to marginzalie
                   -> EliminationOrder dvb-- ^ Ordering of remaining variables
@@ -202,7 +131,7 @@
 
 -- | Compute the prior marginal. All the variables in the
 -- elimination order are conditionning variables ( p( . | conditionning variables) )
-priorMarginal :: (Graph g, Factor f, Show f, BayesianDiscreteVariable dva, BayesianDiscreteVariable dvb) 
+priorMarginal :: (Graph g, IsBucketItem f, Factor f,Show f, BayesianDiscreteVariable dva, BayesianDiscreteVariable dvb) 
               => BayesianNetwork g f -- ^ Bayesian Network
               -> EliminationOrder dva-- ^ Ordering of variables to marginalize
               -> EliminationOrder dvb-- ^ Ordering of remaining to keep in result
diff --git a/Bayes/VariableElimination/Buckets.hs b/Bayes/VariableElimination/Buckets.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/VariableElimination/Buckets.hs
@@ -0,0 +1,111 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{- | Bucket algorithms for variable elimination with enough flexibility to
+also work with influence diagrams.
+    
+-}
+module Bayes.VariableElimination.Buckets(
+    -- * Types 
+      Buckets(..)
+    , EliminationOrder(..)
+    , IsBucketItem(..)
+    -- * Functions 
+    , createBuckets 
+    , getBucket 
+    , updateBucket 
+    , addBucket 
+    , removeFromBucket 
+    , marginalizeOneVariable
+    ) where 
+
+import Bayes.PrivateTypes
+import qualified Data.Map as M
+import Data.List(partition,minimumBy,(\\),find,foldl')
+import Data.Maybe(fromJust)
+
+-- | Elimination order
+type EliminationOrder dv = [dv]
+
+-- | Used for bucket elimination. Factor are organized by their first DV
+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"
+
+-- | Operations needed to process a bucket items
+class IsBucketItem f where 
+    scalarItem :: f -> Bool 
+    itemProduct :: [f] -> f
+    itemProjectOut :: DV -> f -> f
+    itemContainsVariable :: f -> DV  -> Bool
+
+
+
+addDVToBucket :: IsBucketItem f => ([f],M.Map DV [f]) -> DV -> ([f],M.Map DV [f]) 
+addDVToBucket (rf, m) dv  =
+  let (fk,remaining) = partition (flip itemContainsVariable dv) rf
+  in 
+  (remaining, M.insert dv fk m)
+
+createBuckets ::  (IsBucketItem f) 
+              => [f] -- ^ Factor to use for computing the marginal one
+              -> 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
+      -- Which means the function can work with a partial order which is completed with other
+      -- variables by default.
+      theOrder = e ++ r
+      (_,b) = foldl' addDVToBucket (s,M.empty) theOrder
+  in
+  Buckets theOrder b
+
+-- | Get the factors for a bucket
+getBucket :: DV 
+          -> Buckets f 
+          -> [f]
+getBucket dv (Buckets _ m) = fromJust $ M.lookup dv m
+
+-- | Update bucket
+updateBucket :: IsBucketItem f
+             => DV -- ^ Variable that was eliminated
+             -> f -- ^ New factor resulting from this elimination
+             -> Buckets f 
+             -> Buckets f 
+updateBucket dv f b@(Buckets e m) = 
+  if scalarItem f 
+    then 
+      Buckets (remainingVarsToProcess e) (M.insert dv [f] m)
+    else
+      let b' = removeFromBucket dv b
+      in
+      addBucket b' f 
+ where 
+  remainingVarsToProcess [] = []
+  remainingVarsToProcess l = tail l
+
+-- | Add a factor to the right bucket
+addBucket :: IsBucketItem f => Buckets f -> f -> Buckets f
+addBucket (Buckets e b) f = 
+  let inBucket = find (f `itemContainsVariable`) e
+  in 
+  case inBucket of 
+    Nothing -> Buckets e b
+    Just bucket -> Buckets e (M.insertWith' (++) bucket [f] b)
+
+-- | Remove a variable from the bucket
+removeFromBucket :: DV -> Buckets f -> Buckets f 
+removeFromBucket dv (Buckets [] m) = Buckets [] (M.delete dv m) 
+removeFromBucket dv (Buckets e m) = Buckets (tail e) (M.delete dv m) 
+
+marginalizeOneVariable :: IsBucketItem f => Buckets f -> DV -> Buckets f
+marginalizeOneVariable currentBucket dv   = 
+  let fk = getBucket dv currentBucket
+      p = itemProduct fk
+      f' = itemProjectOut dv p
+  in
+  updateBucket dv f' currentBucket
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.3
+Version:             0.4
 
 -- A short (one-line) description of the package.
 Synopsis:            Inference with Discrete Bayesian Networks
@@ -74,6 +74,12 @@
     Bayes.Examples
     Bayes.Examples.Tutorial
     Bayes.Test.ReferencePatterns
+    Bayes.Test.InfluencePatterns
+    Bayes.ImportExport
+    Bayes.BayesianNetwork
+    Bayes.InfluenceDiagram
+    Bayes.VariableElimination.Buckets
+    Bayes.Examples.Influence
   other-modules:
     Paths_hbayes
     Bayes.ImportExport.HuginNet.Splitting
@@ -81,6 +87,7 @@
     Bayes.FactorElimination.JTree
     Bayes.Tools
     Bayes.Factor.PrivateCPT
+    Bayes.Network
 
   GHC-Options: -funbox-strict-fields
   Extensions: CPP
@@ -105,6 +112,7 @@
     parsec,
     filepath,
     directory,
+    binary >= 0.5,
     test-framework-quickcheck2,
     test-framework,
     test-framework-hunit,