Merge branch 'eaton'
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 4 Aug 2009 15:35:49 +0000 (17:35 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 4 Aug 2009 15:35:49 +0000 (17:35 +0200)
Conflicts:

ChangeLog
Makefile
include/dai/daialg.h

19 files changed:
1  2 
ChangeLog
Makefile
Makefile.CYGWIN
Makefile.LINUX
Makefile.MACOSX
Makefile.WINDOWS
README
include/dai/alldai.h
include/dai/bipgraph.h
include/dai/bp.h
include/dai/daialg.h
include/dai/doc.h
include/dai/factorgraph.h
include/dai/lc.h
include/dai/mr.h
include/dai/util.h
src/bp.cpp
src/factorgraph.cpp
tests/testdai.cpp

diff --combined ChangeLog
+++ b/ChangeLog
@@@ -1,15 -1,25 +1,33 @@@
- Branch master:
 -git b6efbfd9e4e76de65c44e4060a9b56f99543b70a (head of branch eaton before merge)
 ---------------------------------------------------------------------------------
 +* Added work-around for bug in Boost Graph Library
 +* Improvements of TFactor<T>:
 +  - Extended functionality of TFactor<T>::operators +,-,+=,-= to handle different VarSets
 +  - Added TFactor<T>::maxMarginal (similar to marginal but with max instead of sum)
 +* Added BipartiteGraph::eraseEdge
 +* Removed all the virtual default constructors *::create(), as they are not used.
 +* Fixed bug in MaxSpanningTree (it wrongly assumed that the tree was not empty).
 +* [Charlie Vaske] Added Expectation Maximization code.
 +* Added MatLab QuickStart to README.
 +* MEX file dai now also returns variable and factor beliefs.
+ * Cleanups and extra documentation of the code contributed by Frederik Eaton
+ * Removed erroneous 'inline' directives in gibbs.cpp
+ * [Frederik Eaton] Added script for automatically generating properties code (used by CBP and BBP)
+ * [Frederik Eaton] Updated doxygen.conf to version 1.5.9
+ * [Frederik Eaton] Added newInfAlgFromString to alldai.h/cpp
+ * [Frederik Eaton] Added FactorGraph::clampVar and FactorGraph::clampFactor
+ * [Frederik Eaton] Changed BP to be able to record order in which messages are sent (for use by BBP)
+ * [Frederik Eaton] Improved documentation of BipartiteGraph::Neighbor
+ * [Frederik Eaton]: Extended InfAlg interface with beliefV() and beliefF() methods;
+ * [Frederik Eaton]: Improved PropertySet class:
+   - Added default constructor
+   - Allow construction from string
+   - Added method allKeys() to produce list of keys
+   - In getStringAs(), check for typeid(ValueType) before typeid(std::string)
+     (this fixes a strange bug for empty string property)
+ * [Frederik Eaton]: Added some utility macros (DAI_PV, DAI_DMSG, DAI_ACCMUT, DAI_IFVERB)
+   and functions (concat(),rnd()) to include/dai/util.h
+ * [Frederik Eaton] Added backwards compatibility layer for edge interface to
+   BipartiteGraph and FactorGraph (which will be obsoleted soon)
+ * [Frederik Eaton] Added BP_dual, BBP and CBP algorithms
  * [Frederik Eaton] Added Gibbs::state() accessors/mutators
  * [Frederik Eaton] Fixed Gibbs::getVarDist(size_t) to return uniform 
    distribution if no state is allowed
@@@ -47,7 -57,7 +65,7 @@@
    errors (thanks to Dan Preston for reporting this)
  * toc() now returns POSIX system time with maximum accuracy of microseconds
  * Exception objects now remember their error code
 -* Added examples/example_springler.cpp
 +* Added examples/example_sprinkler.cpp
  
  
  git 065eae35cbfcc36f1a945ae3053c80c23f366306
diff --combined Makefile
+++ b/Makefile
@@@ -31,7 -31,7 +31,7 @@@ SRC=sr
  LIB=lib
  
  # Define build targets
 -TARGETS=tests utils lib examples testregression
 +TARGETS=tests utils lib examples testregression testem
  ifdef WITH_DOC
    TARGETS:=$(TARGETS) doc 
  endif
@@@ -45,7 -45,7 +45,7 @@@ ifdef WITH_MATLA
  endif
  
  # Define conditional build targets
 -OBJECTS:=exactinf$(OE)
 +OBJECTS:=exactinf$(OE) evidence$(OE) emalg$(OE)
  ifdef WITH_BP
    CCFLAGS:=$(CCFLAGS) -DDAI_WITH_BP
    OBJECTS:=$(OBJECTS) bp$(OE)
@@@ -78,9 -78,14 +78,14 @@@ ifdef WITH_GIBB
    CCFLAGS:=$(CCFLAGS) -DDAI_WITH_GIBBS
    OBJECTS:=$(OBJECTS) gibbs$(OE)
  endif
+ ifdef WITH_CBP
+   CCFLAGS:=$(CCFLAGS) -DDAI_WITH_CBP
+   OBJECTS:=$(OBJECTS) bbp$(OE) cbp$(OE) bp_dual$(OE)
+ endif
  
  # Define standard libDAI header dependencies
- HEADERS=$(INC)/bipgraph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h
+ HEADERS=$(INC)/bipgraph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h $(INC)/util.h
  
  # Setup final command for C++ compiler and MEX
  ifdef DEBUG
@@@ -104,7 -109,7 +109,7 @@@ examples : examples/example$(EE) exampl
  
  matlabs : matlab/dai$(ME) matlab/dai_readfg$(ME) matlab/dai_writefg$(ME) matlab/dai_potstrength$(ME)
  
- tests : tests/testdai$(EE) tests/testem/testem$(EE)
 -tests : tests/testdai$(EE) tests/testbbp$(EE)
++tests : tests/testdai$(EE) tests/testem/testem$(EE) tests/testbbp$(EE)
  
  utils : utils/createfg$(EE) utils/fg2dot$(EE) utils/fginfo$(EE)
  
@@@ -126,6 -131,15 +131,15 @@@ exactinf$(OE) : $(SRC)/exactinf.cpp $(I
  bp$(OE) : $(SRC)/bp.cpp $(INC)/bp.h $(HEADERS)
        $(CC) -c $(SRC)/bp.cpp
  
+ bp_dual$(OE) : $(SRC)/bp_dual.cpp $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/bp_dual.cpp
+ bbp$(OE) : $(SRC)/bbp.cpp $(INC)/bbp.h $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/bbp.cpp
+ cbp$(OE) : $(SRC)/cbp.cpp $(INC)/cbp.h $(INC)/bbp.h $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/cbp.cpp
  lc$(OE) : $(SRC)/lc.cpp $(INC)/lc.h $(HEADERS)
        $(CC) -c $(SRC)/lc.cpp
  
@@@ -162,12 -176,6 +176,12 @@@ mr$(OE) : $(SRC)/mr.cpp $(INC)/mr.h $(H
  gibbs$(OE) : $(SRC)/gibbs.cpp $(INC)/gibbs.h $(HEADERS)
        $(CC) -c $(SRC)/gibbs.cpp
  
 +evidence$(OE) : $(SRC)/evidence.cpp $(INC)/evidence.h $(HEADERS)
 +      $(CC) -c $(SRC)/evidence.cpp
 +
 +emalg$(OE) : $(SRC)/emalg.cpp $(INC)/emalg.h $(INC)/evidence.h $(HEADERS)
 +      $(CC) -c $(SRC)/emalg.cpp
 +
  properties$(OE) : $(SRC)/properties.cpp $(HEADERS)
        $(CC) -c $(SRC)/properties.cpp
  
@@@ -199,9 -207,10 +213,12 @@@ examples/example_sprinkler$(EE) : examp
  
  tests/testdai$(EE) : tests/testdai.cpp $(HEADERS) $(LIB)/libdai$(LE)
        $(CC) $(CCO)tests/testdai$(EE) tests/testdai.cpp $(LIBS) $(BOOSTLIBS)
 +tests/testem/testem$(EE) : tests/testem/testem.cpp $(HEADERS) $(LIB)/libdai$(LE)
 +      $(CC) $(CCO)$@ $< $(LIBS) $(BOOSTLIBS)
  
+ tests/testbbp$(EE) : tests/testbbp.cpp $(HEADERS) $(LIB)/libdai$(LE)
+       $(CC) $(CCO)tests/testbbp$(EE) tests/testbbp.cpp $(LIBS)
  
  # MATLAB INTERFACE
  ###################
@@@ -252,22 -261,16 +269,22 @@@ endi
  # REGRESSION TESTS
  ###################
  
 -ifneq ($(OS),WINDOWS)
  testregression : tests/testdai$(EE)
        @echo Starting regression test...this can take a minute or so!
 +ifneq ($(OS),WINDOWS)
        cd tests && ./testregression && cd ..
  else
 -testregression : tests/testdai$(EE)
 -      @echo Starting regression test...this can take a minute or so!
        cd tests && testregression.bat && cd ..
  endif
  
 +testem : tests/testem/testem$(EE)
 +      @echo Starting EM tests
 +ifneq ($(OS),WINDOWS)
 +      cd tests/testem && ./runtests && cd ../..
 +else
 +      cd tests\testem && runtests && cd ..\..
 +endif
 +
  
  # DOCUMENTATION
  ################
@@@ -289,13 -292,12 +306,12 @@@ clean 
        -rm *$(OE)
        -rm matlab/*$(ME)
        -rm examples/example$(EE) examples/example_bipgraph$(EE) examples/example_varset$(EE) examples/example_sprinkler$(EE)
-       -rm tests/testdai$(EE)
-       -rm tests/testem/testem$(EE)
 -      -rm tests/testdai$(EE) tests/testbbp$(EE)
++      -rm tests/testdai$(EE) tests/testem/testem$(EE) tests/testbbp$(EE)
        -rm utils/fg2dot$(EE) utils/createfg$(EE) utils/fginfo$(EE)
        -rm -R doc
        -rm -R lib
  else
  .PHONY : clean
  clean :
 -      -del *$(OE) *.ilk *.pdb *$(EE) matlab\*$(ME) examples\*$(EE) examples\*.ilk examples\*.pdb tests\testdai$(EE) tests\*.pdb tests\*.ilk utils\*$(EE) utils\*.pdb utils\*.ilk $(LIB)\libdai$(LE)
 +      -del *$(OE) *.ilk *.pdb *$(EE) matlab\*$(ME) examples\*$(EE) examples\*.ilk examples\*.pdb tests\testdai$(EE) tests\testem\testem$(EE) tests\*.pdb tests\*.ilk utils\*$(EE) utils\*.pdb utils\*.ilk $(LIB)\libdai$(LE)
  endif
diff --combined Makefile.CYGWIN
@@@ -15,13 -15,14 +15,14 @@@ WITH_TREEEP=tru
  WITH_JTREE=true
  WITH_MR=true
  WITH_GIBBS=true
+ WITH_CBP=true
  # Build doxygen documentation?
  WITH_DOC=true
  # Build with debug info?
  DEBUG=true
  # Build matlab interface?
  WITH_MATLAB=
 -# New/old matlab version?
 +# MatLab version 7.3 (R2006b) or newer?
  NEW_MATLAB=true
  
  # OPERATING SYSTEM
diff --combined Makefile.LINUX
@@@ -16,13 -16,14 +16,14 @@@ WITH_TREEEP=tru
  WITH_JTREE=true
  WITH_MR=true
  WITH_GIBBS=true
+ WITH_CBP=true
  # Build doxygen documentation?
  WITH_DOC=true
  # Build with debug info?
  DEBUG=true
  # Build matlab interface?
  WITH_MATLAB=
 -# New/old matlab version?
 +# MatLab version 7.3 (R2006b) or newer?
  NEW_MATLAB=true
  
  # OPERATING SYSTEM
diff --combined Makefile.MACOSX
@@@ -15,13 -15,14 +15,14 @@@ WITH_TREEEP=tru
  WITH_JTREE=true
  WITH_MR=true
  WITH_GIBBS=true
+ WITH_CBP=true
  # Build doxygen documentation?
  WITH_DOC=true
  # Build with debug info?
  DEBUG=true
  # Build matlab interface?
  WITH_MATLAB=
 -# New/old matlab version?
 +# MatLab version 7.3 (R2006b) or newer?
  NEW_MATLAB=true
  
  # OPERATING SYSTEM
diff --combined Makefile.WINDOWS
@@@ -16,13 -16,14 +16,14 @@@ WITH_TREEEP=tru
  WITH_JTREE=true
  WITH_MR=true
  WITH_GIBBS=true
+ WITH_CBP=true
  # Build doxygen documentation?
  WITH_DOC=
  # Build with debug info?
  DEBUG=true
  # Build matlab interface?
  WITH_MATLAB=
 -# New/old matlab version?
 +# MatLab version 7.3 (R2006b) or newer?
  NEW_MATLAB=true
  
  # OPERATING SYSTEM
diff --combined README
--- 1/README
--- 2/README
+++ b/README
@@@ -13,7 -13,6 +13,7 @@@ with contributions from
  Martijn Leisink
  Giuseppe Passino
  Frederik Eaton
 +Charlie Vaske
  Bastian Wemmenhove
  Christian Wojek
  Claudio Lima
@@@ -48,7 -47,7 +48,7 @@@ to the literature, to allow replication
  Acknowledgements section.  The appropriate citation is: 
  
  J. M. Mooij (2008) "libDAI 0.2.2: A free/open source C++ library for Discrete 
 -Approximate Inference methods", http://mloss.org/software/view/77/.
 +Approximate Inference methods", http://www.libdai.org
  
  Moreover, as a personal note, I would appreciate it if you would email me with
  citations of papers referencing this work so I can mention them to my funding
@@@ -88,22 -87,20 +88,23 @@@ Currently, libDAI supports the followin
      * Generalized Belief Propagation [YFW05];
      * Double-loop GBP [HAK03];
      * Various variants of Loop Corrected Belief Propagation [MoK07, MoR05];
-     * Gibbs sampler.
+     * Gibbs sampler;
+     * Conditioned BP [EaG09].
  
 +In addition, libDAI supports parameter learning of conditional probability
 +tables by Expectation Maximization.
 +
  
  Why C++?
  --------
  Because libDAI is implemented in C++, it is very fast compared with
  implementations in MatLab (a factor 1000 faster is not uncommon). libDAI does
 -provide a MatLab interface for easy integration with MatLab.
 +provide a (limited) MatLab interface for easy integration with MatLab.
  
  
  Releases
  --------
 -Releases can be obtained from http://mloss.org/software/view/77/
 +Releases can be obtained from www.libdai.org
  License: GNU Public License v2 (or higher).
  
  libDAI-0.2      December 1, 2006
@@@ -168,7 -165,7 +169,7 @@@ from http://www.boost.org/ and compile/
  To build the libDAI source, first copy a template Makefile.* to Makefile.conf
  (for example, copy Makefile.LINUX to Makefile.conf if you use GNU/Linux). 
  Then, edit the Makefile.conf template to adapt it to your local setup.
 -Especially directories may change from system to system. Finally, run
 +Especially directories may differ from system to system. Finally, run
      
      make
  
@@@ -202,42 -199,3 +203,42 @@@ If the build was successful, you can te
  or the more elaborate test program:
  
      tests\testdai --aliases tests\aliases.conf --filename tests\alarm.fg --methods JTREE_HUGIN BP_SEQMAX
 +
 +
 +Quick start (MatLab)
 +--------------------
 +You need:
 +- MatLab
 +- The platform-dependent requirements described above
 +
 +First, you need to build the libDAI source as described above for your
 +platform. By default, the MatLab interface is disabled, so before compiling the
 +source, you have to enable it in the Makefile.conf by setting
 +"WITH_MATLAB=true". Also, you have to configure the MatLab-specific parts of
 +Makefile.conf to match your system (in particular, the Makefile variables ME,
 +MATLABDIR and MEX). The MEX file extension depends on your platform; for a
 +64-bit linux x86_64 system this would be "ME=.mexa64", for a 32-bit linux x86
 +system this would be "ME=.mexglx". If you are unsure about your MEX file
 +extension: it needs to be the same as what the MatLab command "mexext" returns.
 +The required MEX files are built by issuing
 +
 +    make
 +
 +from the command line. The MatLab interface is much less powerful than using
 +libDAI from C++. There are two reasons for this: (i) it is boring to write MEX
 +files; (ii) the large performance penalty paid when large data structures (like
 +factor graphs) have to be converted between their native C++ data structure to
 +something that MatLab understands.
 +
 +A simple example of how to use the MatLab interface is the following (entered
 +at the MatLab prompt), which performs exact inference by the junction tree
 +algorithm and approximate inference by belief propagation on the ALARM network:
 +
 +    cd path_to_libdai/matlab
 +    [psi] = dai_readfg ('../examples/alarm.fg');
 +    [logZ,q,md,qv,qf] = dai (psi, 'JTREE', '[updates=HUGIN,verbose=0]')
 +    [logZ,q,md,qv,qf] = dai (psi, 'BP', '[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0]')
 +
 +where "path_to_libdai" has to be replaced with the directory in which libDAI
 +was installed. For other algorithms and their parameters, see
 +tests/aliases.conf.
diff --combined include/dai/alldai.h
@@@ -33,8 -33,6 +33,8 @@@
  #include <dai/daialg.h>
  #include <dai/properties.h>
  #include <dai/exactinf.h>
 +#include <dai/evidence.h>
 +#include <dai/emalg.h>
  #ifdef DAI_WITH_BP
      #include <dai/bp.h>
  #endif
@@@ -59,6 -57,9 +59,9 @@@
  #ifdef DAI_WITH_GIBBS
      #include <dai/gibbs.h>
  #endif
+ #ifdef DAI_WITH_CBP
+     #include <dai/cbp.h>
+ #endif
  
  
  /// Namespace for libDAI
@@@ -74,6 -75,14 +77,14 @@@ namespace dai 
  InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const PropertySet &opts );
  
  
+ /// Constructs a new approximate inference algorithm.
+ /** \param nameOpts The name and options of the approximate inference algorithm (should be in the format "name[opts]").
+  *  \param fg The FactorGraph that the algorithm should be applied to.
+  *  \return Returns a pointer to the new InfAlg object; it is the responsibility of the caller to delete it later.
+  */
+ InfAlg *newInfAlgFromString( const std::string &nameOpts, const FactorGraph &fg );
  /// Contains the names of all approximate inference algorithms compiled into libDAI.
  static const char* DAINames[] = {
      ExactInf::Name,
  #endif
  #ifdef DAI_WITH_GIBBS
      Gibbs::Name,
+ #endif
+ #ifdef DAI_WITH_CBP
+     CBP::Name,
  #endif
      ""
  };
diff --combined include/dai/bipgraph.h
@@@ -56,14 -56,59 +56,59 @@@ namespace dai 
   */
  class BipartiteGraph {
      public:
-         /// Describes a neighboring node of some other node in a BipartiteGraph.
-         /** A Neighbor structure has three members: \a iter, \a node and \a dual. The \a 
-          *  node member is the most important member: it contains the index of the neighboring node. The \a iter 
-          *  member is useful for iterating over neighbors, and contains the index of this Neighbor entry in the
-          *  corresponding BipartiteGraph::Neighbors variable. The \a dual member is useful to find the dual Neighbor 
-          *  element: a pair of neighboring nodes can be either specified as a node of type 1 and a neighbor of type 
-          *  2, or as a node of type 2 and a neighbor of type 1; the \a dual member contains the index of the dual 
-          *  Neighbor element (see the example for another explanation of the dual member).
+         /// Describes the neighbor relationship of two nodes in a BipartiteGraph.
+         /** Sometimes we want to do an action, such as sending a
+          *  message, for all edges in a graph. However, most graphs
+          *  will be sparse, so we need some way of storing a set of
+          *  the neighbors of a node, which is both fast and
+          *  memory-efficient. We also need to be able to go between
+          *  viewing node \a A as a neighbor of node \a B, and node \a
+          *  B as a neighbor of node \a A. The Neighbor struct solves
+          *  both of these problems. Each node has a list of neighbors,
+          *  stored as a vector<Neighbor>, and extra information is
+          *  included in the Neighbor struct which allows us to access
+          *  a node as a neighbor of its neighbor (the \a dual member).
+          *
+          *  By convention, variable identifiers naming indices into a
+          *  vector of neighbors are prefixed with an underscore ("_"). 
+          *  The neighbor list which they point into is then understood
+          *  from context. For example:
+          *
+          *  \code
+          *  void BP::calcNewMessage( size_t i, size_t _I )
+          *  \endcode
+          *
+          *  Here, \a i is the "absolute" index of node i, but \a _I is
+          *  understood as a "relative" index, giving node I's entry in
+          *  nb1(i). The corresponding Neighbor structure can be
+          *  accessed as nb1(i,_I) or nb1(i)[_I]. The absolute index of
+          *  \a _I, which would be called \a I, can be recovered from
+          *  the \a node member: nb1(i,_I).node. The \a iter member
+          *  gives the relative index \a _I, and the \a dual member
+          *  gives the "dual" relative index, i.e. the index of \a i in
+          *  \a I's neighbor list.
+          *
+          *  \code
+          *  Neighbor n = nb1(i,_I);
+          *  n.node == I &&
+          *  n.iter == _I &&
+          *  nb2(n.node,n.dual).node == i
+          *  \endcode
+          *
+          *  In a FactorGraph, nb1 is called nbV, and nb2 is called
+          *  nbF.
+          *
+          *  There is no easy way to transform a pair of absolute node
+          *  indices \a i and \a I into a Neighbor structure relative
+          *  to one of the nodes. Such a feature has never yet been
+          *  found to be necessary. Iteration over edges can always be
+          *  accomplished using the Neighbor lists, and by writing
+          *  functions that accept relative indices:
+          *  \code
+          *  for( size_t i = 0; i < nrVars(); ++i )
+          *      foreach( const Neighbor &I, nbV(i) )
+          *          calcNewMessage( i, I.iter );
+          *  \endcode
           */
          struct Neighbor {
              /// Corresponds to the index of this Neighbor entry in the vector of neighbors
              std::vector<size_t> ind2;       // indices of nodes of type 2
          };
  
+         /// @name Backwards compatibility layer (to be removed soon)
+         //@{
+         /// Enable backwards compatibility layer?
+         bool _edge_indexed;
+         /// Call indexEdges() first to initialize these members
+         std::vector<Edge> _edges;
+         /// Call indexEdges() first to initialize these members
+         hash_map<Edge,size_t> _vv2e;
+         //}@
      public:
          /// Default constructor (creates an empty bipartite graph)
-         BipartiteGraph() : _nb1(), _nb2() {}
+         BipartiteGraph() : _nb1(), _nb2(), _edge_indexed(false) {}
  
          /// Constructs BipartiteGraph from a range of edges.
          /** \tparam EdgeInputIterator Iterator that iterates over instances of BipartiteGraph::Edge.
           *  \param end Points just beyond the last edge.
           */
          template<typename EdgeInputIterator>
-         BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ) {
+         BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ), _edge_indexed(false) {
              construct( nr1, nr2, begin, end );
          }
  
          /// Removes node n2 of type 2 and all incident edges.
          void erase2( size_t n2 );
  
 +        /// Removes edge between node n1 of type 1 and node n2 of type 2.
 +        void eraseEdge( size_t n1, size_t n2 ) {
 +            assert( n1 < nr1() );
 +            assert( n2 < nr2() );
 +            for( Neighbors::iterator i1 = _nb1[n1].begin(); i1 != _nb1[n1].end(); i1++ )
 +                if( i1->node == n2 ) {
 +                    _nb1[n1].erase( i1 );
 +                    break;
 +                }
 +            for( Neighbors::iterator i2 = _nb2[n2].begin(); i2 != _nb2[n2].end(); i2++ )
 +                if( i2->node == n1 ) {
 +                    _nb2[n2].erase( i2 );
 +                    break;
 +                }
 +        }
 +
          /// Adds an edge between node n1 of type 1 and node n2 of type 2.
          /** If check == true, only adds the edge if it does not exist already.
           */
          /// Writes this BipartiteGraph to an output stream in GraphViz .dot syntax
          void printDot( std::ostream& os ) const;
  
+         /// @name Backwards compatibility layer (to be removed soon)
+         //@{
+         void indexEdges() {
+             std::cerr << "Warning: this BipartiteGraph edge interface is obsolete!" << std::endl;
+             _edges.clear();
+             _vv2e.clear();
+             size_t i=0;
+             foreach(const Neighbors &nb1s, _nb1) {
+                 foreach(const Neighbor &n2, nb1s) {
+                     Edge e(i, n2.node);
+                     _edges.push_back(e);
+                 }
+                 i++;
+             }
+             sort(_edges.begin(), _edges.end()); // unnecessary?
+             i=0;
+             foreach(const Edge& e, _edges) {
+                 _vv2e[e] = i++;
+             }
+             _edge_indexed = true;
+         }
+         
+         const Edge& edge(size_t e) const {
+             assert(_edge_indexed);
+             return _edges[e];
+         }
+         
+         const std::vector<Edge>& edges() const {
+             return _edges;
+         }
+         
+         size_t VV2E(size_t n1, size_t n2) const {
+             assert(_edge_indexed);
+             Edge e(n1,n2);
+             hash_map<Edge,size_t>::const_iterator i = _vv2e.find(e);
+             assert(i != _vv2e.end());
+             return i->second;
+         }
+         size_t nr_edges() const {
+             assert(_edge_indexed);
+             return _edges.size();
+         }
+         //}@
      private:
          /// Checks internal consistency
          void check() const;
diff --combined include/dai/bp.h
@@@ -57,15 -57,17 +57,17 @@@ class BP : public DAIAlgFG 
          double _maxdiff;
          /// Number of iterations needed
          size_t _iters;
+         /// The history of message updates (only recorded if recordSentMessages is true)
+         std::vector<std::pair<std::size_t, std::size_t> > _sentMessages;
      
      public:
          /// Parameters of this inference algorithm
          struct Properties {
              /// Enumeration of possible update schedules
-             DAI_ENUM(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL)
+             DAI_ENUM(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL);
  
              /// Enumeration of inference variants
-             DAI_ENUM(InfType,SUMPROD,MAXPROD)
+             DAI_ENUM(InfType,SUMPROD,MAXPROD);
          
              /// Verbosity
              size_t verbose;
          /// Name of this inference algorithm
          static const char *Name;
  
+         /// Specifies whether the history of message updates should be recorded
+         bool recordSentMessages;
      public:
          /// Default constructor
-         BP() : DAIAlgFG(), _edges(), _edge2lut(), _lut(), _maxdiff(0.0), _iters(0U), props() {}
+         BP() : DAIAlgFG(), _edges(), _edge2lut(), _lut(), _maxdiff(0.0), _iters(0U), _sentMessages(), props(), recordSentMessages(false) {}
  
          /// Copy constructor
-         BP( const BP &x ) : DAIAlgFG(x), _edges(x._edges), _edge2lut(x._edge2lut), _lut(x._lut), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {
+         BP( const BP &x ) : DAIAlgFG(x), _edges(x._edges), _edge2lut(x._edge2lut),
+             _lut(x._lut), _maxdiff(x._maxdiff), _iters(x._iters), _sentMessages(x._sentMessages), 
+             props(x.props), recordSentMessages(x.recordSentMessages) 
+         {
              for( LutType::iterator l = _lut.begin(); l != _lut.end(); ++l )
                  _edge2lut[l->second.first][l->second.second] = l;
          }
                      _edge2lut[l->second.first][l->second.second] = l;
                  _maxdiff = x._maxdiff;
                  _iters = x._iters;
+                 _sentMessages = x._sentMessages;
                  props = x.props;
+                 recordSentMessages = x.recordSentMessages;
              }
              return *this;
          }
  
          /// Construct from FactorGraph fg and PropertySet opts
-         BP( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _edges(), _maxdiff(0.0), _iters(0U), props() {
+         BP( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _edges(), _maxdiff(0.0), _iters(0U), _sentMessages(), props(), recordSentMessages(false) {
              setProperties( opts );
              construct();
          }
          /// @name General InfAlg interface
          //@{
          virtual BP* clone() const { return new BP(*this); }
 -        virtual BP* create() const { return new BP(); }
          virtual std::string identify() const;
          virtual Factor belief( const Var &n ) const;
          virtual Factor belief( const VarSet &ns ) const;
           */
          std::vector<std::size_t> findMaximum() const;
  
+         /// Returns history of sent messages
+         const std::vector<std::pair<std::size_t, std::size_t> >& getSentMessages() const {
+             return _sentMessages;
+         }
+         /// Clears history of sent messages
+         void clearSentMessages() {
+             _sentMessages.clear();
+         }
      private:
          const Prob & message(size_t i, size_t _I) const { return _edges[i][_I].message; }
          Prob & message(size_t i, size_t _I) { return _edges[i][_I].message; }
diff --combined include/dai/daialg.h
@@@ -54,6 -54,9 +54,6 @@@ class InfAlg 
          /// Returns a pointer to a new, cloned copy of *this (i.e., virtual copy constructor)
          virtual InfAlg* clone() const = 0;
  
 -        /// Returns a pointer to a newly constructed object *this (i.e., virtual default constructor)
 -        virtual InfAlg* create() const { DAI_THROW(NOT_IMPLEMENTED); }
 -        
          /// Identifies itself for logging purposes
          virtual std::string identify() const = 0;
  
          /// Returns the "belief" (i.e., approximate marginal probability distribution) of a set of variables
          virtual Factor belief( const VarSet &n ) const = 0;
  
+         /// Returns marginal for a variable.
+         /** Sometimes preferred to belief() for performance reasons.
+           * Faster implementations exist in e.g. BP.
+           */
+         virtual Factor beliefV( size_t i ) const { return belief( fg().var(i) ); }
+         /// Returns marginal for a factor.
+         /** Sometimes preferred to belief() for performance reasons.
+           * Faster implementations exist in e.g. BP.
+           */
+         virtual Factor beliefF( size_t I ) const { return belief( fg().factor(I).vars() ); }
          /// Returns all "beliefs" (i.e., approximate marginal probability distribution) calculated by the algorithm
          virtual std::vector<Factor> beliefs() const = 0;
  
diff --combined include/dai/doc.h
   *  - Double-loop GBP [\ref HAK03];
   *  - Various variants of Loop Corrected Belief Propagation
   *    [\ref MoK07, \ref MoR05];
-  *  - Gibbs sampler.
+  *  - Gibbs sampler;
+  *  - Conditioned BP [\ref EaG09];
   *
 + *  In addition, libDAI supports parameter learning of conditional probability
 + *  tables by Expectation Maximization.
 + *
   *  \section language Why C++?
   *  Because libDAI is implemented in C++, it is very fast compared with
   *  implementations in MatLab (a factor 1000 faster is not uncommon).
   */
  
  
 -/** \page biboliography Bibliography
 +/** \page bibliography Bibliography
   *  \section Bibliograpy
   *  \anchor KFL01 \ref KFL01
   *  F. R. Kschischang and B. J. Frey and H.-A. Loeliger (2001):
   *  "Sufficient Conditions for Convergence of the Sum-Product Algorithm",
   *  <em>IEEE Transactions on Information Theory</em> 53(12):4422-4437.
   *  http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4385778
+  *
+  *  \anchor EaG09 \ref EaG09
+  *  F. Eaton and Z. Ghahramani (2009):
+  *  "Choosing a Variable to Clamp",
+  *  <em>Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics (AISTATS 2009)</em> 5:145-152
+  *  http://jmlr.csail.mit.edu/proceedings/papers/v5/eaton09a/eaton09a.pdf
   */
  
  
@@@ -113,6 -113,9 +113,6 @@@ class FactorGraph 
          /// Clone *this (virtual copy constructor)
          virtual FactorGraph* clone() const { return new FactorGraph(); }
  
 -        /// Create (virtual default constructor)
 -        virtual FactorGraph* create() const { return new FactorGraph(*this); }
 -
          /// Returns const reference to i'th variable
          const Var & var(size_t i) const { return _vars[i]; }
          /// Returns const reference to all factors
              }
          }
  
-         /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$);
-         /// If backup == true, make a backup of all factors that are changed
+         /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$)
+         /** If backup == true, make a backup of all factors that are changed
+          */
          virtual void clamp( const Var & n, size_t i, bool backup = false );
  
+         /// Clamp a variable in a factor graph to have one out of a list of values
+         /** If backup == true, make a backup of all factors that are changed
+          */
+         void clampVar( size_t i, const std::vector<size_t> &xis, bool backup = false );
+         /// Clamp a factor in a factor graph to have one out of a list of values
+         /** If backup == true, make a backup of all factors that are changed
+          */
+         void clampFactor( size_t I, const std::vector<size_t> &xIs, bool backup = false );
          /// Set all factors interacting with the i'th variable 1
          virtual void makeCavity( unsigned i, bool backup = false );
  
          void ReadFromFile(const char *filename);
  
          /// Writes a FactorGraph to a file
 -        void WriteToFile(const char *filename) const;
 +        void WriteToFile(const char *filename, size_t precision=15) const;
  
          /// Writes a FactorGraph to a GraphViz .dot file
          void printDot( std::ostream& os ) const;
          friend std::ostream& operator << (std::ostream& os, const FactorGraph& fg);
          friend std::istream& operator >> (std::istream& is, FactorGraph& fg);
  
+         /// @name Backwards compatibility layer (to be removed soon)
+         //@{
+         size_t VV2E(size_t n1, size_t n2) const { return G.VV2E(n1,n2); }
+         const Edge& edge(size_t e) const { return G.edge(e); }
+         void indexEdges() { G.indexEdges(); }
+         size_t nr_edges() const { return G.nr_edges(); }
+         const std::vector<Edge>& edges() const { return G.edges(); }
+         //}@
      private:
          /// Part of constructors (creates edges, neighbors and adjacency matrix)
          void constructGraph( size_t nrEdges );
diff --combined include/dai/lc.h
@@@ -60,10 -60,10 +60,10 @@@ class LC : public DAIAlgFG 
          /// Parameters of this inference algorithm
          struct Properties {
              /// Enumeration of possible ways to initialize the cavities
-             DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM)
+             DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM);
  
              /// Enumeration of different update schedules
-             DAI_ENUM(UpdateType,SEQFIX,SEQRND,NONE)
+             DAI_ENUM(UpdateType,SEQFIX,SEQRND,NONE);
  
              /// Verbosity
              size_t verbose;
          /// @name General InfAlg interface
          //@{
          virtual LC* clone() const { return new LC(*this); }
 -        virtual LC* create() const { return new LC(); }
          virtual std::string identify() const;
          virtual Factor belief( const Var &n ) const { return( _beliefs[findVar(n)] ); }
          virtual Factor belief( const VarSet &/*ns*/ ) const { DAI_THROW(NOT_IMPLEMENTED); return Factor(); }
diff --combined include/dai/mr.h
@@@ -72,10 -72,10 +72,10 @@@ class MR : public DAIAlgFG 
          /// Parameters of this inference algorithm
          struct Properties {
              /// Enumeration of different types of update equations
-             DAI_ENUM(UpdateType,FULL,LINEAR)
+             DAI_ENUM(UpdateType,FULL,LINEAR);
  
              /// Enumeration of different ways of initializing the cavity correlations
-             DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT)
+             DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT);
  
              /// Verbosity
              size_t verbose;
          /// @name General InfAlg interface
          //@{
          virtual MR* clone() const { return new MR(*this); }
 -        virtual MR* create() const { return new MR(); }
          virtual std::string identify() const;
          virtual Factor belief( const Var &n ) const;
          virtual Factor belief( const VarSet &/*ns*/ ) const { DAI_THROW(NOT_IMPLEMENTED); return Factor(); }
diff --combined include/dai/util.h
@@@ -35,6 -35,7 +35,7 @@@
  #include <iostream>
  #include <cstdio>
  #include <boost/foreach.hpp>
+ #include <boost/functional/hash.hpp>
  #include <algorithm>
  
  
  /// An alias to the BOOST_FOREACH macro from the boost::foreach library
  #define foreach BOOST_FOREACH
  
+ #ifdef DAI_DEBUG
+ /// \brief "Print variable". Prints the text of an expression, followed by its value (only if DAI_DEBUG is defined)
+ /**
+  *  Useful debugging macro to see what your code is doing. 
+  *  Example: \code DAI_PV(3+4) \endcode
+  *  Output: \code 3+4= 7 \endcode
+  */
+ #define DAI_PV(x) do {std::cerr << #x "= " << (x) << std::endl;} while(0)
+ /// "Debugging message": Prints a message (only if DAI_DEBUG is defined)
+ #define DAI_DMSG(str) do {std::cerr << str << std::endl;} while(0)
+ #else
+ #define DAI_PV(x) do {} while(0)
+ #define DAI_DMSG(str) do {} while(0)
+ #endif
+ /// Produces accessor and mutator methods according to the common pattern.
+ /** Example:
+  *  \code DAI_ACCMUT(size_t& maxIter(), { return props.maxiter; }); \endcode
+  *  \todo At the moment, only the mutator appears in doxygen documentation.
+  */
+ #define DAI_ACCMUT(x,y)                     \
+       x y;                                  \
+       const x const y;
+ /// Macro to give error message \a stmt if props.verbose>=\a n
+ #define DAI_IFVERB(n, stmt) if(props.verbose>=n) { cerr << stmt; }
  
  /// Real number (alias for double, which could be changed to long double if necessary)
  typedef double Real;
@@@ -77,14 -105,14 +105,14 @@@ namespace dai 
      /// hash_map is an alias for std::map.
      /** Since there is no TR1 unordered_map implementation available yet, we fall back on std::map.
       */
-     template <typename T, typename U>
+     template <typename T, typename U, typename H = boost::hash<T> >
          class hash_map : public std::map<T,U> {};
  #else
      /// hash_map is an alias for std::tr1::unordered_map.
      /** We use the (experimental) TR1 unordered_map implementation included in modern GCC distributions.
       */
-     template <typename T, typename U>
-         class hash_map : public std::tr1::unordered_map<T,U> {};
+     template <typename T, typename U, typename H = boost::hash<T> >
+         class hash_map : public std::tr1::unordered_map<T,U,H> {};
  #endif
  
  
@@@ -104,6 -132,11 +132,11 @@@ double rnd_stdnormal()
  /// Returns a random integer in interval [min, max]
  int rnd_int( int min, int max );
  
+ /// Returns a random integer in the half-open interval \f$[0,n)\f$
+ inline int rnd( int n) {
+     return rnd_int( 0, n-1 );
+ }
  
  /// Writes a std::vector to a std::ostream
  template<class T> 
@@@ -142,6 -175,17 +175,17 @@@ std::ostream& operator << (std::ostream
      return os;
  }
  
+ /// Concatenate two vectors
+ template<class T>
+ std::vector<T> concat( const std::vector<T>& u, const std::vector<T>& v ) {
+     std::vector<T> w;
+     w.reserve( u.size() + v.size() );
+     for( size_t i = 0; i < u.size(); i++ )
+         w.push_back( u[i] );
+     for( size_t i = 0; i < v.size(); i++ )
+         w.push_back( v[i] );
+     return w;
+ }
  
  /// Used to keep track of the progress made by iterative algorithms
  class Diffs : public std::vector<double> {
  };
  
  
 +/// Split a string into tokens
 +void tokenizeString( const std::string& s,
 +                  std::vector<std::string>& outTokens,
 +                  const std::string& delim="\t\n" );
 +
 +
  } // end of namespace dai
  
  
diff --combined src/bp.cpp
@@@ -145,6 -145,18 +145,6 @@@ void BP::init() 
  
  
  void BP::findMaxResidual( size_t &i, size_t &_I ) {
 -/*
 -    i = 0;
 -    _I = 0;
 -    double maxres = residual( i, _I );
 -    for( size_t j = 0; j < nrVars(); ++j )
 -        foreach( const Neighbor &I, nbV(j) )
 -            if( residual( j, I.iter ) > maxres ) {
 -                i = j;
 -                _I = I.iter;
 -                maxres = residual( i, _I );
 -            }
 -*/
      assert( !_lut.empty() );
      LutType::const_iterator largestEl = _lut.end();
      --largestEl;
@@@ -157,36 -169,41 +157,36 @@@ void BP::calcNewMessage( size_t i, size
      // calculate updated message I->i
      size_t I = nbV(i,_I);
  
 -    if( !DAI_BP_FAST ) {
 -        /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
 -        Factor prod( factor( I ) );
 -        foreach( const Neighbor &j, nbF(I) )
 -            if( j != i ) {     // for all j in I \ i
 -                foreach( const Neighbor &J, nbV(j) )
 -                    if( J != I ) {     // for all J in nb(j) \ I 
 -                        prod *= Factor( var(j), message(j, J.iter) );
 -                    }
 -            }
 -        newMessage(i,_I) = prod.marginal( var(i) ).p();
 -    } else {
 -        /* OPTIMIZED VERSION */
 -        Prob prod( factor(I).p() );
 -        if( props.logdomain ) 
 -            prod.takeLog();
 +    Factor Fprod( factor(I) );
 +    Prob &prod = Fprod.p();
 +    if( props.logdomain ) 
 +        prod.takeLog();
 +
 +    // Calculate product of incoming messages and factor I
 +    foreach( const Neighbor &j, nbF(I) )
 +        if( j != i ) { // for all j in I \ i
 +            // prod_j will be the product of messages coming into j
 +            Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 ); 
 +            foreach( const Neighbor &J, nbV(j) )
 +                if( J != I ) { // for all J in nb(j) \ I 
 +                    if( props.logdomain )
 +                        prod_j += message( j, J.iter );
 +                    else
 +                        prod_j *= message( j, J.iter );
 +                }
  
 -        // Calculate product of incoming messages and factor I
 -        foreach( const Neighbor &j, nbF(I) ) {
 -            if( j != i ) {     // for all j in I \ i
 +            // multiply prod with prod_j
 +            if( !DAI_BP_FAST ) {
 +                /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
 +                if( props.logdomain )
 +                    Fprod += Factor( var(j), prod_j );
 +                else
 +                    Fprod *= Factor( var(j), prod_j );
 +            } else {
 +                /* OPTIMIZED VERSION */
                  size_t _I = j.dual;
                  // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
                  const ind_t &ind = index(j, _I);
 -
 -                // prod_j will be the product of messages coming into j
 -                Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 ); 
 -                foreach( const Neighbor &J, nbV(j) )
 -                    if( J != I ) { // for all J in nb(j) \ I 
 -                        if( props.logdomain )
 -                            prod_j += message( j, J.iter );
 -                        else
 -                            prod_j *= message( j, J.iter );
 -                    }
 -
 -                // multiply prod with prod_j
                  for( size_t r = 0; r < prod.size(); ++r )
                      if( props.logdomain )
                          prod[r] += prod_j[ind[r]];
                          prod[r] *= prod_j[ind[r]];
              }
          }
 -        if( props.logdomain ) {
 -            prod -= prod.max();
 -            prod.takeExp();
 -        }
  
 -        // Marginalize onto i
 -        Prob marg( var(i).states(), 0.0 );
 +    if( props.logdomain ) {
 +        prod -= prod.max();
 +        prod.takeExp();
 +    }
 +
 +    // Marginalize onto i
 +    Prob marg;
 +    if( !DAI_BP_FAST ) {
 +        /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
 +        if( props.inference == Properties::InfType::SUMPROD ) 
 +            marg = Fprod.marginal( var(i) ).p();
 +        else
 +            marg = Fprod.maxMarginal( var(i) ).p();
 +    } else {
 +        /* OPTIMIZED VERSION */
 +        marg = Prob( var(i).states(), 0.0 );
          // ind is the precalculated IndexFor(i,I) i.e. to x_I == k corresponds x_i == ind[k]
          const ind_t ind = index(i,_I);
          if( props.inference == Properties::InfType::SUMPROD ) 
                  if( prod[r] > marg[ind[r]] ) 
                      marg[ind[r]] = prod[r];
          marg.normalize();
 -
 -        // Store result
 -        if( props.logdomain )
 -            newMessage(i,_I) = marg.log();
 -        else
 -            newMessage(i,_I) = marg;
      }
  
 +    // Store result
 +    if( props.logdomain )
 +        newMessage(i,_I) = marg.log();
 +    else
 +        newMessage(i,_I) = marg;
 +
      // Update the residual if necessary
      if( props.updates == Properties::UpdateType::SEQMAX )
          updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), Prob::DISTLINF ) );
@@@ -353,46 -360,34 +353,46 @@@ void BP::calcBeliefV( size_t i, Prob &
  
  
  void BP::calcBeliefF( size_t I, Prob &p ) const {
 -    p = factor(I).p();
 -    if( props.logdomain )
 -        p.takeLog();
 +    Factor Fprod( factor( I ) );
 +    Prob &prod = Fprod.p();
  
 -    foreach( const Neighbor &j, nbF(I) ) {
 -        size_t _I = j.dual;
 -        // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
 -        const ind_t & ind = index(j, _I);
 +    if( props.logdomain ) 
 +        prod.takeLog();
  
 +    foreach( const Neighbor &j, nbF(I) ) {
          // prod_j will be the product of messages coming into j
 -        Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 ); 
 -        foreach( const Neighbor &J, nbV(j) ) {
 +        Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
 +        foreach( const Neighbor &J, nbV(j) )
              if( J != I ) { // for all J in nb(j) \ I 
                  if( props.logdomain )
                      prod_j += newMessage( j, J.iter );
                  else
                      prod_j *= newMessage( j, J.iter );
              }
 -        }
  
 -        // multiply p with prod_j
 -        for( size_t r = 0; r < p.size(); ++r ) {
 +        // multiply prod with prod_j
 +        if( !DAI_BP_FAST ) {
 +            /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
              if( props.logdomain )
 -                p[r] += prod_j[ind[r]];
 +                Fprod += Factor( var(j), prod_j );
              else
 -                p[r] *= prod_j[ind[r]];
 +                Fprod *= Factor( var(j), prod_j );
 +        } else {
 +            /* OPTIMIZED VERSION */
 +            size_t _I = j.dual;
 +            // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
 +            const ind_t & ind = index(j, _I);
 +
 +            for( size_t r = 0; r < prod.size(); ++r ) {
 +                if( props.logdomain )
 +                    prod[r] += prod_j[ind[r]];
 +                else
 +                    prod[r] *= prod_j[ind[r]];
 +            }
          }
      }
 +
 +    p = prod;
  }
  
  
@@@ -404,26 -399,12 +404,26 @@@ Factor BP::beliefV( size_t i ) const 
          p -= p.max();
          p.takeExp();
      }
 -
      p.normalize();
 +
      return( Factor( var(i), p ) );
  }
  
  
 +Factor BP::beliefF( size_t I ) const {
 +    Prob p;
 +    calcBeliefF( I, p );
 +
 +    if( props.logdomain ) {
 +        p -= p.max();
 +        p.takeExp();
 +    }
 +    p.normalize();
 +
 +    return( Factor( factor(I).vars(), p ) );
 +}
 +
 +
  Factor BP::belief( const Var &n ) const {
      return( beliefV( findVar( n ) ) );
  }
@@@ -453,6 -434,37 +453,6 @@@ Factor BP::belief( const VarSet &ns ) c
  }
  
  
 -Factor BP::beliefF( size_t I ) const {
 -    if( !DAI_BP_FAST ) {
 -        /*  UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
 -
 -        Factor prod( factor(I) );
 -        foreach( const Neighbor &j, nbF(I) ) {
 -            foreach( const Neighbor &J, nbV(j) ) {
 -                if( J != I )  // for all J in nb(j) \ I
 -                    prod *= Factor( var(j), newMessage(j, J.iter) );
 -            }
 -        }
 -        return prod.normalized();
 -    } else {
 -        /* OPTIMIZED VERSION */
 -
 -        Prob prod;
 -        calcBeliefF( I, prod );
 -
 -        if( props.logdomain ) {
 -            prod -= prod.max();
 -            prod.takeExp();
 -        }
 -        prod.normalize();
 -
 -        Factor result( factor(I).vars(), prod );
 -
 -        return( result );
 -    }
 -}
 -
 -
  Real BP::logZ() const {
      Real sum = 0.0;
      for(size_t i = 0; i < nrVars(); ++i )
@@@ -483,6 -495,8 +483,8 @@@ void BP::init( const VarSet &ns ) 
  
  
  void BP::updateMessage( size_t i, size_t _I ) {
+     if( recordSentMessages )
+         _sentMessages.push_back(make_pair(i,_I));
      if( props.damping == 0.0 ) {
          message(i,_I) = newMessage(i,_I);
          if( props.updates == Properties::UpdateType::SEQMAX )
diff --combined src/factorgraph.cpp
@@@ -21,7 -21,6 +21,7 @@@
  
  
  #include <iostream>
 +#include <iomanip>
  #include <iterator>
  #include <map>
  #include <set>
@@@ -101,8 -100,11 +101,8 @@@ ostream& operator << (ostream& os, cons
                  nr_nonzeros++;
          os << nr_nonzeros << endl;
          for( size_t k = 0; k < fg.factor(I).states(); k++ )
 -            if( fg.factor(I)[k] != 0.0 ) {
 -                char buf[20];
 -                sprintf(buf,"%18.14g", fg.factor(I)[k]);
 -                os << k << " " << buf << endl;
 -            }
 +            if( fg.factor(I)[k] != 0.0 )
 +                os << k << " " << setw(os.precision()+4) << fg.factor(I)[k] << endl;
      }
  
      return(os);
@@@ -266,11 -268,10 +266,11 @@@ void FactorGraph::ReadFromFile( const c
  }
  
  
 -void FactorGraph::WriteToFile( const char *filename ) const {
 +void FactorGraph::WriteToFile( const char *filename, size_t precision ) const {
      ofstream outfile;
      outfile.open( filename );
      if( outfile.is_open() ) {
 +        outfile.precision( precision );
          outfile << *this;
          outfile.close();
      } else
@@@ -330,6 -331,37 +330,37 @@@ void FactorGraph::clamp( const Var & n
  }
  
  
+ void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
+     Var n = var(i);
+     Factor mask_n( n, 0.0 );
+     foreach( size_t i, is ) {
+         assert( i <= n.states() );
+         mask_n[i] = 1.0;
+     }
+     map<size_t, Factor> newFacs;
+     for( size_t I = 0; I < nrFactors(); I++ ) 
+         if( factor(I).vars().contains( n ) ) {
+             newFacs[I] = factor(I) * mask_n;
+         }
+     setFactors( newFacs, backup );
+ }
+ void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup ) {
+     size_t st = factor(I).states();
+     Factor newF( factor(I).vars(), 0.0 );
+     foreach( size_t i, is ) { 
+         assert( i <= st ); 
+         newF[i] = factor(I)[i];
+     }
+     setFactor( I, newF, backup );
+ }
  void FactorGraph::backupFactor( size_t I ) {
      map<size_t,Factor>::iterator it = _backup.find( I );
      if( it != _backup.end() )
@@@ -372,6 -404,7 +403,7 @@@ void FactorGraph::restoreFactors() 
      _backup.clear();
  }
  
  void FactorGraph::backupFactors( const std::set<size_t> & facs ) {
      for( std::set<size_t>::const_iterator fac = facs.begin(); fac != facs.end(); fac++ )
          backupFactor( *fac );
diff --combined tests/testdai.cpp
@@@ -197,7 -197,7 +197,7 @@@ pair<string, PropertySet> parseMethod( 
          if( ps.first == DAINames[n] )
              break;
      if( strlen( DAINames[n] ) == 0 && (ps.first != "LDPC") )
-         throw string("Unknown DAI algorithm \"") + ps.first + string("\" in \"") + _s + string("\"");
+         throw std::runtime_error(string("Unknown DAI algorithm \"") + ps.first + string("\" in \"") + _s + string("\""));
  
      return ps;
  }
@@@ -253,132 -253,124 +253,132 @@@ int main( int argc, char *argv[] ) 
          cout << "error and relative logZ error (comparing with the results of" << endl;
          cout << "<method0>, the base method, for which one can use JTREE_HUGIN)." << endl << endl;
          cout << opts_required << opts_optional << endl;
 +#ifdef DAI_DEBUG
 +        cout << "This is a debugging (unoptimised) build of libDAI." << endl;
 +#endif
          return 1;
      }
  
 -    // Read aliases
 -    map<string,string> Aliases;
 -    if( !aliases.empty() ) {
 -        ifstream infile;
 -        infile.open (aliases.c_str());
 -        if (infile.is_open()) {
 -            while( true ) {
 -                string line;
 -                getline(infile,line);
 -                if( infile.fail() )
 -                    break;
 -                if( (!line.empty()) && (line[0] != '#') ) {
 -                    string::size_type pos = line.find(':',0);
 -                    if( pos == string::npos )
 -                        throw "Invalid alias";
 -                    else {
 -                        string::size_type posl = line.substr(0, pos).find_last_not_of(" \t");
 -                        string key = line.substr(0, posl + 1);
 -                        string::size_type posr = line.substr(pos + 1, line.length()).find_first_not_of(" \t");
 -                        string val = line.substr(pos + 1 + posr, line.length());
 -                        Aliases[key] = val;
 +    try {
 +        // Read aliases
 +        map<string,string> Aliases;
 +        if( !aliases.empty() ) {
 +            ifstream infile;
 +            infile.open (aliases.c_str());
 +            if (infile.is_open()) {
 +                while( true ) {
 +                    string line;
 +                    getline(infile,line);
 +                    if( infile.fail() )
 +                        break;
 +                    if( (!line.empty()) && (line[0] != '#') ) {
 +                        string::size_type pos = line.find(':',0);
 +                        if( pos == string::npos )
 +                            throw string("Invalid alias");
 +                        else {
 +                            string::size_type posl = line.substr(0, pos).find_last_not_of(" \t");
 +                            string key = line.substr(0, posl + 1);
 +                            string::size_type posr = line.substr(pos + 1, line.length()).find_first_not_of(" \t");
 +                            string val = line.substr(pos + 1 + posr, line.length());
 +                            Aliases[key] = val;
 +                        }
                      }
                  }
 -            }
 -            infile.close();
 -        } else
 -            throw "Error opening aliases file";
 -    }
 -
 -    FactorGraph fg;
 -    fg.ReadFromFile( filename.c_str() );
 -
 -    vector<Factor> q0;
 -    double logZ0 = 0.0;
 -
 -    cout.setf( ios_base::scientific );
 -    cout.precision( 3 );
 -
 -    cout << "# " << filename << endl;
 -    cout.width( 39 );
 -    cout << left << "# METHOD" << "\t";
 -    if( report_time )
 -        cout << right << "SECONDS  " << "\t";
 -    if( report_iters )
 -        cout << "ITERS" << "\t";
 -    cout << "MAX ERROR" << "\t";
 -    cout << "AVG ERROR" << "\t";
 -    cout << "LOGZ ERROR" << "\t";
 -    cout << "MAXDIFF" << "\t";
 -    cout << endl;
 -
 -    for( size_t m = 0; m < methods.size(); m++ ) {
 -        pair<string, PropertySet> meth = parseMethod( methods[m], Aliases );
 -
 -        if( vm.count("tol") )
 -            meth.second.Set("tol",tol);
 -        if( vm.count("maxiter") )
 -            meth.second.Set("maxiter",maxiter);
 -        if( vm.count("verbose") )
 -            meth.second.Set("verbose",verbose);
 -        TestDAI piet(fg, meth.first, meth.second );
 -        piet.doDAI();
 -        if( m == 0 ) {
 -            q0 = piet.q;
 -            logZ0 = piet.logZ;
 +                infile.close();
 +            } else
 +                throw string("Error opening aliases file");
          }
 -        piet.calcErrs(q0);
  
 +        FactorGraph fg;
 +        fg.ReadFromFile( filename.c_str() );
 +
 +        vector<Factor> q0;
 +        double logZ0 = 0.0;
 +
 +        cout.setf( ios_base::scientific );
 +        cout.precision( 3 );
 +
 +        cout << "# " << filename << endl;
          cout.width( 39 );
 -        cout << left << methods[m] << "\t";
 +        cout << left << "# METHOD" << "\t";
          if( report_time )
 -            cout << right << piet.time << "\t";
 -        if( report_iters ) {
 -            if( piet.has_iters ) {
 -                cout << piet.iters << "\t";
 -            } else {
 -                cout << "N/A  \t";
 +            cout << right << "SECONDS  " << "\t";
 +        if( report_iters )
 +            cout << "ITERS" << "\t";
 +        cout << "MAX ERROR" << "\t";
 +        cout << "AVG ERROR" << "\t";
 +        cout << "LOGZ ERROR" << "\t";
 +        cout << "MAXDIFF" << "\t";
 +        cout << endl;
 +
 +        for( size_t m = 0; m < methods.size(); m++ ) {
 +            pair<string, PropertySet> meth = parseMethod( methods[m], Aliases );
 +
 +            if( vm.count("tol") )
 +                meth.second.Set("tol",tol);
 +            if( vm.count("maxiter") )
 +                meth.second.Set("maxiter",maxiter);
 +            if( vm.count("verbose") )
 +                meth.second.Set("verbose",verbose);
 +            TestDAI piet(fg, meth.first, meth.second );
 +            piet.doDAI();
 +            if( m == 0 ) {
 +                q0 = piet.q;
 +                logZ0 = piet.logZ;
 +            }
 +            piet.calcErrs(q0);
 +
 +            cout.width( 39 );
 +            cout << left << methods[m] << "\t";
 +            if( report_time )
 +                cout << right << piet.time << "\t";
 +            if( report_iters ) {
 +                if( piet.has_iters ) {
 +                    cout << piet.iters << "\t";
 +                } else {
 +                    cout << "N/A  \t";
 +                }
              }
 -        }
  
 -        if( m > 0 ) {
 -            cout.setf( ios_base::scientific );
 -            cout.precision( 3 );
 -            
 -            double me = clipdouble( piet.maxErr(), 1e-9 );
 -            cout << me << "\t";
 -            
 -            double ae = clipdouble( piet.avgErr(), 1e-9 );
 -            cout << ae << "\t";
 -            
 -            if( piet.has_logZ ) {
 -                cout.setf( ios::showpos );
 -                double le = clipdouble( piet.logZ / logZ0 - 1.0, 1e-9 );
 -                cout << le << "\t";
 -                cout.unsetf( ios::showpos );
 -            } else
 -                cout << "N/A       \t";
 -
 -            if( piet.has_maxdiff ) {
 -                double md = clipdouble( piet.maxdiff, 1e-9 );
 -                if( isnan( me ) )
 -                    md = me;
 -                if( isnan( ae ) )
 -                    md = ae;
 -                cout << md << "\t";
 -            } else
 -                cout << "N/A    \t";
 -        }
 -        cout << endl;
 +            if( m > 0 ) {
 +                cout.setf( ios_base::scientific );
 +                cout.precision( 3 );
 +                
 +                double me = clipdouble( piet.maxErr(), 1e-9 );
 +                cout << me << "\t";
 +                
 +                double ae = clipdouble( piet.avgErr(), 1e-9 );
 +                cout << ae << "\t";
 +                
 +                if( piet.has_logZ ) {
 +                    cout.setf( ios::showpos );
 +                    double le = clipdouble( piet.logZ / logZ0 - 1.0, 1e-9 );
 +                    cout << le << "\t";
 +                    cout.unsetf( ios::showpos );
 +                } else
 +                    cout << "N/A       \t";
 +
 +                if( piet.has_maxdiff ) {
 +                    double md = clipdouble( piet.maxdiff, 1e-9 );
 +                    if( isnan( me ) )
 +                        md = me;
 +                    if( isnan( ae ) )
 +                        md = ae;
 +                    cout << md << "\t";
 +                } else
 +                    cout << "N/A    \t";
 +            }
 +            cout << endl;
  
 -        if( marginals ) {
 -            for( size_t i = 0; i < piet.q.size(); i++ )
 -                cout << "# " << piet.q[i] << endl;
 +            if( marginals ) {
 +                for( size_t i = 0; i < piet.q.size(); i++ )
 +                    cout << "# " << piet.q[i] << endl;
 +            }
          }
 -    }
  
 -    return 0;
 +        return 0;
 +    } catch( string &s ) {
 +        cerr << "Exception: " << s << endl;
 +        return 2;
 +    }
  }