Merge branch 'no-edges2'
authorJoris Mooij <jorism@marvin.jorismooij.nl>
Mon, 8 Sep 2008 14:18:16 +0000 (16:18 +0200)
committerJoris Mooij <jorism@marvin.jorismooij.nl>
Mon, 8 Sep 2008 14:18:16 +0000 (16:18 +0200)
27 files changed:
Makefile
include/dai/alldai.h
include/dai/bipgraph.h
include/dai/bp.h
include/dai/daialg.h
include/dai/factorgraph.h
include/dai/hak.h
include/dai/jtree.h
include/dai/lc.h
include/dai/mf.h
include/dai/regiongraph.h
include/dai/treeep.h
include/dai/util.h
src/alldai.cpp
src/bp.cpp
src/factorgraph.cpp
src/hak.cpp
src/jtree.cpp
src/lc.cpp
src/mf.cpp
src/mr.cpp
src/properties.cpp
src/regiongraph.cpp
src/treeep.cpp
utils/createfg.cpp
utils/fg2dot.cpp
utils/fginfo.cpp

index f845617..a3e5c35 100644 (file)
--- a/Makefile
+++ b/Makefile
 #   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 
 
+# Enable/disable various approximate inference methods
+WITH_BP = true
+WITH_MF = true
+WITH_HAK = true
+WITH_LC = true
+WITH_TREEEP = true
+WITH_JTREE = true
+WITH_MR = true
+DEBUG = true
+NEW_MATLAB = true
+WITH_MATLAB =
+
 # Directories
 INC = include/dai
 SRC = src
@@ -30,24 +42,62 @@ BOOSTFLAGS = -lboost_program_options
 CC = g++
 
 # Flags for the C++ compiler
-CCFLAGS = -Wall -W -Wextra -fpic -I./include -Llib -O3 -g -DDAI_DEBUG #-static #-pg #-DVERBOSE
+CCFLAGS = -Wall -W -Wextra -fpic -I./include -Llib -O3 #-pg #-static -DVERBOSE
+ifdef DEBUG
+CCFLAGS := $(CCFLAGS) -g -DDAI_DEBUG
+else
+CCFLAGS := $(CCFLAGS)
+endif
+
+ifdef WITH_BP
+CCFLAGS := $(CCFLAGS) -DWITH_BP
+OBJECTS := $(OBJECTS) bp.o
+endif
+ifdef WITH_MF
+CCFLAGS := $(CCFLAGS) -DWITH_MF
+OBJECTS := $(OBJECTS) mf.o
+endif
+ifdef WITH_HAK
+CCFLAGS := $(CCFLAGS) -DWITH_HAK
+OBJECTS := $(OBJECTS) hak.o
+endif
+ifdef WITH_LC
+CCFLAGS := $(CCFLAGS) -DWITH_LC
+OBJECTS := $(OBJECTS) lc.o
+endif
+ifdef WITH_TREEEP
+CCFLAGS := $(CCFLAGS) -DWITH_TREEEP
+OBJECTS := $(OBJECTS) treeep.o
+endif
+ifdef WITH_JTREE
+CCFLAGS := $(CCFLAGS) -DWITH_JTREE
+OBJECTS := $(OBJECTS) jtree.o
+endif
+ifdef WITH_MR
+CCFLAGS := $(CCFLAGS) -DWITH_MR
+OBJECTS := $(OBJECTS) mr.o
+endif
 
-# To enable the Matlab interface, define WITH_MATLAB = yes
-WITH_MATLAB = 
 ifdef WITH_MATLAB
 # Replace the following by the directory where Matlab has been installed
 MATLABDIR = /opt/matlab/bin
-MEX = $(MATLABDIR)/mex
-MEXFLAGS = -g -I. -DDAI_DEBUG -largeArrayDims #-g means debugging
-endif
-
 # Replace the following with the extension of compiled MEX files on this platform, e.g. .mexglx for x86
 MEXEXT = .mexglx
+MEX = $(MATLABDIR)/mex
+MEXFLAGS = -I.
+ifdef DEBUG
+MEXFLAGS := $(MEXFLAGS) -g -DDAI_DEBUG
+endif
+ifdef NEW_MATLAB
+MEXFLAGS := $(MEXFLAGS) -largeArrayDims
+else
+MEXFLAGS := $(MEXFLAGS) -DSMALLMEM
+endif
+endif
 
 HEADERS = $(INC)/bipgraph.h $(INC)/diffs.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/x2x.h
 
-# target matlabs is disabled by default since it only compiles with a very recent MatLab version
-TARGETS = tests utils $(LIB)/libdai.a example testregression
+TARGETS = tests utils $(LIB)/libdai.a example testregression doc
 ifdef WITH_MATLAB
 TARGETS := $(TARGETS) matlabs
 endif
@@ -56,8 +106,8 @@ all : $(TARGETS)
 
 matlabs : matlab/dai.$(MEXEXT) matlab/dai_readfg.$(MEXEXT) matlab/dai_writefg.$(MEXEXT) matlab/dai_removeshortloops.$(MEXEXT) matlab/dai_potstrength.$(MEXEXT)
 
-$(LIB)/libdai.a : daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
-       ar rcs $(LIB)/libdai.a daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
+$(LIB)/libdai.a : daialg.o alldai.o clustergraph.o factorgraph.o properties.o regiongraph.o util.o weightedgraph.o x2x.o $(OBJECTS)
+       ar rcs $(LIB)/libdai.a daialg.o alldai.o clustergraph.o factorgraph.o properties.o regiongraph.o util.o weightedgraph.o x2x.o $(OBJECTS)
 
 tests : tests/test
 
@@ -65,13 +115,13 @@ utils : utils/createfg utils/fg2dot utils/remove_short_loops utils/fginfo
 
 testregression : tests/test
        echo Testing...this can take a while...
-       cd tests; ./testregression; cd ..
+       cd tests; time ./testregression; cd ..
 
 doc : $(INC)/*.h $(SRC)/*.cpp doxygen.conf
        doxygen doxygen.conf
 
 clean :
-       rm *.o *.$(MEXEXT) example matlab/*.$(MEXEXT) matlab/*.o tests/test utils/fg2dot utils/createfg utils/remove_short_loops utils/fginfo $(LIB)/libdai.a; echo
+       rm *.o example matlab/*.$(MEXEXT) matlab/*.o tests/test utils/fg2dot utils/createfg utils/remove_short_loops utils/fginfo $(LIB)/libdai.a; echo
        rm -R doc; echo
 
 
@@ -132,15 +182,15 @@ example : $(SRC)/example.cpp $(HEADERS) $(LIB)/libdai.a
 # TESTS
 ########
 
-tests/test : tests/test.cpp $(HEADERS) lib/libdai.a
+tests/test : tests/test.cpp $(HEADERS) $(LIB)/libdai.a
        $(CC) $(CCFLAGS) -o tests/test tests/test.cpp -ldai $(BOOSTFLAGS)
 
 
 # MATLAB INTERFACE
 ###################
 
-matlab/dai.$(MEXEXT) : matlab/dai.cpp $(HEADERS) matlab/matlab.o daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
-       $(MEX) $(MEXFLAGS) -o matlab/dai matlab/dai.cpp matlab/matlab.o daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
+matlab/dai.$(MEXEXT) : matlab/dai.cpp $(HEADERS) matlab/matlab.o $(LIB)/libdai.a
+       $(MEX) $(MEXFLAGS) -o matlab/dai matlab/dai.cpp matlab/matlab.o $(LIB)/libdai.a
 
 matlab/dai_readfg.$(MEXEXT) : matlab/dai_readfg.cpp $(HEADERS) factorgraph.o matlab/matlab.o
        $(MEX) $(MEXFLAGS) -o matlab/dai_readfg matlab/dai_readfg.cpp factorgraph.o matlab/matlab.o
index 111a9b6..87bb1d8 100644 (file)
 
 #include <string>
 #include <dai/daialg.h>
-#include <dai/bp.h>
-#include <dai/lc.h>
-#include <dai/hak.h>
-#include <dai/mf.h>
-#include <dai/jtree.h>
-#include <dai/treeep.h>
-#include <dai/mr.h>
+#ifdef WITH_BP
+    #include <dai/bp.h>
+#endif
+#ifdef WITH_MF
+    #include <dai/mf.h>
+#endif
+#ifdef WITH_HAK
+    #include <dai/hak.h>
+#endif
+#ifdef WITH_LC
+    #include <dai/lc.h>
+#endif
+#ifdef WITH_TREEEP
+    #include <dai/treeep.h>
+#endif
+#ifdef WITH_JTREE
+    #include <dai/jtree.h>
+#endif
+#ifdef WITH_MR
+    #include <dai/mr.h>
+#endif
 
 
 namespace dai {
@@ -43,8 +57,32 @@ namespace dai {
 InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const Properties &opts );
 
 
-/// AINames contains the names of all approximate inference algorithms
-static const char* DAINames[] = {BP::Name, MF::Name, HAK::Name, LC::Name, TreeEP::Name, MR::Name, JTree::Name};
+/// DAINames contains the names of all approximate inference algorithms
+
+static const char* DAINames[] = {
+#ifdef WITH_BP
+    BP::Name, 
+#endif
+#ifdef WITH_MF
+    MF::Name,
+#endif
+#ifdef WITH_HAK
+    HAK::Name,
+#endif
+#ifdef WITH_LC
+    LC::Name,
+#endif
+#ifdef WITH_TREEEP
+    TreeEP::Name,
+#endif
+#ifdef WITH_JTREE
+    JTree::Name,
+#endif
+#ifdef WITH_MR
+    MR::Name,
+#endif
+    ""
+};
 
 
 } // end of namespace dai
index d7decba..ef8e4cb 100644 (file)
 
 #include <vector>
 #include <algorithm>
-#include <boost/numeric/ublas/vector_sparse.hpp>
+#include <dai/util.h>
 
 
 namespace dai {
 
 
-/// A BipartiteGraph represents a graph with two types of nodes
-template <class T1, class T2> class BipartiteGraph {
+/// A BipartiteGraph represents a bipartite graph, with two types of nodes (both are numbered
+/// as 0,1,2,...), with edges only between nodes of different type. The edges are stored as 
+/// lists of adjacent nodes for each node.
+class BipartiteGraph {
     public:
-        typedef std::pair<size_t,size_t>    _edge_t;
-        typedef std::vector<size_t>         _nb_t;
-        typedef _nb_t::const_iterator       _nb_cit;
-        
+        /// A Neighbor describes a neighboring node of some other node.
+        /** Iterating over all neighbors of some node i can be done in the following way:
+         *  \code
+         *      foreach( const BipartiteGraph::Neighbor &I, nb1(i) ) {
+         *          size_t _I = I.iter;
+         *          size_t _i = I.dual;
+         *          // I == I.node;
+         *          // The _I'th neighbor of i is I, and the _i'th neighbor of I is i:
+         *          // nb1(i)[_I] == I, nb2(I)[_i] == i
+         *      }
+         *  \endcode
+         */
+        struct Neighbor {
+            /// iter corresponds to the index of this Neighbor entry in the list of neighbors
+            unsigned iter;
+            /// node contains the number of the neighboring node
+            unsigned node;
+            /// dual contains the "dual" iter
+            unsigned dual;
+            /// cast to unsigned returns node
+            operator unsigned () const { return node; };
+        };
+
+        /// Neighbors is a vector of Neigbor entries
+        typedef std::vector<Neighbor> Neighbors;
 
     private:
-        /// Vertices of first type
-        std::vector<T1>                     _V1;
-
-        /// Vertices of second type
-        std::vector<T2>                     _V2;
-        
-        /// Edges, which are pairs (v1,v2) with v1 in _V1 and v2 in _V2
-        std::vector<_edge_t>                _E12;
-
-        /// Conversion matrix that computes which index of _E12 corresponds to a vertex index pair (v1,v2)
-        std::vector<boost::numeric::ublas::mapped_vector<size_t> > _E12ind;
-
-        /// Neighbour indices of vertices of first type
-        std::vector<_nb_t>                  _nb1;
-
-        /// Neighbour indices of vertices of second type
-        std::vector<_nb_t>                  _nb2;
-
+        /// _nb1 contains for each node of the first kind a list of its neighbors
+        std::vector<Neighbors> _nb1;
+        /// _nb2 contains for each node of the second kind a list of its neighbors
+        std::vector<Neighbors> _nb2;
 
     public:
         /// Default constructor
-        BipartiteGraph<T1,T2> () {};
+        BipartiteGraph() : _nb1(), _nb2() {}
 
         /// Copy constructor
-        BipartiteGraph<T1,T2> ( const BipartiteGraph<T1,T2> & x ) : _V1(x._V1), _V2(x._V2), _E12(x._E12), _E12ind(x._E12ind), _nb1(x._nb1), _nb2(x._nb2) {};
+        BipartiteGraph( const BipartiteGraph & x ) : _nb1(x._nb1), _nb2(x._nb2) {}
 
         /// Assignment operator
-        BipartiteGraph<T1,T2> & operator=(const BipartiteGraph<T1,T2> & x) {
+        BipartiteGraph & operator=( const BipartiteGraph & x ) {
             if( this != &x ) {
-                _V1 =       x._V1;
-                _V2 =       x._V2;
-                _E12 =      x._E12;
-                _E12ind =   x._E12ind;
-                _nb1 =      x._nb1;
-                _nb2 =      x._nb2;
+                _nb1 = x._nb1;
+                _nb2 = x._nb2;
             }
             return *this;
         }
+
+        /// Create bipartite graph from a range of edges, encoded as pairs of node numbers
+        /// (more precisely, a std::pair<unsigned, unsigned> where the first integer corresponds
+        /// to the node of the first type, and the second integer corresponds to the node of the
+        /// second type). nr1 is the number of nodes of the first type, nr2 the number of nodes
+        /// of the second type.
+        template<typename EdgeInputIterator>
+        void create( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end );
+
+        /// Construct bipartite graph from a range of edges, encoded as pairs of node numbers
+        /// (more precisely, a std::pair<unsigned, unsigned> where the first integer corresponds
+        /// to the node of the first type, and the second integer corresponds to the node of the
+        /// second type). nr1 is the number of nodes of the first type, nr2 the number of nodes
+        /// of the second type.
+        template<typename EdgeInputIterator>
+        BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ) {
+            create( nr1, nr2, begin, end );
+        }
+
+        /// Returns constant reference to the _i2'th neighbor of node i1 of first type
+        const Neighbor & nb1( size_t i1, size_t _i2 ) const { return _nb1[i1][_i2]; }
+        /// Returns reference to the _i2'th neighbor of node i1 of first type
+        Neighbor & nb1( size_t i1, size_t _i2 ) { return _nb1[i1][_i2]; }
+
+        /// Returns constant reference to the _i1'th neighbor of node i2 of second type
+        const Neighbor & nb2( size_t i2, size_t _i1 ) const { return _nb2[i2][_i1]; }
+        /// Returns reference to the _i1'th neighbor of node i2 of second type
+        Neighbor & nb2( size_t i2, size_t _i1 ) { return _nb2[i2][_i1]; }
+
+        /// Returns constant reference to all neighbors of node of first type
+        const Neighbors & nb1( size_t i1 ) const { return _nb1[i1]; }
+        /// Returns reference to all neighbors of node of first type
+        Neighbors & nb1( size_t i1 ) { return _nb1[i1]; }
+
+        /// Returns constant reference to all neighbors of node of second type
+        const Neighbors & nb2( size_t i2 ) const { return _nb2[i2]; }
+        /// Returns reference to all neighbors of node of second type
+        Neighbors & nb2( size_t i2 ) { return _nb2[i2]; }
+
+        /// Returns number of nodes of first type
+        size_t nr1() const { return _nb1.size(); }
+        /// Returns number of nodes of second type
+        size_t nr2() const { return _nb2.size(); }
         
-        /// Provides read access to node of first type
-        const T1 & V1( size_t i1 ) const { return _V1[i1]; }
-        /// Provides full access to node of first type
-        T1 & V1( size_t i1 ) { return _V1[i1]; }
-        /// Provides read access to all nodes of first type
-        const std::vector<T1> & V1s() const { return _V1; }
-        /// Provides full access to all nodes of first type
-        std::vector<T1> & V1s() { return _V1; }
-
-        /// Provides read access to node of second type
-        const T2 & V2( size_t i2 ) const { return _V2[i2]; }
-        /// Provides full access to node of second type
-        T2 & V2( size_t i2 ) { return _V2[i2]; }
-        /// Provides read access to all nodes of second type
-        const std::vector<T2> & V2s() const { return _V2; }
-        /// Provides full access to all nodes of second type
-        std::vector<T2> & V2s() { return _V2; }
-
-        /// Provides read access to edge
-        const _edge_t & edge(size_t ind) const { return _E12[ind]; }
-        /// Provides full access to edge
-        _edge_t & edge(size_t ind) { return _E12[ind]; }
-        /// Provides read access to all edges
-        const std::vector<_edge_t> & edges() const { return _E12; }
-        /// Provides full access to all edges
-        std::vector<_edge_t> & edges() { return _E12; }
-        /// Returns number of edges
-        size_t nr_edges() const { return _E12.size(); }
-
-        /// Provides read access to neighbours of node of first type
-        const _nb_t & nb1( size_t i1 ) const { return _nb1[i1]; }
-        /// Provides full access to neighbours of node of first type
-        _nb_t & nb1( size_t i1 ) { return _nb1[i1]; }
-
-        /// Provides read access to neighbours of node of second type
-        const _nb_t & nb2( size_t i2 ) const { return _nb2[i2]; }
-        /// Provides full access to neighbours of node of second type
-        _nb_t & nb2( size_t i2 ) { return _nb2[i2]; }
-
-        /// Converts the pair of indices (i1,i2) to the corresponding edge index
-        size_t VV2E( const size_t i1, const size_t i2 ) const { return _E12ind[i1](i2); }
-
-
-        /// Regenerates internal structures (_E12ind, _nb1 and _nb2) based upon _V1, _V2 and _E12
-        void Regenerate() {
-            // Calculate _nb1 and _nb2
-
-            // Start with empty vectors
-            _nb1.clear();
-            _nb1.resize(_V1.size());
-            // Start with empty vectors
-            _nb2.clear();
-            _nb2.resize(_V2.size());
-            // Each edge yields a neighbour pair
-            for( std::vector<_edge_t>::const_iterator e = _E12.begin(); e != _E12.end(); e++ ) {
-                _nb1[e->first].push_back(e->second);
-                _nb2[e->second].push_back(e->first);
-            }
-            // Remove duplicates from _nb1
-            for( size_t i1 = 0; i1 < _V1.size(); i1++ ) {
-                _nb_t::iterator new_end = unique(_nb1[i1].begin(), _nb1[i1].end());
-                _nb1[i1].erase( new_end, _nb1[i1].end() );
-            }
-            // Remove duplicates from _nb2
-            for( size_t i2 = 0; i2 < _V2.size(); i2++ ) {
-                _nb_t::iterator new_end = unique(_nb2[i2].begin(), _nb2[i2].end());
-                _nb2[i2].erase( new_end, _nb2[i2].end() );
-            }
+        /// Calculates the number of edges
+        size_t nrEdges() const {
+            size_t sum = 0;
+            for( size_t i1 = 0; i1 < nr1(); i1++ )
+                sum += nb1(i1).size();
+            return sum;
+        }
 
-            // Calculate _E12ind
-            _E12ind = std::vector<boost::numeric::ublas::mapped_vector<size_t> >(_V1.size(), boost::numeric::ublas::mapped_vector<size_t>(_V2.size()) );            
-            for( size_t k = 0; k < _E12.size(); k++ )
-                _E12ind[_E12[k].first](_E12[k].second) = k;
-                    
+        /// Returns true if the graph is connected
+        /// FIXME: this should be optimized
+        bool isConnected() const {
+            if( nr1() == 0 ) {
+                return true;
+            } else {
+                std::vector<bool> incomponent1( nr1(), false );
+                std::vector<bool> incomponent2( nr2(), false );
+
+                incomponent1[0] = true;
+                bool found_new_nodes;
+                do {
+                    found_new_nodes = false;
+
+                    // For all nodes of second type, check if they are connected with the (growing) component
+                    for( size_t n2 = 0; n2 < nr2(); n2++ )
+                        if( !incomponent2[n2] ) {
+                            foreach( const Neighbor &n1, nb2(n2) ) {
+                                if( incomponent1[n1] ) {
+                                    found_new_nodes = true;
+                                    incomponent2[n2] = true;
+                                    break;
+                                }
+                            }
+                        }
+
+                    // For all nodes of first type, check if they are connected with the (growing) component
+                    for( size_t n1 = 0; n1 < nr1(); n1++ )
+                        if( !incomponent1[n1] ) {
+                            foreach( const Neighbor &n2, nb1(n1) ) {
+                                if( incomponent2[n2] ) {
+                                    found_new_nodes = true;
+                                    incomponent1[n1] = true;
+                                    break;
+                                }
+                            }
+                        }
+                } while( found_new_nodes );
+
+                // Check if there are remaining nodes (not in the component)
+                bool all_connected = true;
+                for( size_t n1 = 0; (n1 < nr1()) && all_connected; n1++ )
+                    if( !incomponent1[n1] )
+                        all_connected = false;
+                for( size_t n2 = 0; (n2 < nr2()) && all_connected; n2++ )
+                    if( !incomponent2[n2] )
+                        all_connected = false;
+
+                return all_connected;
+            }
         }
+
 };
 
 
+template<typename EdgeInputIterator>
+void BipartiteGraph::create( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) {
+    _nb1.clear();
+    _nb1.resize( nr1 );
+    _nb2.clear();
+    _nb2.resize( nr2 );
+
+    for( EdgeInputIterator e = begin; e != end; e++ ) {
+        // Each edge yields a neighbor pair
+        Neighbor nb_1;
+        nb_1.iter = _nb1[e->first].size();
+        nb_1.node = e->second;
+        nb_1.dual = _nb2[e->second].size();
+
+        Neighbor nb_2;
+        nb_2.iter = nb_1.dual;
+        nb_2.node = e->first;
+        nb_2.dual = nb_1.iter;
+
+        _nb1[e->first].push_back( nb_1 );
+        _nb2[e->second].push_back( nb_2 );
+    }
+}
+
+
 } // end of namespace dai
 
 
index 7b66aef..7ec76f0 100644 (file)
@@ -34,10 +34,15 @@ namespace dai {
 
 class BP : public DAIAlgFG {
     protected:
-        typedef std::vector<size_t>  _ind_t;
+        typedef std::vector<size_t> ind_t;
+        struct EdgeProp {
+            ind_t  index;
+            Prob   message;
+            Prob   newMessage;
+            double residual;
+        };
 
-        std::vector<_ind_t>          _indices;
-        std::vector<Prob>            _messages, _newmessages;
+        std::vector<std::vector<EdgeProp> > edges;
 
     public:
         ENUM4(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL)
@@ -46,7 +51,7 @@ class BP : public DAIAlgFG {
         // default constructor
         BP() : DAIAlgFG() {};
         // copy constructor
-        BP(const BP & x) : DAIAlgFG(x), _indices(x._indices), _messages(x._messages), _newmessages(x._newmessages) {};
+        BP(const BP & x) : DAIAlgFG(x), edges(x.edges) {};
         BP* clone() const { return new BP(*this); }
         // construct BP object from FactorGraph
         BP(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
@@ -57,29 +62,30 @@ class BP : public DAIAlgFG {
         BP & operator=(const BP & x) {
             if(this!=&x) {
                 DAIAlgFG::operator=(x);
-                _messages = x._messages;
-                _newmessages = x._newmessages;
-                _indices = x._indices;
+                edges = x.edges;
             }
             return *this;
         }
 
         static const char *Name;
 
-        Prob & message(size_t i1, size_t i2) { return( _messages[VV2E(i1,i2)] ); }  
-        const Prob & message(size_t i1, size_t i2) const { return( _messages[VV2E(i1,i2)] ); }  
-        Prob & newMessage(size_t i1, size_t i2) { return( _newmessages[VV2E(i1,i2)] ); }    
-        const Prob & newMessage(size_t i1, size_t i2) const { return( _newmessages[VV2E(i1,i2)] ); }    
-        _ind_t & index(size_t i1, size_t i2) { return( _indices[VV2E(i1,i2)] ); }
-        const _ind_t & index(size_t i1, size_t i2) const { return( _indices[VV2E(i1,i2)] ); }
+        Prob & message(size_t i, size_t _I) { return edges[i][_I].message; }
+        const Prob & message(size_t i, size_t _I) const { return edges[i][_I].message; }
+        Prob & newMessage(size_t i, size_t _I) { return edges[i][_I].newMessage; }
+        const Prob & newMessage(size_t i, size_t _I) const { return edges[i][_I].newMessage; }
+        ind_t & index(size_t i, size_t _I) { return edges[i][_I].index; }
+        const ind_t & index(size_t i, size_t _I) const { return edges[i][_I].index; }
+        double & residual(size_t i, size_t _I) { return edges[i][_I].residual; }
+        const double & residual(size_t i, size_t _I) const { return edges[i][_I].residual; }
+        void findMaxResidual( size_t &i, size_t &_I );
 
         std::string identify() const;
         void Regenerate();
         void init();
-        void calcNewMessage(size_t iI);
+        void calcNewMessage( size_t i, size_t _I );
         double run();
-        Factor belief1 (size_t i) const;
-        Factor belief2 (size_t I) const;
+        Factor beliefV (size_t i) const;
+        Factor beliefF (size_t I) const;
         Factor belief (const Var &n) const;
         Factor belief (const VarSet &n) const;
         std::vector<Factor> beliefs() const;
index 46cce56..97d0dc0 100644 (file)
@@ -153,14 +153,11 @@ class InfAlg {
         /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$)
         virtual void clamp( const Var & n, size_t i ) = 0;
 
-        /// Return all variables that interact with n
-        virtual VarSet delta( const Var & n ) const = 0;
+        /// Return all variables that interact with var(i)
+        virtual VarSet delta( size_t i ) const = 0;
 
-        /// Set all factors interacting with n to 1
-        virtual void makeCavity( const Var & n ) = 0;
-
-        /// Set factor I to 1
-        virtual void makeFactorCavity( size_t I ) = 0;
+        /// Set all factors interacting with var(i) to 1
+        virtual void makeCavity( size_t i ) = 0;
 
         /// Get index of variable n
         virtual size_t findVar( const Var & n ) const = 0;
@@ -223,14 +220,11 @@ class DAIAlg : public InfAlg, public T {
         /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$) (using T::clamp)
         void clamp( const Var & n, size_t i ) { T::clamp( n, i ); }
 
-        /// Return all variables that interact with n (using T::delta)
-        VarSet delta( const Var & n ) const { return T::delta(n); }
-
-        /// Set all factors interacting with n to 1 (using T::makeCavity)
-        void makeCavity( const Var & n ) { T::makeCavity(n); }
+        /// Return all variables that interact with var(i) (using T::delta)
+        VarSet delta( size_t i ) const { return T::delta( i ); }
 
-        /// Set factor I to 1 (using T::makeFactorCavity)
-        void makeFactorCavity( size_t I ) { T::makeFactorCavity(I); }
+        /// Set all factors interacting with var(i) to 1 (using T::makeCavity)
+        void makeCavity( size_t i ) { T::makeCavity( i ); }
 
         /// Get index of variable n (using T::findVar)
         size_t findVar( const Var & n ) const { return T::findVar(n); }
index f68ea46..adb6722 100644 (file)
@@ -37,16 +37,23 @@ bool hasShortLoops( const std::vector<Factor> &P );
 void RemoveShortLoops( std::vector<Factor> &P );
 
 
-class FactorGraph : public BipartiteGraph<Var,Factor> {
+class FactorGraph {
+    public:
+        BipartiteGraph         G;
+        std::vector<Var>       vars;
+        std::vector<Factor>    factors;
+        typedef BipartiteGraph::Neighbor  Neighbor;
+        typedef BipartiteGraph::Neighbors Neighbors;
+
     protected:
-        std::map<size_t,Prob>    _undoProbs;
-        Prob::NormType           _normtype;
+        std::map<size_t,Prob>  _undoProbs;
+        Prob::NormType         _normtype;
 
     public:
         /// Default constructor
-        FactorGraph() : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {};
+        FactorGraph() : G(), vars(), factors(), _undoProbs(), _normtype(Prob::NORMPROB) {};
         /// Copy constructor
-        FactorGraph(const FactorGraph & x) : BipartiteGraph<Var,Factor>(x), _undoProbs(), _normtype(x._normtype) {};
+        FactorGraph(const FactorGraph & x) : G(x.G), vars(x.vars), factors(x.factors), _undoProbs(x._undoProbs), _normtype(x._normtype) {};
         /// Construct FactorGraph from vector of Factors
         FactorGraph(const std::vector<Factor> &P);
         // Construct a FactorGraph from given factor and variable iterators
@@ -55,38 +62,46 @@ class FactorGraph : public BipartiteGraph<Var,Factor> {
         
         /// Assignment operator
         FactorGraph & operator=(const FactorGraph & x) {
-            if(this!=&x) {
-                BipartiteGraph<Var,Factor>::operator=(x);
-                _undoProbs      = x._undoProbs;
-                _normtype       = x._normtype;
+            if( this != &x ) {
+                G          = x.G;
+                vars       = x.vars;
+                factors    = x.factors;
+                _undoProbs = x._undoProbs;
+                _normtype  = x._normtype;
             }
             return *this;
         }
         virtual ~FactorGraph() {}
 
         // aliases
-        Var & var(size_t i) { return V1(i); }
-        const Var & var(size_t i) const { return V1(i); }
-        const std::vector<Var> & vars() const { return V1s(); }
-        std::vector<Var> & vars() { return V1s(); }
-        size_t nrVars() const { return V1s().size(); }
-        Factor & factor(size_t I) { return V2(I); }
-        const Factor & factor(size_t I) const { return V2(I); }
-        const std::vector<Factor> & factors() const { return V2s(); }
-        std::vector<Factor> & factors() { return V2s(); }
-        size_t nrFactors() const { return V2s().size(); }
-
-        /// Provides read access to neighbours of variable
-        const _nb_t & nbV( size_t i1 ) const { return nb1(i1); }
-        /// Provides full access to neighbours of variable
-        _nb_t & nbV( size_t i1 ) { return nb1(i1); }
-        /// Provides read access to neighbours of factor
-        const _nb_t & nbF( size_t i2 ) const { return nb2(i2); }
-        /// Provides full access to neighbours of factor
-        _nb_t & nbF( size_t i2 ) { return nb2(i2); }
+        Var & var(size_t i) { return vars[i]; }
+        const Var & var(size_t i) const { return vars[i]; }
+        Factor & factor(size_t I) { return factors[I]; }
+        const Factor & factor(size_t I) const { return factors[I]; }
+
+        size_t nrVars() const { return vars.size(); }
+        size_t nrFactors() const { return factors.size(); }
+        size_t nrEdges() const { return G.nrEdges(); }
+
+        /// Provides read access to neighbors of variable
+        const Neighbors & nbV( size_t i ) const { return G.nb1(i); }
+        /// Provides full access to neighbors of variable
+        Neighbors & nbV( size_t i ) { return G.nb1(i); }
+        /// Provides read access to neighbors of factor
+        const Neighbors & nbF( size_t I ) const { return G.nb2(I); }
+        /// Provides full access to neighbors of factor
+        Neighbors & nbF( size_t I ) { return G.nb2(I); }
+        /// Provides read access to neighbor of variable
+        const Neighbor & nbV( size_t i, size_t _I ) const { return G.nb1(i)[_I]; }
+        /// Provides full access to neighbor of variable
+        Neighbor & nbV( size_t i, size_t _I ) { return G.nb1(i)[_I]; }
+        /// Provides read access to neighbor of factor
+        const Neighbor & nbF( size_t I, size_t _i ) const { return G.nb2(I)[_i]; }
+        /// Provides full access to neighbor of factor
+        Neighbor & nbF( size_t I, size_t _i ) { return G.nb2(I)[_i]; }
 
         size_t findVar(const Var & n) const {
-            size_t i = find( vars().begin(), vars().end(), n ) - vars().begin();
+            size_t i = find( vars.begin(), vars.end(), n ) - vars.begin();
             assert( i != nrVars() );
             return i;
         }
@@ -102,18 +117,14 @@ class FactorGraph : public BipartiteGraph<Var,Factor> {
         friend std::ostream& operator << (std::ostream& os, const FactorGraph& fg);
         friend std::istream& operator >> (std::istream& is, FactorGraph& fg);
 
-        VarSet delta(const Var & n) const;
-        VarSet Delta(const Var & n) const;
-        virtual void makeFactorCavity(size_t I);
-        virtual void makeCavity(const Var & n);
+        VarSet delta( unsigned i ) const;
+        VarSet Delta( unsigned i ) const;
+        virtual void makeCavity( unsigned i );
 
         long ReadFromFile(const char *filename);
         long WriteToFile(const char *filename) const;
         long WriteToDotFile(const char *filename) const;
 
-        Factor ExactMarginal(const VarSet & x) const;
-        Real ExactlogZ() const;
-
         virtual void clamp( const Var & n, size_t i );
         
         bool hasNegatives() const;
@@ -126,31 +137,29 @@ class FactorGraph : public BipartiteGraph<Var,Factor> {
         virtual void undoProb( size_t I );
         void saveProb( size_t I );
 
-        bool isConnected() const;
-
         virtual void updatedFactor( size_t /*I*/ ) {};
 
     private:
-        /// Part of constructors (creates edges, neighbours and adjacency matrix)
+        /// Part of constructors (creates edges, neighbors and adjacency matrix)
         void createGraph( size_t nrEdges );
 };
 
 
 // assumes that the set of variables in [var_begin,var_end) is the union of the variables in the factors in [fact_begin, fact_end)
 template<typename FactorInputIterator, typename VarInputIterator>
-FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint, size_t nr_var_hint ) : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {
+FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint, size_t nr_var_hint ) : G(), _undoProbs(), _normtype(Prob::NORMPROB) {
     // add factors
     size_t nrEdges = 0;
-    V2s().reserve( nr_fact_hint );
+    factors.reserve( nr_fact_hint );
     for( FactorInputIterator p2 = fact_begin; p2 != fact_end; ++p2 ) {
-        V2s().push_back( *p2 );
-       nrEdges += p2->vars().size();
+        factors.push_back( *p2 );
+        nrEdges += p2->vars().size();
     }
  
     // add variables
-    V1s().reserve( nr_var_hint );
+    vars.reserve( nr_var_hint );
     for( VarInputIterator p1 = var_begin; p1 != var_end; ++p1 )
-       V1s().push_back( *p1 );
+        vars.push_back( *p1 );
 
     // create graph structure
     createGraph( nrEdges );
index fac40ba..f1adfc4 100644 (file)
@@ -35,10 +35,10 @@ namespace dai {
 /// HAK provides an implementation of the single and double-loop algorithms by Heskes, Albers and Kappen
 class HAK : public DAIAlgRG {
     protected:
-        std::vector<Factor>          _Qa;
-        std::vector<Factor>          _Qb;
-        std::vector<Factor>          _muab;
-        std::vector<Factor>          _muba;
+        std::vector<Factor>                _Qa;
+        std::vector<Factor>                _Qb;
+        std::vector<std::vector<Factor> >  _muab;
+        std::vector<std::vector<Factor> >  _muba;
         
     public:
         /// Default constructor
@@ -75,12 +75,11 @@ class HAK : public DAIAlgRG {
         bool DoubleLoop() { return GetPropertyAs<bool>("doubleloop"); }
         size_t LoopDepth() { return GetPropertyAs<size_t>("loopdepth"); }
 
-        Factor & muab( size_t alpha, size_t beta ) { return _muab[ORIR2E(alpha,beta)]; }
-        Factor & muba( size_t beta, size_t alpha ) { return _muba[ORIR2E(alpha,beta)]; }
+        Factor & muab( size_t alpha, size_t _beta ) { return _muab[alpha][_beta]; }
+        Factor & muba( size_t alpha, size_t _beta ) { return _muba[alpha][_beta]; }
         const Factor& Qa( size_t alpha ) const { return _Qa[alpha]; };
         const Factor& Qb( size_t beta ) const { return _Qb[beta]; };
 
-//      void Regenerate();
         double doGBP();
         double doDoubleLoop();
         double run();
index 48b8b47..e6bf12f 100644 (file)
@@ -42,7 +42,7 @@ class JTree : public DAIAlgRG {
         DEdgeVec             _RTree;     // rooted tree
         std::vector<Factor>  _Qa;
         std::vector<Factor>  _Qb;
-        std::vector<Factor>  _mes;
+        std::vector<std::vector<Factor> >  _mes;
         double               _logZ;
 
 
@@ -67,15 +67,12 @@ class JTree : public DAIAlgRG {
         JTree( const FactorGraph &fg, const Properties &opts, bool automatic=true );
         void GenerateJT( const std::vector<VarSet> &Cliques );
 
-        Factor & message(size_t i1, size_t i2) { return( _mes[ORIR2E(i1,i2)] ); }   
-        const Factor & message(size_t i1, size_t i2) const { return( _mes[ORIR2E(i1,i2)] ); }   
+        Factor & message( size_t alpha, size_t _beta ) { return _mes[alpha][_beta]; }   
+        const Factor & message( size_t alpha, size_t _beta ) const { return _mes[alpha][_beta]; }   
 
         static const char *Name;
         std::string identify() const;
-//      void Regenerate();
-        void init() {
-            assert( checkProperties() );
-        }
+        void init() { assert( checkProperties() ); }
         void runHUGIN();
         void runShaferShenoy();
         double run();
index b7d5fee..231b043 100644 (file)
@@ -34,19 +34,14 @@ namespace dai {
 
 class LC : public DAIAlgFG {
     protected:
-        typedef struct { size_t i; size_t I; } _iI_type;
-
         std::vector<Factor>      _pancakes;      // used by all LC types (psi_I is stored in the pancake)
         std::vector<Factor>      _cavitydists;   // used by all LC types to store the approximate cavity distribution
-        /// _phis[VV2E(i,I)] corresponds to \f$ \phi^{\setminus i}_I(x_{I \setminus i}) \f$
-        std::vector<Factor>      _phis;
+        /// _phis[i][_I] corresponds to \f$ \phi^{\setminus i}_I(x_{I \setminus i}) \f$
+        std::vector<std::vector<Factor> >      _phis;
 
         /// Single variable beliefs
         std::vector<Factor>      _beliefs;
 
-        /// For each pair (i,j) with j in delta(i), store i and the common factor I
-        std::vector<_iI_type>    _iI;
-
     public:
         ENUM6(CavityType,FULL,PAIR,PAIR2,PAIRINT,PAIRCUM,UNIFORM)
         ENUM3(UpdateType,SEQFIX,SEQRND,NONE)
@@ -58,7 +53,7 @@ class LC : public DAIAlgFG {
         /// Default constructor
         LC() : DAIAlgFG() {};
         /// Copy constructor
-        LC(const LC & x) : DAIAlgFG(x), _pancakes(x._pancakes), _cavitydists(x._cavitydists), _phis(x._phis), _beliefs(x._beliefs), _iI(x._iI) {};
+        LC(const LC & x) : DAIAlgFG(x), _pancakes(x._pancakes), _cavitydists(x._cavitydists), _phis(x._phis), _beliefs(x._beliefs) {};
         /// Clone function
         LC* clone() const { return new LC(*this); }
         /// Construct LC object from a FactorGraph and parameters
@@ -71,7 +66,6 @@ class LC : public DAIAlgFG {
                 _cavitydists    = x._cavitydists;
                 _phis           = x._phis;
                 _beliefs        = x._beliefs;
-                _iI             = x._iI;
             }
             return *this;
         }
@@ -82,7 +76,7 @@ class LC : public DAIAlgFG {
         long SetCavityDists( std::vector<Factor> &Q );
 
         void init();
-        Factor NewPancake (size_t iI, bool & hasNaNs);
+        Factor NewPancake (size_t i, size_t _I, bool & hasNaNs);
         double run();
 
         std::string identify() const;
@@ -94,12 +88,10 @@ class LC : public DAIAlgFG {
         const Factor & belief (size_t i) const { return _beliefs[i]; };
         const Factor & pancake (size_t i) const { return _pancakes[i]; };
         const Factor & cavitydist (size_t i) const { return _cavitydists[i]; };
-        size_t nr_iI() const { return _iI.size(); };
 
         void clamp( const Var &/*n*/, size_t /*i*/ ) { assert( 0 == 1 ); }
         void undoProbs( const VarSet &/*ns*/ ) { assert( 0 == 1 ); }
         void saveProbs( const VarSet &/*ns*/ ) { assert( 0 == 1 ); }
-        void makeFactorCavity(size_t /*I*/) { assert( 0 == 1 ); }
         virtual void makeCavity(const Var & /*n*/) { assert( 0 == 1 ); }
         bool checkProperties();
 };
index d44cf0e..29aa907 100644 (file)
@@ -60,7 +60,7 @@ class MF : public DAIAlgFG {
         void Regenerate();
         void init();
         double run();
-        Factor belief1 (size_t i) const;
+        Factor beliefV (size_t i) const;
         Factor belief (const Var &n) const;
         Factor belief (const VarSet &ns) const;
         std::vector<Factor> beliefs() const;
index b7ab9e7..09e17a5 100644 (file)
@@ -96,98 +96,73 @@ class FRegion : public Factor {
 };
 
 
-typedef BipartiteGraph<FRegion,Region> BipRegGraph;
-
-
 /// A RegionGraph is a bipartite graph consisting of outer regions (type FRegion) and inner regions (type Region)
-class RegionGraph : public FactorGraph, BipRegGraph {
+class RegionGraph : public FactorGraph {
     public:
-        typedef BipRegGraph::_nb_t      R_nb_t;
-        typedef R_nb_t::const_iterator  R_nb_cit;
-        typedef BipRegGraph::_edge_t    R_edge_t;
+        BipartiteGraph          G;
+        std::vector<FRegion>    ORs;
+        std::vector<Region>     IRs;
 
-        
     protected:
         /// Give back the OR index that corresponds to a factor index
-        typedef std::map<size_t,size_t>::const_iterator fac2OR_cit;
-        std::map<size_t,size_t>         _fac2OR;
+        std::vector<size_t>     fac2OR;
 
 
     public:
         /// Default constructor
-        RegionGraph() : FactorGraph(), BipRegGraph(), _fac2OR() {}
+        RegionGraph() : FactorGraph(), G(), ORs(), IRs(), fac2OR() {}
 
         /// Constructs a RegionGraph from a FactorGraph
-        RegionGraph(const FactorGraph & fg) : FactorGraph(fg), BipRegGraph(), _fac2OR() {}
+        RegionGraph( const FactorGraph &fg ) : FactorGraph(fg), G(), ORs(), IRs(), fac2OR() {}
 
         /// Constructs a RegionGraph from a FactorGraph, a vector of outer regions, a vector of inner regions and a vector of edges
-        RegionGraph(const FactorGraph & fg, const std::vector<Region> & ors, const std::vector<Region> & irs, const std::vector<R_edge_t> & edges);
+        RegionGraph( const FactorGraph &fg, const std::vector<Region> &ors, const std::vector<Region> &irs, const std::vector<std::pair<size_t,size_t> > &edges );
         
         /// Constructs a RegionGraph from a FactorGraph and a vector of outer VarSets (CVM style)
-        RegionGraph(const FactorGraph & fg, const std::vector<VarSet> & cl);
+        RegionGraph( const FactorGraph &fg, const std::vector<VarSet> &cl );
 
         /// Copy constructor
-        RegionGraph(const RegionGraph & x) : FactorGraph(x), BipRegGraph(x), _fac2OR(x._fac2OR) {}
+        RegionGraph( const RegionGraph &x ) : FactorGraph(x), G(x.G), ORs(x.ORs), IRs(x.IRs), fac2OR(x.fac2OR) {}
 
         /// Assignment operator
-        RegionGraph & operator=(const RegionGraph & x) {
+        RegionGraph & operator=( const RegionGraph &x ) {
             if( this != &x ) {
-                FactorGraph::operator=(x);
-                BipRegGraph::operator=(x);
-                _fac2OR = x._fac2OR;
+                FactorGraph::operator=( x );
+                G = x.G;
+                ORs = x.ORs;
+                IRs = x.IRs;
+                fac2OR = x.fac2OR;
             }
             return *this;
         }
 
         /// Provides read access to outer region
-        const FRegion & OR(long alpha) const { return BipRegGraph::V1(alpha); }
+        const FRegion & OR(size_t alpha) const { return ORs[alpha]; }
         /// Provides access to outer region
-        FRegion & OR(long alpha) { return BipRegGraph::V1(alpha); }
-        /// Provides read access to all outer regions
-        const std::vector<FRegion> & ORs() const { return BipRegGraph::V1s(); }
-        /// Provides access to all outer regions
-        std::vector<FRegion> &ORs() { return BipRegGraph::V1s(); }
-        /// Returns number of outer regions
-        size_t nr_ORs() const { return BipRegGraph::V1s().size(); }
+        FRegion & OR(size_t alpha) { return ORs[alpha]; }
 
         /// Provides read access to inner region
-        const Region & IR(long beta) const { return BipRegGraph::V2(beta); }
+        const Region & IR(size_t beta) const { return IRs[beta]; }
         /// Provides access to inner region
-        Region & IR(long beta) { return BipRegGraph::V2(beta); }
-        /// Provides read access to all inner regions
-        const std::vector<Region> & IRs() const { return BipRegGraph::V2s(); }
-        /// Provides access to all inner regions
-        std::vector<Region> & IRs() { return BipRegGraph::V2s(); }
+        Region & IR(size_t beta) { return IRs[beta]; }
+
+        /// Returns number of outer regions
+        size_t nrORs() const { return ORs.size(); }
         /// Returns number of inner regions
-        size_t nr_IRs() const { return BipRegGraph::V2s().size(); }
-
-        /// Provides read access to edge
-        const R_edge_t & Redge(size_t ind) const { return BipRegGraph::edge(ind); }
-        /// Provides full access to edge
-        R_edge_t & Redge(size_t ind) { return BipRegGraph::edge(ind); }
-        /// Provides read access to all edges
-        const std::vector<R_edge_t> & Redges() const { return BipRegGraph::edges(); }
-        /// Provides full access to all edges
-        std::vector<R_edge_t> & Redges() { return BipRegGraph::edges(); }
-        /// Returns number of edges
-        size_t nr_Redges() const { return BipRegGraph::edges().size(); }
-
-        /// Provides read access to neighbours of outer region
-        const R_nb_t & nbOR( size_t i1 ) const { return BipRegGraph::nb1(i1); }
-        /// Provides full access to neighbours of outer region
-        R_nb_t & nbOR( size_t i1 ) { return BipRegGraph::nb1(i1); }
-
-        /// Provides read access to neighbours of inner region
-        const R_nb_t & nbIR( size_t i2 ) const { return BipRegGraph::nb2(i2); }
-        /// Provides full access to neighbours of inner region
-        R_nb_t & nbIR( size_t i2 ) { return BipRegGraph::nb2(i2); }
-
-        /// Converts the pair of outer/inner region indices (i1,i2) to the corresponding edge index
-        size_t ORIR2E( const size_t i1, const size_t i2 ) const { return BipRegGraph::VV2E(i1, i2); }
-
-        void Regenerate() { BipRegGraph::Regenerate(); }
-        
-        
+        size_t nrIRs() const { return IRs.size(); }
+
+
+        /// Provides read access to neighbors of outer region
+        const Neighbors & nbOR( size_t alpha ) const { return G.nb1(alpha); }
+        /// Provides full access to neighbors of outer region
+        Neighbors & nbOR( size_t alpha ) { return G.nb1(alpha); }
+
+        /// Provides read access to neighbors of inner region
+        const Neighbors & nbIR( size_t beta ) const { return G.nb2(beta); }
+        /// Provides full access to neighbors of inner region
+        Neighbors & nbIR( size_t beta ) { return G.nb2(beta); }
+
+
         /// Calculates counting numbers of inner regions based upon counting numbers of outer regions
         void Calc_Counting_Numbers();
         /// Check whether the counting numbers are valid
@@ -206,10 +181,7 @@ class RegionGraph : public FactorGraph, BipRegGraph {
         void clamp( const Var &n, size_t i ) { FactorGraph::clamp( n, i ); RecomputeORs( n ); }
 
         /// We have to overload FactorGraph::makeCavity because the corresponding outer regions have to be recomputed
-        void makeCavity( const Var &n ) { FactorGraph::makeCavity( n ); RecomputeORs( n ); }
-
-        /// We have to overload FactorGraph::makeFactorCavity because the corresponding outer regions have to be recomputed
-        void makeFactorCavity( size_t I ) { FactorGraph::makeFactorCavity( I ); RecomputeOR( I ); }
+        void makeCavity( size_t i ) { FactorGraph::makeCavity( i ); RecomputeORs( var(i) ); }
 
         /// We have to overload FactorGraph::undoProbs because the corresponding outer regions have to be recomputed
         void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs( ns ); RecomputeORs( ns ); }
index 900e906..4d3e417 100644 (file)
@@ -43,8 +43,8 @@ class TreeEPSubTree {
         std::vector<Factor>  _Qa;
         std::vector<Factor>  _Qb;
         DEdgeVec             _RTree;
-        std::vector<size_t>  _a;             // _Qa[alpha]  <->  superTree._Qa[_a[alpha]]
-        std::vector<size_t>  _b;             // _Qb[beta]   <->  superTree._Qb[_b[beta]]
+        std::vector<size_t>  _a;        // _Qa[alpha]  <->  superTree._Qa[_a[alpha]]
+        std::vector<size_t>  _b;        // _Qb[beta]   <->  superTree._Qb[_b[beta]]
                                         // _Qb[beta]   <->  _RTree[beta]    
         const Factor *       _I;
         VarSet               _ns;
@@ -113,7 +113,7 @@ class TreeEP : public JTree {
         double run();
         Complex logZ() const;
 
-        bool offtree(size_t I) const { return !_fac2OR.count(I); }
+        bool offtree( size_t I ) const { return (fac2OR[I] == -1U); }
 
         void init( const VarSet &/*ns*/ ) { init(); }
         void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
index 9356e8c..7dfea58 100644 (file)
 
 #include <sys/times.h>
 #include <cstdio>
+#include <boost/foreach.hpp>
+
+
+#define foreach BOOST_FOREACH
 
 
 namespace dai {
index 872a7b5..b8b8b24 100644 (file)
@@ -30,20 +30,34 @@ using namespace std;
 
 
 InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const Properties &opts ) {
+#ifdef WITH_BP
     if( name == BP::Name ) 
         return new BP (fg, opts);
+#endif
+#ifdef WITH_MF
     else if( name == MF::Name ) 
         return new MF (fg, opts);
+#endif
+#ifdef WITH_HAK
     else if( name == HAK::Name ) 
         return new HAK (fg, opts);
+#endif
+#ifdef WITH_LC
     else if( name == LC::Name )
         return new LC (fg, opts);
+#endif
+#ifdef WITH_TREEEP
     else if( name == TreeEP::Name )
         return new TreeEP (fg, opts);
-    else if( name == MR::Name )
-        return new MR (fg, opts);
+#endif
+#ifdef WITH_JTREE
     else if( name == JTree::Name )
         return new JTree (fg, opts);
+#endif
+#ifdef WITH_MR
+    else if( name == MR::Name )
+        return new MR (fg, opts);
+#endif
     else
         throw "Unknown inference algorithm";
 }
index 3b7477b..0d65992 100644 (file)
@@ -59,44 +59,56 @@ bool BP::checkProperties() {
 
 
 void BP::Regenerate() {
-    DAIAlgFG::Regenerate();
-    
-    // clear messages
-    _messages.clear();
-    _messages.reserve(nr_edges());
-
-    // clear indices
-    _indices.clear();
-    _indices.reserve(nr_edges());
-
-    // create messages and indices
-    for( vector<_edge_t>::const_iterator iI=edges().begin(); iI!=edges().end(); ++iI ) {
-        _messages.push_back( Prob( var(iI->first).states() ) );
-
-        vector<size_t> ind( factor(iI->second).stateSpace(), 0 );
-        Index i (var(iI->first), factor(iI->second).vars() );
-        for( size_t j = 0; i >= 0; ++i,++j )
-            ind[j] = i; 
-        _indices.push_back( ind );
+    // create edge properties
+    edges.clear();
+    edges.reserve( nrVars() );
+    for( size_t i = 0; i < nrVars(); ++i ) {
+        edges.push_back( vector<EdgeProp>() );
+        edges[i].reserve( nbV(i).size() ); 
+        foreach( const Neighbor &I, nbV(i) ) {
+            EdgeProp newEP;
+            newEP.message = Prob( var(i).states() );
+            newEP.newMessage = Prob( var(i).states() );
+
+            newEP.index.reserve( factor(I).stateSpace() );
+            for( Index k( var(i), factor(I).vars() ); k >= 0; ++k )
+                newEP.index.push_back( k );
+
+            newEP.residual = 0.0;
+            edges[i].push_back( newEP );
+        }
     }
-
-    // create new_messages
-    _newmessages = _messages;
 }
 
 
 void BP::init() {
     assert( checkProperties() );
-    for( vector<Prob>::iterator mij = _messages.begin(); mij != _messages.end(); ++mij )
-        mij->fill(1.0 / mij->size());
-    _newmessages = _messages;
+    for( size_t i = 0; i < nrVars(); ++i ) {
+        foreach( const Neighbor &I, nbV(i) ) {
+            message( i, I.iter ).fill( 1.0 );
+            newMessage( i, I.iter ).fill( 1.0 );
+        }
+    }
+}
+
+
+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 );
+            }
 }
 
 
-void BP::calcNewMessage (size_t iI) {
+void BP::calcNewMessage( size_t i, size_t _I ) {
     // calculate updated message I->i
-    size_t i = edge(iI).first;
-    size_t I = edge(iI).second;
+    size_t I = nbV(i,_I);
 
 /*  UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
 
@@ -112,32 +124,35 @@ void BP::calcNewMessage (size_t iI) {
     Prob prod( factor(I).p() );
 
     // Calculate product of incoming messages and factor I
-    for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); ++j )
-        if( *j != i ) {     // for all j in I \ i
+    foreach( const Neighbor &j, nbF(I) ) {
+        if( j != i ) {     // for all j in I \ i
+            size_t _I = j.dual;
             // ind is the precalculated Index(j,I) i.e. to x_I == k corresponds x_j == ind[k]
-            _ind_t* ind = &(index(*j,I));
+            const ind_t & ind = index(j, _I);
 
             // prod_j will be the product of messages coming into j
-            Prob prod_j( var(*j).states() ); 
-            for( _nb_cit J = nb1(*j).begin(); J != nb1(*j).end(); ++J )
-                if( *J != I )   // for all J in nb(j) \ I 
-                    prod_j *= message(*j,*J);
+            Prob prod_j( var(j).states() ); 
+            foreach( const Neighbor &J, nbV(j) ) {
+                if( J != I )   // for all J in nb(j) \ I 
+                    prod_j *= message( j, J.iter );
+            }
 
             // multiply prod with prod_j
             for( size_t r = 0; r < prod.size(); ++r )
-                prod[r] *= prod_j[(*ind)[r]];
+                prod[r] *= prod_j[ind[r]];
         }
+    }
 
     // Marginalize onto i
     Prob marg( var(i).states(), 0.0 );
     // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
-    _ind_t* ind = &(index(i,I));
+    const ind_t ind = index(i,_I);
     for( size_t r = 0; r < prod.size(); ++r )
-        marg[(*ind)[r]] += prod[r];
+        marg[ind[r]] += prod[r];
     marg.normalize( _normtype );
     
     // Store result
-    _newmessages[iI] = marg;
+    newMessage(i,_I) = marg;
 }
 
 
@@ -152,29 +167,30 @@ double BP::run() {
     clock_t tic = toc();
     Diffs diffs(nrVars(), 1.0);
     
-    vector<size_t> edge_seq;
-    vector<double> residuals;
+    typedef pair<size_t,size_t> Edge;
+    vector<Edge> update_seq;
 
     vector<Factor> old_beliefs;
     old_beliefs.reserve( nrVars() );
     for( size_t i = 0; i < nrVars(); ++i )
-        old_beliefs.push_back(belief1(i));
+        old_beliefs.push_back( beliefV(i) );
 
     size_t iter = 0;
+    size_t nredges = nrEdges();
 
     if( Updates() == UpdateType::SEQMAX ) {
         // do the first pass
-        for(size_t iI = 0; iI < nr_edges(); ++iI ) 
-            calcNewMessage(iI);
-
-        // calculate initial residuals
-        residuals.reserve(nr_edges());
-        for( size_t iI = 0; iI < nr_edges(); ++iI )
-            residuals.push_back( dist( _newmessages[iI], _messages[iI], Prob::DISTLINF ) );
+        for( size_t i = 0; i < nrVars(); ++i )
+            foreach( const Neighbor &I, nbV(i) ) {
+                calcNewMessage( i, I.iter );
+                // calculate initial residuals
+                residual( i, I.iter ) = dist( newMessage( i, I.iter ), message( i, I.iter ), Prob::DISTLINF );
+            }
     } else {
-        edge_seq.reserve( nr_edges() );
-        for( size_t i = 0; i < nr_edges(); ++i )
-            edge_seq.push_back( i );
+        update_seq.reserve( nredges );
+        for( size_t i = 0; i < nrVars(); ++i )
+            foreach( const Neighbor &I, nbV(i) )
+                update_seq.push_back( Edge( i, I.iter ) );
     }
 
     // do several passes over the network until maximum number of iterations has
@@ -182,47 +198,51 @@ double BP::run() {
     for( iter=0; iter < MaxIter() && diffs.max() > Tol(); ++iter ) {
         if( Updates() == UpdateType::SEQMAX ) {
             // Residuals-BP by Koller et al.
-            for( size_t t = 0; t < nr_edges(); ++t ) {
+            for( size_t t = 0; t < nredges; ++t ) {
                 // update the message with the largest residual
-                size_t iI = max_element(residuals.begin(), residuals.end()) - residuals.begin();
-                _messages[iI] = _newmessages[iI];
-                residuals[iI] = 0;
+
+                size_t i, _I;
+                findMaxResidual( i, _I );
+                message( i, _I ) = newMessage( i, _I );
+                residual( i, _I ) = 0.0;
 
                 // I->i has been updated, which means that residuals for all
                 // J->j with J in nb[i]\I and j in nb[J]\i have to be updated
-                size_t i = edge(iI).first;
-                size_t I = edge(iI).second;
-                for( _nb_cit J = nb1(i).begin(); J != nb1(i).end(); ++J ) 
-                    if( *J != I )
-                        for( _nb_cit j = nb2(*J).begin(); j != nb2(*J).end(); ++j )
-                            if( *j != i ) {
-                                size_t jJ = VV2E(*j,*J);
-                                calcNewMessage(jJ);
-                                residuals[jJ] = dist( _newmessages[jJ], _messages[jJ], Prob::DISTLINF );
+                foreach( const Neighbor &J, nbV(i) ) {
+                    if( J.iter != _I ) {
+                        foreach( const Neighbor &j, nbF(J) ) {
+                            size_t _J = j.dual;
+                            if( j != i ) {
+                                calcNewMessage( j, _J );
+                                residual( j, _J ) = dist( newMessage( j, _J ), message( j, _J ), Prob::DISTLINF );
                             }
+                        }
+                    }
+                }
             }
         } else if( Updates() == UpdateType::PARALL ) {
             // Parallel updates 
-            for( size_t t = 0; t < nr_edges(); ++t )
-                calcNewMessage(t);
+            for( size_t i = 0; i < nrVars(); ++i )
+                foreach( const Neighbor &I, nbV(i) )
+                    calcNewMessage( i, I.iter );
 
-            for( size_t t = 0; t < nr_edges(); ++t )
-                _messages[t] = _newmessages[t];
+            for( size_t i = 0; i < nrVars(); ++i )
+                foreach( const Neighbor &I, nbV(i) )
+                    message( i, I.iter ) = newMessage( i, I.iter );
         } else {
             // Sequential updates
             if( Updates() == UpdateType::SEQRND )
-                random_shuffle( edge_seq.begin(), edge_seq.end() );
+                random_shuffle( update_seq.begin(), update_seq.end() );
             
-            for( size_t t = 0; t < nr_edges(); ++t ) {
-                size_t k = edge_seq[t];
-                calcNewMessage(k);
-                _messages[k] = _newmessages[k];
+            foreach( const Edge &e, update_seq ) {
+                calcNewMessage( e.first, e.second );
+                message( e.first, e.second ) = newMessage( e.first, e.second );
             }
         }
 
         // calculate new beliefs and compare with old ones
         for( size_t i = 0; i < nrVars(); ++i ) {
-            Factor nb( belief1(i) );
+            Factor nb( beliefV(i) );
             diffs.push( dist( nb, old_beliefs[i], Prob::DISTLINF ) );
             old_beliefs[i] = nb;
         }
@@ -249,10 +269,10 @@ double BP::run() {
 }
 
 
-Factor BP::belief1( size_t i ) const {
+Factor BP::beliefV( size_t i ) const {
     Prob prod( var(i).states() ); 
-    for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); ++I ) 
-        prod *= newMessage(i,*I);
+    foreach( const Neighbor &I, nbV(i) )
+        prod *= newMessage( i, I.iter );
 
     prod.normalize( Prob::NORMPROB );
     return( Factor( var(i), prod ) );
@@ -260,16 +280,16 @@ Factor BP::belief1( size_t i ) const {
 
 
 Factor BP::belief (const Var &n) const {
-    return( belief1( findVar( n ) ) );
+    return( beliefV( findVar( n ) ) );
 }
 
 
 vector<Factor> BP::beliefs() const {
     vector<Factor> result;
     for( size_t i = 0; i < nrVars(); ++i )
-        result.push_back( belief1(i) );
+        result.push_back( beliefV(i) );
     for( size_t I = 0; I < nrFactors(); ++I )
-        result.push_back( belief2(I) );
+        result.push_back( beliefF(I) );
     return result;
 }
 
@@ -283,27 +303,29 @@ Factor BP::belief( const VarSet &ns ) const {
             if( factor(I).vars() >> ns )
                 break;
         assert( I != nrFactors() );
-        return belief2(I).marginal(ns);
+        return beliefF(I).marginal(ns);
     }
 }
 
 
-Factor BP::belief2 (size_t I) const {
+Factor BP::beliefF (size_t I) const {
     Prob prod( factor(I).p() );
 
-    for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); ++j ) {
+    foreach( const Neighbor &j, nbF(I) ) {
+        size_t _I = j.dual;
         // ind is the precalculated Index(j,I) i.e. to x_I == k corresponds x_j == ind[k]
-        const _ind_t *ind = &(index(*j, I));
+        const ind_t & ind = index(j, _I);
 
         // prod_j will be the product of messages coming into j
-        Prob prod_j( var(*j).states() ); 
-        for( _nb_cit J = nb1(*j).begin(); J != nb1(*j).end(); ++J ) 
-            if( *J != I )   // for all J in nb(j) \ I 
-                prod_j *= newMessage(*j,*J);
+        Prob prod_j( var(j).states() ); 
+        foreach( const Neighbor &J, nbV(j) ) {
+            if( J != I )   // for all J in nb(j) \ I 
+                prod_j *= newMessage( j, J.iter );
+        }
 
         // multiply prod with prod_j
         for( size_t r = 0; r < prod.size(); ++r )
-            prod[r] *= prod_j[(*ind)[r]];
+            prod[r] *= prod_j[ind[r]];
     }
 
     Factor result( factor(I).vars(), prod );
@@ -326,9 +348,9 @@ Factor BP::belief2 (size_t I) const {
 Complex BP::logZ() const {
     Complex sum = 0.0;
     for(size_t i = 0; i < nrVars(); ++i )
-        sum += Complex(1.0 - nb1(i).size()) * belief1(i).entropy();
+        sum += Complex(1.0 - nbV(i).size()) * beliefV(i).entropy();
     for( size_t I = 0; I < nrFactors(); ++I )
-        sum -= KL_dist( belief2(I), factor(I) );
+        sum -= KL_dist( beliefF(I), factor(I) );
     return sum;
 }
 
@@ -343,8 +365,8 @@ string BP::identify() const {
 void BP::init( const VarSet &ns ) {
     for( VarSet::const_iterator n = ns.begin(); n != ns.end(); ++n ) {
         size_t ni = findVar( *n );
-        for( _nb_cit I = nb1(ni).begin(); I != nb1(ni).end(); ++I )
-            message(ni,*I).fill( 1.0 );
+        foreach( const Neighbor &I, nbV( ni ) )
+            message( ni, I.iter ).fill( 1.0 );
     }
 }
 
index fcbf626..1d5f1a2 100644 (file)
@@ -29,6 +29,7 @@
 #include <functional>
 #include <tr1/unordered_map>
 #include <dai/factorgraph.h>
+#include <dai/util.h>
 
 
 namespace dai {
@@ -37,21 +38,21 @@ namespace dai {
 using namespace std;
 
 
-FactorGraph::FactorGraph( const vector<Factor> &P ) : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {
+FactorGraph::FactorGraph( const vector<Factor> &P ) : G(), _undoProbs(), _normtype(Prob::NORMPROB) {
     // add factors, obtain variables
     set<Var> _vars;
-    V2s().reserve( P.size() );
+    factors.reserve( P.size() );
     size_t nrEdges = 0;
     for( vector<Factor>::const_iterator p2 = P.begin(); p2 != P.end(); p2++ ) {
-        V2s().push_back( *p2 );
-       copy( p2->vars().begin(), p2->vars().end(), inserter( _vars, _vars.begin() ) );
-       nrEdges += p2->vars().size();
+        factors.push_back( *p2 );
+        copy( p2->vars().begin(), p2->vars().end(), inserter( _vars, _vars.begin() ) );
+        nrEdges += p2->vars().size();
     }
 
     // add _vars
-    V1s().reserve( _vars.size() );
+    vars.reserve( _vars.size() );
     for( set<Var>::const_iterator p1 = _vars.begin(); p1 != _vars.end(); p1++ )
-        V1s().push_back( *p1 );
+        vars.push_back( *p1 );
     
     // create graph structure
     createGraph( nrEdges );
@@ -63,39 +64,24 @@ void FactorGraph::createGraph( size_t nrEdges ) {
     // create a mapping for indices
     std::tr1::unordered_map<size_t, size_t> hashmap;
     
-    for( size_t i = 0; i < vars().size(); i++ )        
-       hashmap[vars()[i].label()] = i;
+    for( size_t i = 0; i < vars.size(); i++ )
+        hashmap[var(i).label()] = i;
     
-    // create edges
-    edges().reserve( nrEdges );
+    // create edge list
+    typedef pair<unsigned,unsigned> Edge;
+    vector<Edge> edges;
+    edges.reserve( nrEdges );
     for( size_t i2 = 0; i2 < nrFactors(); i2++ ) {
         const VarSet& ns = factor(i2).vars();
         for( VarSet::const_iterator q = ns.begin(); q != ns.end(); q++ )
-            edges().push_back(_edge_t(hashmap[q->label()], i2));               
+            edges.push_back( Edge(hashmap[q->label()], i2) );
     }
 
-    // calc neighbours and adjacency matrix
-    Regenerate();
+    // create bipartite graph
+    G.create( nrVars(), nrFactors(), edges.begin(), edges.end() );
 }
 
 
-/*FactorGraph& FactorGraph::addFactor( const Factor &I ) {
-    // add Factor
-    _V2.push_back( I );
-
-    // add new vars in Factor
-    for( VarSet::const_iterator i = I.vars().begin(); i != I.vars().end(); i++ ) {
-        size_t i_ind = find(vars().begin(), vars().end(), *i) - vars().begin();
-        if( i_ind == vars().size() )
-            _V1.push_back( *i );
-        _E12.push_back( _edge_t( i_ind, nrFactors() - 1 ) );
-    }
-
-    Regenerate();
-    return(*this);
-}*/
-
-
 ostream& operator << (ostream& os, const FactorGraph& fg) {
     os << fg.nrFactors() << endl;
 
@@ -268,37 +254,27 @@ istream& operator >> (istream& is, FactorGraph& fg) {
 }
 
 
-VarSet FactorGraph::delta(const Var & n) const {
+VarSet FactorGraph::delta( unsigned i ) const {
     // calculate Markov Blanket
-    size_t i = findVar( n );
-
     VarSet del;
-    for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); ++I )
-        for( _nb_cit j = nb2(*I).begin(); j != nb2(*I).end(); ++j )
-            if( *j != i )
-                del |= var(*j);
+    foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        foreach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
+            if( j != i )
+                del |= var(j);
 
     return del;
 }
 
 
-VarSet FactorGraph::Delta(const Var & n) const {
-    return( delta(n) | n );
+VarSet FactorGraph::Delta( unsigned i ) const {
+    return( delta(i) | var(i) );
 }
 
 
-void FactorGraph::makeFactorCavity(size_t I) {
-    // fill Factor I with ones
-    factor(I).fill(1.0);
-}
-
-
-void FactorGraph::makeCavity(const Var & n) {
-    // fills all Factors that include Var n with ones
-    size_t i = findVar( n );
-
-    for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); ++I )
-        factor(*I).fill(1.0);
+void FactorGraph::makeCavity( unsigned i ) {
+    // fills all Factors that include var(i) with ones
+    foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        factor(I).fill( 1.0 );
 }
 
 
@@ -306,62 +282,11 @@ bool FactorGraph::hasNegatives() const {
     bool result = false;
     for( size_t I = 0; I < nrFactors() && !result; I++ )
         if( factor(I).hasNegatives() )
-           result = true;
+            result = true;
     return result;
 }
  
 
-/*FactorGraph & FactorGraph::DeleteFactor(size_t I) {
-    // Go through all edges
-    for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
-        if( edge->second >= I ) {
-            if( edge->second == I )
-                edge->second = -1UL;
-            else 
-                (edge->second)--;
-        }
-    // Remove all edges containing I
-    for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
-        if( edge->second == -1UL )
-            edge = _E12.erase( edge );
-//  vector<_edge_t>::iterator new_end = _E12.remove_if( _E12.begin(), _E12.end(), compose1( bind2nd(equal_to<size_t>(), -1), select2nd<_edge_t>() ) );
-//  _E12.erase( new_end, _E12.end() );
-
-    // Erase the factor
-    _V2.erase( _V2.begin() + I );
-    
-    Regenerate();
-
-    return *this;
-}
-
-
-FactorGraph & FactorGraph::DeleteVar(size_t i) {
-    // Go through all edges
-    for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
-        if( edge->first >= i ) {
-            if( edge->first == i )
-                edge->first = -1UL;
-            else 
-                (edge->first)--;
-        }
-    // Remove all edges containing i
-    for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
-        if( edge->first == -1UL )
-            edge = _E12.erase( edge );
-                
-//  vector<_edge_t>::iterator new_end = _E12.remove_if( _E12.begin(), _E12.end(), compose1( bind2nd(equal_to<size_t>(), -1), select1st<_edge_t>() ) );
-//  _E12.erase( new_end, _E12.end() );
-
-    // Erase the variable
-    _V1.erase( _V1.begin() + i );
-    
-    Regenerate();
-
-    return *this;
-}*/
-
-
 long FactorGraph::ReadFromFile(const char *filename) {
     ifstream infile;
     infile.open (filename);
@@ -408,8 +333,9 @@ long FactorGraph::WriteToDotFile(const char *filename) const {
             outfile << "node[shape=box,style=filled,color=lightgrey,width=0.3,height=0.3,fixedsize=true];" << endl;
             for( size_t I = 0; I < nrFactors(); I++ )
                 outfile << "\tp" << I << ";" << endl;
-            for( size_t iI = 0; iI < nr_edges(); iI++ )
-                outfile << "\tx" << var(edge(iI).first).label() << " -- p" << edge(iI).second << ";" << endl;
+            for( size_t i = 0; i < nrVars(); i++ )
+                foreach( const Neighbor &I, nbV(i) )  // for all neighboring factors I of i
+                    outfile << "\tx" << var(i).label() << " -- p" << I << ";" << endl;
             outfile << "}" << endl;
         } catch (char *e) {
             cout << e << endl;
@@ -467,22 +393,6 @@ void RemoveShortLoops(vector<Factor> &P) {
 }
 
 
-Factor FactorGraph::ExactMarginal(const VarSet & x) const {
-    Factor P;
-    for( size_t I = 0; I < nrFactors(); I++ )
-        P *= factor(I);
-    return P.marginal(x);
-}
-
-
-Real FactorGraph::ExactlogZ() const {
-    Factor P;
-    for( size_t I = 0; I < nrFactors(); I++ )
-        P *= factor(I);
-    return std::log(P.totalSum());
-}
-
-
 vector<VarSet> FactorGraph::Cliques() const {
     vector<VarSet> result;
     
@@ -595,36 +505,4 @@ void FactorGraph::undoProbs( const VarSet &ns ) {
 }
 
 
-bool FactorGraph::isConnected() const {
-    if( nrVars() == 0 )
-        return false;
-    else {
-        Var n = var( 0 );
-
-        VarSet component = n;
-
-        VarSet remaining;
-        for( size_t i = 1; i < nrVars(); i++ )
-            remaining |= var(i);
-
-        bool found_new_vars = true;
-        while( found_new_vars ) {
-            VarSet new_vars;
-            for( VarSet::const_iterator m = remaining.begin(); m != remaining.end(); m++ )
-                if( delta(*m) && component )
-                    new_vars |= *m;
-
-            if( new_vars.empty() )
-                found_new_vars = false;
-            else 
-                found_new_vars = true;
-
-            component |= new_vars;
-            remaining /= new_vars;
-        };
-        return remaining.empty();
-    }
-}
-
-
 } // end of namespace dai
index c141fb1..aaa673b 100644 (file)
@@ -64,24 +64,30 @@ bool HAK::checkProperties() {
 void HAK::constructMessages() {
     // Create outer beliefs
     _Qa.clear();
-    _Qa.reserve(nr_ORs());
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    _Qa.reserve(nrORs());
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         _Qa.push_back( Factor( OR(alpha).vars() ) );
 
     // Create inner beliefs
     _Qb.clear();
-    _Qb.reserve(nr_IRs());
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    _Qb.reserve(nrIRs());
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         _Qb.push_back( Factor( IR(beta) ) );
     
     // Create messages
     _muab.clear();
-    _muab.reserve(nr_Redges());
+    _muab.reserve( nrORs() );
     _muba.clear();
-    _muba.reserve(nr_Redges());
-    for( vector<R_edge_t>::const_iterator ab = Redges().begin(); ab != Redges().end(); ab++ ) {
-        _muab.push_back( Factor( IR(ab->second) ) );
-        _muba.push_back( Factor( IR(ab->second) ) );
+    _muba.reserve( nrORs() );
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
+        _muab.push_back( vector<Factor>() );
+        _muba.push_back( vector<Factor>() );
+        _muab[alpha].reserve( nbOR(alpha).size() );
+        _muba[alpha].reserve( nbOR(alpha).size() );
+        foreach( const Neighbor &beta, nbOR(alpha) ) {
+            _muab[alpha].push_back( Factor( IR(beta) ) );
+            _muba[alpha].push_back( Factor( IR(beta) ) );
+        }
     }
 }
 
@@ -95,7 +101,7 @@ HAK::HAK(const RegionGraph & rg, const Properties &opts) : DAIAlgRG(rg, opts) {
 
 void HAK::findLoopClusters( const FactorGraph & fg, set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
     for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
-        VarSet ind = fg.delta( *in );
+        VarSet ind = fg.delta( fg.findVar( *in ) );
         if( (newcl.size()) >= 2 && (ind >> root) ) {
             allcl.insert( newcl | *in );
         }
@@ -113,14 +119,14 @@ HAK::HAK(const FactorGraph & fg, const Properties &opts) : DAIAlgRG(opts) {
         cl = fg.Cliques();
     } else if( Clusters() == ClustersType::DELTA ) {
         for( size_t i = 0; i < fg.nrVars(); i++ )
-            cl.push_back(fg.Delta(fg.var(i))); 
+            cl.push_back(fg.Delta(i)); 
     } else if( Clusters() == ClustersType::LOOP ) {
         cl = fg.Cliques();
         set<VarSet> scl;
-        for( vector<Var>::const_iterator i0 = fg.vars().begin(); i0 != fg.vars().end(); i0++ ) {
-            VarSet i0d = fg.delta(*i0);
+        for( size_t i0 = 0; i0 < fg.nrVars(); i0++ ) {
+            VarSet i0d = fg.delta(i0);
             if( LoopDepth() > 1 )
-                findLoopClusters( fg, scl, *i0, *i0, LoopDepth() - 1, fg.delta(*i0) );
+                findLoopClusters( fg, scl, fg.var(i0), fg.var(i0), LoopDepth() - 1, fg.delta(i0) );
         }
         for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
             cl.push_back(*c);
@@ -153,12 +159,13 @@ void HAK::init( const VarSet &ns ) {
         if( alpha->vars() && ns )
             alpha->fill( 1.0 / alpha->stateSpace() );
 
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         if( IR(beta) && ns ) {
             _Qb[beta].fill( 1.0 / IR(beta).stateSpace() );
-            for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ ) {
-                muab(*alpha,beta).fill( 1.0 / IR(beta).stateSpace() );
-                muba(beta,*alpha).fill( 1.0 / IR(beta).stateSpace() );
+            foreach( const Neighbor &alpha, nbIR(beta) ) {
+                size_t _beta = alpha.dual;
+                muab( alpha, _beta ).fill( 1.0 / IR(beta).stateSpace() );
+                muba( alpha, _beta ).fill( 1.0 / IR(beta).stateSpace() );
             }
         }
 }
@@ -173,10 +180,12 @@ void HAK::init() {
     for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
         beta->fill( 1.0 / beta->stateSpace() );
 
-    for( size_t ab = 0; ab < nr_Redges(); ab++ ) {
-        _muab[ab].fill( 1.0 / _muab[ab].stateSpace() );
-        _muba[ab].fill( 1.0 / _muba[ab].stateSpace() );
-    }
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
+        foreach( const Neighbor &beta, nbOR(alpha) ) {
+            size_t _beta = beta.iter;
+            muab( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).stateSpace() );
+            muba( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).stateSpace() );
+        }
 }
 
 
@@ -189,7 +198,7 @@ double HAK::doGBP() {
     clock_t tic = toc();
 
     // Check whether counting numbers won't lead to problems
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         assert( nbIR(beta).size() + IR(beta).c() != 0.0 );
 
     // Keep old beliefs to check convergence
@@ -205,13 +214,18 @@ double HAK::doGBP() {
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
     for( iter = 0; iter < MaxIter() && diffs.max() > Tol(); iter++ ) {
-        for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
-            for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ )
-                muab(*alpha,beta) = _Qa[*alpha].marginal(IR(beta)).divided_by( muba(beta,*alpha) );
+        for( size_t beta = 0; beta < nrIRs(); beta++ ) {
+            foreach( const Neighbor &alpha, nbIR(beta) ) {
+                size_t _beta = alpha.dual;
+                muab( alpha, _beta ) = _Qa[alpha].marginal(IR(beta)).divided_by( muba(alpha,_beta) );
+            }
 
             Factor Qb_new;
-            for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ )
-                Qb_new *= muab(*alpha,beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
+            foreach( const Neighbor &alpha, nbIR(beta) ) {
+                size_t _beta = alpha.dual;
+                Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
+            }
+
             Qb_new.normalize( _normtype );
             if( Qb_new.hasNaNs() ) {
                 cout << "HAK::doGBP:  Qb_new has NaNs!" << endl;
@@ -220,20 +234,22 @@ double HAK::doGBP() {
 //          _Qb[beta] = Qb_new.makeZero(1e-100);    // damping?
             _Qb[beta] = Qb_new;
 
-            for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ ) {
-                muba(beta,*alpha) = _Qb[beta].divided_by( muab(*alpha,beta) );
+            foreach( const Neighbor &alpha, nbIR(beta) ) {
+                size_t _beta = alpha.dual;
+
+                muba(alpha,_beta) = _Qb[beta].divided_by( muab(alpha,_beta) );
 
-                Factor Qa_new = OR(*alpha);
-                for( R_nb_cit gamma = nbOR(*alpha).begin(); gamma != nbOR(*alpha).end(); gamma++ )
-                    Qa_new *= muba(*gamma,*alpha);
-                Qa_new ^= (1.0 / OR(*alpha).c());
+                Factor Qa_new = OR(alpha);
+                foreach( const Neighbor &gamma, nbOR(alpha) )
+                    Qa_new *= muba(alpha,gamma.iter);
+                Qa_new ^= (1.0 / OR(alpha).c());
                 Qa_new.normalize( _normtype );
                 if( Qa_new.hasNaNs() ) {
                     cout << "HAK::doGBP:  Qa_new has NaNs!" << endl;
                     return NAN;
                 }
-//              _Qa[*alpha] = Qa_new.makeZero(1e-100); // damping?
-                _Qa[*alpha] = Qa_new;
+//              _Qa[alpha] = Qa_new.makeZero(1e-100); // damping?
+                _Qa[alpha] = Qa_new;
             }
         }
 
@@ -275,11 +291,11 @@ double HAK::doDoubleLoop() {
     clock_t tic = toc();
 
     // Save original outer regions
-    vector<FRegion> org_ORs = ORs();
+    vector<FRegion> org_ORs = ORs;
 
     // Save original inner counting numbers and set negative counting numbers to zero
-    vector<double> org_IR_cs( nr_IRs(), 0.0 );
-    for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
+    vector<double> org_IR_cs( nrIRs(), 0.0 );
+    for( size_t beta = 0; beta < nrIRs(); beta++ ) {
         org_IR_cs[beta] = IR(beta).c();
         if( IR(beta).c() < 0.0 )
             IR(beta).c() = 0.0;
@@ -294,9 +310,9 @@ double HAK::doDoubleLoop() {
     // Differences in single node beliefs
     Diffs diffs(nrVars(), 1.0);
 
-    size_t   outer_maxiter   = MaxIter();
+    size_t  outer_maxiter   = MaxIter();
     double  outer_tol       = Tol();
-    size_t   outer_verbose   = Verbose();
+    size_t  outer_verbose   = Verbose();
     double  org_maxdiff     = MaxDiff();
 
     // Set parameters for inner loop
@@ -306,10 +322,10 @@ double HAK::doDoubleLoop() {
     size_t outer_iter = 0;
     for( outer_iter = 0; outer_iter < outer_maxiter && diffs.max() > outer_tol; outer_iter++ ) {
         // Calculate new outer regions
-        for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+        for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
             OR(alpha) = org_ORs[alpha];
-            for( R_nb_cit beta = nbOR(alpha).begin(); beta != nbOR(alpha).end(); beta++ )
-                OR(alpha) *= _Qb[*beta] ^ ((IR(*beta).c() - org_IR_cs[*beta]) / nbIR(*beta).size());
+            foreach( const Neighbor &beta, nbOR(alpha) )
+                OR(alpha) *= _Qb[beta] ^ ((IR(beta).c() - org_IR_cs[beta]) / nbIR(beta).size());
         }
 
         // Inner loop
@@ -335,10 +351,10 @@ double HAK::doDoubleLoop() {
     updateMaxDiff( diffs.max() );
 
     // Restore original outer regions
-    ORs() = org_ORs;
+    ORs = org_ORs;
 
     // Restore original inner counting numbers
-    for( size_t beta = 0; beta < nr_IRs(); ++beta )
+    for( size_t beta = 0; beta < nrIRs(); ++beta )
         IR(beta).c() = org_IR_cs[beta];
 
     if( Verbose() >= 1 ) {
@@ -390,9 +406,9 @@ Factor HAK::belief( const Var &n ) const {
 
 vector<Factor> HAK::beliefs() const {
     vector<Factor> result;
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         result.push_back( Qb(beta) );
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         result.push_back( Qa(alpha) );
     return result;
 }
@@ -400,9 +416,9 @@ vector<Factor> HAK::beliefs() const {
 
 Complex HAK::logZ() const {
     Complex sum = 0.0;
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         sum += Complex(IR(beta).c()) * Qb(beta).entropy();
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
         sum += Complex(OR(alpha).c()) * Qa(alpha).entropy();
         sum += (OR(alpha).log0() * Qa(alpha)).totalSum();
     }
index b82b0e5..7fb6d4b 100644 (file)
@@ -45,7 +45,7 @@ bool JTree::checkProperties() {
 }
 
 
-JTree::JTree( const FactorGraph &fg, const Properties &opts, bool automatic) : DAIAlgRG(fg, opts), _RTree(), _Qa(), _Qb(), _mes(), _logZ() {
+JTree::JTree( const FactorGraph &fg, const Properties &opts, bool automatic ) : DAIAlgRG(fg, opts), _RTree(), _Qa(), _Qb(), _mes(), _logZ() {
     assert( checkProperties() );
 
     if( automatic ) {
@@ -95,52 +95,58 @@ void JTree::GenerateJT( const vector<VarSet> &Cliques ) {
     // Construct corresponding region graph
 
     // Create outer regions
-    ORs().reserve( Cliques.size() );
+    ORs.reserve( Cliques.size() );
     for( size_t i = 0; i < Cliques.size(); i++ )
-        ORs().push_back( FRegion( Factor(Cliques[i], 1.0), 1.0 ) );
+        ORs.push_back( FRegion( Factor(Cliques[i], 1.0), 1.0 ) );
 
     // For each factor, find an outer region that subsumes that factor.
     // Then, multiply the outer region with that factor.
     for( size_t I = 0; I < nrFactors(); I++ ) {
         size_t alpha;
-        for( alpha = 0; alpha < nr_ORs(); alpha++ )
+        for( alpha = 0; alpha < nrORs(); alpha++ )
             if( OR(alpha).vars() >> factor(I).vars() ) {
 //              OR(alpha) *= factor(I);
-                _fac2OR[I] = alpha;
+                fac2OR.push_back( alpha );
                 break;
             }
-        assert( alpha != nr_ORs() );
+        assert( alpha != nrORs() );
     }
     RecomputeORs();
 
     // Create inner regions and edges
-    IRs().reserve( _RTree.size() );
-    Redges().reserve( 2 * _RTree.size() );
+    IRs.reserve( _RTree.size() );
+    typedef pair<size_t,size_t> Edge;
+    vector<Edge> edges;
+    edges.reserve( 2 * _RTree.size() );
     for( size_t i = 0; i < _RTree.size(); i++ ) {
-        Redges().push_back( R_edge_t( _RTree[i].n1, IRs().size() ) );
-        Redges().push_back( R_edge_t( _RTree[i].n2, IRs().size() ) );
+        edges.push_back( Edge( _RTree[i].n1, nrIRs() ) );
+        edges.push_back( Edge( _RTree[i].n2, nrIRs() ) );
         // inner clusters have counting number -1
-        IRs().push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
+        IRs.push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
     }
 
-    // Regenerate BipartiteGraph internals
-    Regenerate();
+    // create bipartite graph
+    G.create( nrORs(), nrIRs(), edges.begin(), edges.end() );
 
     // Create messages and beliefs
     _Qa.clear();
-    _Qa.reserve( nr_ORs() );
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    _Qa.reserve( nrORs() );
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         _Qa.push_back( OR(alpha) );
 
     _Qb.clear();
-    _Qb.reserve( nr_IRs() );
-    for( size_t beta = 0; beta < nr_IRs(); beta++ ) 
+    _Qb.reserve( nrIRs() );
+    for( size_t beta = 0; beta < nrIRs(); beta++ ) 
         _Qb.push_back( Factor( IR(beta), 1.0 ) );
 
     _mes.clear();
-    _mes.reserve( nr_Redges() );
-    for( size_t e = 0; e < nr_Redges(); e++ )
-        _mes.push_back( Factor( IR(Redge(e).second), 1.0 ) );
+    _mes.reserve( nrORs() );
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
+        _mes.push_back( vector<Factor>() );
+        _mes[alpha].reserve( nbOR(alpha).size() );
+        foreach( const Neighbor &beta, nbOR(alpha) )
+            _mes[alpha].push_back( Factor( IR(beta), 1.0 ) );
+    }
 
     // Check counting numbers
     Check_Counting_Numbers();
@@ -178,9 +184,9 @@ Factor JTree::belief( const VarSet &ns ) const {
 
 vector<Factor> JTree::beliefs() const {
     vector<Factor> result;
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         result.push_back( _Qb[beta] );
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         result.push_back( _Qa[alpha] );
     return result;
 }
@@ -193,10 +199,10 @@ Factor JTree::belief( const Var &n ) const {
 
 // Needs no init
 void JTree::runHUGIN() {
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         _Qa[alpha] = OR(alpha);
 
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         _Qb[beta].fill( 1.0 );
 
     // CollectEvidence
@@ -224,7 +230,7 @@ void JTree::runHUGIN() {
     }
 
     // Normalize
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         _Qa[alpha].normalize( Prob::NORMPROB );
 }
 
@@ -233,42 +239,44 @@ void JTree::runHUGIN() {
 void JTree::runShaferShenoy() {
     // First pass
     _logZ = 0.0;
-    for( size_t e = _RTree.size(); (e--) != 0; ) {
+    for( size_t e = nrIRs(); (e--) != 0; ) {
         // send a message from _RTree[e].n2 to _RTree[e].n1
         // or, actually, from the seperator IR(e) to _RTree[e].n1
 
-        size_t i = _RTree[e].n2;
-        size_t j = _RTree[e].n1;
+        size_t i = nbIR(e)[1].node; // = _RTree[e].n2
+        size_t j = nbIR(e)[0].node; // = _RTree[e].n1
+        size_t _e = nbIR(e)[0].dual;
         
         Factor piet = OR(i);
-        for( R_nb_cit k = nbOR(i).begin(); k != nbOR(i).end(); k++ )
-            if( *k != e )
-                piet *= message( i, *k );
-        message( j, e ) = piet.part_sum( IR(e) );
-        _logZ += log( message(j,e).normalize( Prob::NORMPROB ) );
+        foreach( const Neighbor &k, nbOR(i) )
+            if( k != e ) 
+                piet *= message( i, k.iter );
+        message( j, _e ) = piet.part_sum( IR(e) );
+        _logZ += log( message(j,_e).normalize( Prob::NORMPROB ) );
     }
 
     // Second pass
-    for( size_t e = 0; e < _RTree.size(); e++ ) {
-        size_t i = _RTree[e].n1;
-        size_t j = _RTree[e].n2;
+    for( size_t e = 0; e < nrIRs(); e++ ) {
+        size_t i = nbIR(e)[0].node; // = _RTree[e].n1
+        size_t j = nbIR(e)[1].node; // = _RTree[e].n2
+        size_t _e = nbIR(e)[1].dual;
         
         Factor piet = OR(i);
-        for( R_nb_cit k = nbOR(i).begin(); k != nbOR(i).end(); k++ )
-            if( *k != e )
-                piet *= message( i, *k );
-        message( j, e ) = piet.marginal( IR(e) );
+        foreach( const Neighbor &k, nbOR(i) )
+            if(  k != e )
+                piet *= message( i, k.iter );
+        message( j, _e ) = piet.marginal( IR(e) );
     }
 
     // Calculate beliefs
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
         Factor piet = OR(alpha);
-        for( R_nb_cit k = nbOR(alpha).begin(); k != nbOR(alpha).end(); k++ )
-            piet *= message( alpha, *k );
-        if( _RTree.empty() ) {
+        foreach( const Neighbor &k, nbOR(alpha) )
+            piet *= message( alpha, k.iter );
+        if( nrIRs() == 0 ) {
             _logZ += log( piet.normalize( Prob::NORMPROB ) );
             _Qa[alpha] = piet;
-        } else if( alpha == _RTree[0].n1 ) {
+        } else if( alpha == nbIR(0)[0].node /*_RTree[0].n1*/ ) {
             _logZ += log( piet.normalize( Prob::NORMPROB ) );
             _Qa[alpha] = piet;
         } else
@@ -276,8 +284,8 @@ void JTree::runShaferShenoy() {
     }
 
     // Only for logZ (and for belief)...
-    for( size_t beta = 0; beta < nr_IRs(); beta++ ) 
-        _Qb[beta] = _Qa[nbIR(beta)[0]].marginal( IR(beta) );
+    for( size_t beta = 0; beta < nrIRs(); beta++ ) 
+        _Qb[beta] = _Qa[nbIR(beta)[0].node].marginal( IR(beta) );
 }
 
 
@@ -292,9 +300,9 @@ double JTree::run() {
 
 Complex JTree::logZ() const {
     Complex sum = 0.0;
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         sum += Complex(IR(beta).c()) * _Qb[beta].entropy();
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
         sum += Complex(OR(alpha).c()) * _Qa[alpha].entropy();
         sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
     }
@@ -306,7 +314,7 @@ Complex JTree::logZ() const {
 size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t PreviousRoot ) const {
     // find new root clique (the one with maximal statespace overlap with ns)
     size_t maxval = 0, maxalpha = 0;
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
         size_t val = (ns & OR(alpha).vars()).stateSpace();
         if( val > maxval ) {
             maxval = val;
@@ -453,10 +461,10 @@ Factor JTree::calcMarginal( const VarSet& ns ) {
                 size_t alpha1 = T[i].n1;
                 size_t alpha2 = T[i].n2;
                 size_t beta;
-                for( beta = 0; beta < nr_IRs(); beta++ )
+                for( beta = 0; beta < nrIRs(); beta++ )
                     if( UEdge( _RTree[beta].n1, _RTree[beta].n2 ) == UEdge( alpha1, alpha2 ) )
                         break;
-                assert( beta != nr_IRs() );
+                assert( beta != nrIRs() );
                 b[i] = beta;
 
                 if( !_Qa_old.count( alpha1 ) )
index fb125da..2748550 100644 (file)
@@ -71,30 +71,20 @@ bool LC::checkProperties() {
 LC::LC(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
     assert( checkProperties() );
 
-    // calc iI
-    for( size_t i=0; i < nrVars(); i++ ) {
-        for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ ) {
-                _iI_type _iI_entry;
-                _iI_entry.i = i;
-                _iI_entry.I = *I;
-
-                _iI.push_back(_iI_entry);
-            }
-    }
-
     // create pancakes
     _pancakes.resize(nrVars());
    
     // create cavitydists
     for( size_t i=0; i < nrVars(); i++ )
-        _cavitydists.push_back(Factor(delta(var(i))));
+        _cavitydists.push_back(Factor(delta(i)));
 
     // create phis
-    _phis.reserve(nr_edges());
-    for( size_t iI = 0; iI < nr_edges(); iI++ ) {
-        size_t i = edge(iI).first;
-        size_t I = edge(iI).second;
-        _phis.push_back( Factor( factor(I).vars() / var(i) ) );
+    _phis.reserve( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ ) {
+        _phis.push_back( vector<Factor>() );
+        _phis[i].reserve( nbV(i).size() );
+        foreach( const Neighbor &I, nbV(i) )
+            _phis[i].push_back( Factor( factor(I).vars() / var(i) ) );
     }
 
     // create beliefs
@@ -120,27 +110,27 @@ double LC::CalcCavityDist (size_t i, const string &name, const Properties &opts)
     double maxdiff = 0;
 
     if( Verbose() >= 2 )
-        cout << "Initing cavity " << var(i) << "(" << delta(var(i)).size() << " vars, " << delta(var(i)).stateSpace() << " states)" << endl;
+        cout << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).stateSpace() << " states)" << endl;
 
     if( Cavity() == CavityType::UNIFORM )
-        Bi = Factor(delta(var(i)));
+        Bi = Factor(delta(i));
     else {
         InfAlg *cav = newInfAlg( name, *this, opts );
-        cav->makeCavity( var(i) );
+        cav->makeCavity( i );
 
         if( Cavity() == CavityType::FULL )
-            Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
+            Bi = calcMarginal( *cav, cav->delta(i), reInit() );
         else if( Cavity() == CavityType::PAIR )
-            Bi = calcMarginal2ndO( *cav, cav->delta(var(i)), reInit() );
+            Bi = calcMarginal2ndO( *cav, cav->delta(i), reInit() );
         else if( Cavity() == CavityType::PAIR2 ) {
-            vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->delta(var(i)), reInit() );
+            vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->delta(i), reInit() );
             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
                 Bi *= pairbeliefs[ij];
         } else if( Cavity() == CavityType::PAIRINT ) {
-            Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
+            Bi = calcMarginal( *cav, cav->delta(i), reInit() );
             
             // Set interactions of order > 2 to zero
-            size_t N = delta(var(i)).size();
+            size_t N = delta(i).size();
             Real *p = &(*Bi.p().p().begin());
             x2x::p2logp (N, p);
             x2x::logp2w (N, p);
@@ -149,10 +139,10 @@ double LC::CalcCavityDist (size_t i, const string &name, const Properties &opts)
 //            x2x::logpnorm (N, p);
             x2x::logp2p (N, p);
         } else if( Cavity() == CavityType::PAIRCUM ) {
-            Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
+            Bi = calcMarginal( *cav, cav->delta(i), reInit() );
             
             // Set cumulants of order > 2 to zero
-            size_t N = delta(var(i)).size();
+            size_t N = delta(i).size();
             Real *p = &(*Bi.p().p().begin());
             x2x::p2m (N, p);
             x2x::m2c (N, p, N);
@@ -218,19 +208,19 @@ long LC::SetCavityDists (vector<Factor> &Q) {
 
 
 void LC::init() {
-    for( size_t iI = 0; iI < nr_edges(); iI++ ) {
-        if( Updates() == UpdateType::SEQRND )
-            _phis[iI].randomize();
-        else
-            _phis[iI].fill(1.0);
-    }
+    for( size_t i = 0; i < nrVars(); ++i )
+        foreach( const Neighbor &I, nbV(i) )
+            if( Updates() == UpdateType::SEQRND )
+                _phis[i][I.iter].randomize();
+            else
+                _phis[i][I.iter].fill(1.0);
     for( size_t i = 0; i < nrVars(); i++ ) {
         _pancakes[i] = _cavitydists[i];
         
-        for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ ) {
-            _pancakes[i] *= factor(*I);
+        foreach( const Neighbor &I, nbV(i) ) {
+            _pancakes[i] *= factor(I);
             if( Updates() == UpdateType::SEQRND )
-              _pancakes[i] *= _phis[VV2E(i,*I)];
+              _pancakes[i] *= _phis[i][I.iter];
         }
         
         _pancakes[i].normalize( _normtype );
@@ -240,11 +230,8 @@ void LC::init() {
 }
 
 
-Factor LC::NewPancake (size_t iI, bool & hasNaNs) {
-    size_t i = _iI[iI].i;
-    size_t I = _iI[iI].I;
-    iI = VV2E(i, I);
-
+Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
+    size_t I = nbV(i)[_I];
     Factor piet = _pancakes[i];
 
     // recalculate _pancake[i]
@@ -255,16 +242,16 @@ Factor LC::NewPancake (size_t iI, bool & hasNaNs) {
             A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).part_sum( Ivars / var(i) );
     if( Ivars.size() > 1 )
         A_I ^= (1.0 / (Ivars.size() - 1));
-    Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[iI].inverse()).part_sum( Ivars / var(i) );
+    Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).part_sum( Ivars / var(i) );
     Factor quot = A_I.divided_by(A_Ii);
 
-    piet *= quot.divided_by( _phis[iI] ).normalized( _normtype );
-    _phis[iI] = quot.normalized( _normtype );
+    piet *= quot.divided_by( _phis[i][_I] ).normalized( _normtype );
+    _phis[i][_I] = quot.normalized( _normtype );
 
     piet.normalize( _normtype );
 
     if( piet.hasNaNs() ) {
-        cout << "LC::NewPancake(" << iI << "):  has NaNs!" << endl;
+        cout << "LC::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
         hasNaNs = true;
     }
 
@@ -298,11 +285,15 @@ double LC::run() {
         return NAN;
     }
 
-    vector<long> update_seq(nr_iI(),0);
-    for( size_t k=0; k < nr_iI(); k++ )
-        update_seq[k] = k;
+    size_t nredges = nrEdges();
+    typedef pair<size_t,size_t> Edge;
+    vector<Edge> update_seq;
+    update_seq.reserve( nredges );
+    for( size_t i = 0; i < nrVars(); ++i )
+        foreach( const Neighbor &I, nbV(i) )
+            update_seq.push_back( Edge( i, I.iter ) );
 
-    size_t iter=0;
+    size_t iter = 0;
 
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
@@ -311,13 +302,13 @@ double LC::run() {
         if( Updates() == UpdateType::SEQRND )
             random_shuffle( update_seq.begin(), update_seq.end() );
         
-        for( size_t t=0; t < nr_iI(); t++ ) {
-            long iI = update_seq[t];
-            long i = _iI[iI].i;
-            _pancakes[i] = NewPancake(iI, hasNaNs);
+        for( size_t t=0; t < nredges; t++ ) {
+            size_t i = update_seq[t].first;
+            size_t _I = update_seq[t].second;
+            _pancakes[i] = NewPancake( i, _I, hasNaNs);
             if( hasNaNs )
                 return NAN;
-            CalcBelief(i);
+            CalcBelief( i );
         }
 
         // compare new beliefs with old ones
index e95006e..7fe6958 100644 (file)
@@ -54,15 +54,15 @@ bool MF::checkProperties() {
 
 
 void MF::Regenerate() {
-    DAIAlgFG::Regenerate();
+//    DAIAlgFG::Regenerate();
 
     // clear beliefs
     _beliefs.clear();
     _beliefs.reserve( nrVars() );
 
     // create beliefs
-    for( vector<Var>::const_iterator i = vars().begin(); i != vars().end(); i++ ) 
-        _beliefs.push_back(Factor(*i));
+    for( size_t i = 0; i < nrVars(); ++i )
+        _beliefs.push_back(Factor(var(i)));
 }
 
 
@@ -97,19 +97,18 @@ double MF::run() {
 
         Factor jan;
         Factor piet;
-        for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ ) {
-
+        foreach( const Neighbor &I, nbV(i) ) {
             Factor henk;
-            for( _nb_cit j = nb2(*I).begin(); j != nb2(*I).end(); j++ ) // for all j in I \ i
-                if( *j != i )
-                    henk *= _beliefs[*j];
-            piet = factor(*I).log0();
+            foreach( const Neighbor &j, nbF(I) ) // for all j in I \ i
+                if( j != i )
+                    henk *= _beliefs[j];
+            piet = factor(I).log0();
             piet *= henk;
             piet = piet.part_sum(var(i));
             piet = piet.exp();
             jan *= piet; 
         }
-            
+
         jan.normalize( _normtype );
 
         if( jan.hasNaNs() ) {
@@ -140,7 +139,7 @@ double MF::run() {
 }
 
 
-Factor MF::belief1 (size_t i) const {
+Factor MF::beliefV (size_t i) const {
     Factor piet;
     piet = _beliefs[i];
     piet.normalize( Prob::NORMPROB );
@@ -159,14 +158,14 @@ Factor MF::belief (const VarSet &ns) const {
 
 
 Factor MF::belief (const Var &n) const {
-    return( belief1( findVar( n) ) );
+    return( beliefV( findVar( n ) ) );
 }
 
 
 vector<Factor> MF::beliefs() const {
     vector<Factor> result;
     for( size_t i = 0; i < nrVars(); i++ )
-        result.push_back( belief1(i) );
+        result.push_back( beliefV(i) );
     return result;
 }
 
@@ -175,11 +174,11 @@ Complex MF::logZ() const {
     Complex sum = 0.0;
     
     for(size_t i=0; i < nrVars(); i++ )
-        sum -= belief1(i).entropy();
+        sum -= beliefV(i).entropy();
     for(size_t I=0; I < nrFactors(); I++ ) {
         Factor henk;
-        for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); j++ )   // for all j in I
-            henk *= _beliefs[*j];
+        foreach( const Neighbor &j, nbF(I) )  // for all j in I
+            henk *= _beliefs[j];
         henk.normalize( Prob::NORMPROB );
         Factor piet;
         piet = factor(I).log0();
index 5c6efee..cf80165 100644 (file)
@@ -492,12 +492,12 @@ void MR::init_cor() {
         vector<Factor> pairq;
         if( Inits() == InitType::CLAMPING ) {
             BP bpcav(*this, Properties()("updates",string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("1000UL"))("verbose", string("0UL")));
-            bpcav.makeCavity( var(i) );
-            pairq = calcPairBeliefs( bpcav, delta(var(i)), false );
+            bpcav.makeCavity( i );
+            pairq = calcPairBeliefs( bpcav, delta(i), false );
         } else if( Inits() == InitType::EXACT ) {
             JTree jtcav(*this, Properties()("updates",string("HUGIN"))("verbose", string("0UL")) );
-            jtcav.makeCavity( var(i) );
-            pairq = calcPairBeliefs( jtcav, delta(var(i)), false );
+            jtcav.makeCavity( i );
+            pairq = calcPairBeliefs( jtcav, delta(i), false );
         }
         for( size_t jk = 0; jk < pairq.size(); jk++ ) {
             VarSet::const_iterator kit = pairq[jk].vars().begin();
@@ -605,7 +605,7 @@ MR::MR( const FactorGraph &fg, const Properties &opts ) : DAIAlgFG(fg, opts), su
     // check whether all vars in fg are binary
     // check whether connectivity is <= kmax
     for( size_t i = 0; i < fg.nrVars(); i++ )
-        if( (fg.var(i).states() > 2) || (fg.delta(fg.var(i)).size() > kmax) ) {
+        if( (fg.var(i).states() > 2) || (fg.delta(i).size() > kmax) ) {
             supported = false;
             break;
         }
index 466cc99..14a2714 100644 (file)
@@ -40,22 +40,34 @@ std::ostream& operator<< (std::ostream & os, const Property & p) {
         os << boost::any_cast<bool>(p.second);
     else if( p.second.type() == typeid(Properties) )
         os << boost::any_cast<Properties>(p.second);
+#ifdef WITH_BP
     else if( p.second.type() == typeid(BP::UpdateType) )
         os << boost::any_cast<BP::UpdateType>(p.second);
+#endif
+#ifdef WITH_HAK
     else if( p.second.type() == typeid(HAK::ClustersType) )
         os << boost::any_cast<HAK::ClustersType>(p.second);
+#endif
+#ifdef WITH_JTREE
     else if( p.second.type() == typeid(JTree::UpdateType) )
         os << boost::any_cast<JTree::UpdateType>(p.second);
+#endif
+#ifdef WITH_MR
     else if( p.second.type() == typeid(MR::UpdateType) )
         os << boost::any_cast<MR::UpdateType>(p.second);
     else if( p.second.type() == typeid(MR::InitType) )
         os << boost::any_cast<MR::InitType>(p.second);
+#endif
+#ifdef WITH_TREEEP
     else if( p.second.type() == typeid(TreeEP::TypeType) )
         os << boost::any_cast<TreeEP::TypeType>(p.second);
+#endif
+#ifdef WITH_LC
     else if( p.second.type() == typeid(LC::CavityType) )
         os << boost::any_cast<LC::CavityType>(p.second);
     else if( p.second.type() == typeid(LC::UpdateType) )
         os << boost::any_cast<LC::UpdateType>(p.second);
+#endif
     else
         throw "Unknown property type";
     return( os );
index 9a3e624..e841554 100644 (file)
@@ -32,110 +32,62 @@ namespace dai {
 using namespace std;
 
 
-RegionGraph::RegionGraph(const FactorGraph & fg, const vector<Region> & ors, const vector<Region> & irs, const vector<R_edge_t> & edges) : FactorGraph(fg), BipRegGraph(), _fac2OR() {
+RegionGraph::RegionGraph( const FactorGraph &fg, const vector<Region> &ors, const vector<Region> &irs, const vector<pair<size_t,size_t> > &edges) : FactorGraph(fg), G(), ORs(), IRs(irs), fac2OR() {
     // Copy outer regions (giving them counting number 1.0)
-    ORs().reserve( ors.size() );
+    ORs.reserve( ors.size() );
     for( vector<Region>::const_iterator alpha = ors.begin(); alpha != ors.end(); alpha++ )
-        ORs().push_back( FRegion(Factor(*alpha, 1.0), 1.0) );
-
-    // Incorporate factors into outer regions
-/*  if( !_hasNegatives ) {
-        // For each factor, find the outer regions that subsume that factor.
-        // Then, "divide" the factor over its subsuming outer regions.
-        for( size_t I = 0; I < nrFactors(); I++ ) {
-            vector<size_t> subsuming_alphas;
-
-            for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
-                if( OR(alpha).vars() >> factor(I).vars() )
-                    subsuming_alphas.push_back( alpha );
-            }
-
-            assert( subsuming_alphas.size() >= 1 );
+        ORs.push_back( FRegion(Factor(*alpha, 1.0), 1.0) );
 
-            for( vector<size_t>::const_iterator alpha = subsuming_alphas.begin(); alpha != subsuming_alphas.end(); alpha++ )
-                OR(*alpha) *= (factor(I) ^ (1.0 / subsuming_alphas.size()));
-        }
-    } else {
-        cout << "Careful composition of outer regions" << endl;
-*/
     // For each factor, find an outer regions that subsumes that factor.
     // Then, multiply the outer region with that factor.
+    fac2OR.reserve( nrFactors() );
     for( size_t I = 0; I < nrFactors(); I++ ) {
         size_t alpha;
-        for( alpha = 0; alpha < nr_ORs(); alpha++ )
+        for( alpha = 0; alpha < nrORs(); alpha++ )
             if( OR(alpha).vars() >> factor(I).vars() ) {
-                _fac2OR[I] = alpha;
+                fac2OR.push_back( alpha );
 //              OR(alpha) *= factor(I);
                 break;
             }
-        assert( alpha != nr_ORs() );
+        assert( alpha != nrORs() );
     }
     RecomputeORs();
-//  }
-    
-    // Copy inner regions
-    IRs().reserve( irs.size() );
-    for( vector<Region>::const_iterator beta = irs.begin(); beta != irs.end(); beta++ )
-        IRs().push_back( *beta );
     
-    // Copy edges
-    Redges().reserve( edges.size() );
-    for( vector<R_edge_t>::const_iterator e = edges.begin(); e != edges.end(); e++ )
-        Redges().push_back( *e );
-    
-    // Regenerate BipartiteGraph internals
-    BipRegGraph::Regenerate();
+    // create bipartite graph
+    G.create( nrORs(), nrIRs(), edges.begin(), edges.end() );
 
     // Check counting numbers
+#ifdef DAI_DEBUG
     Check_Counting_Numbers();
+#endif
 }
 
 
 // CVM style
-RegionGraph::RegionGraph(const FactorGraph & fg, const vector<VarSet> & cl) : FactorGraph(fg), BipRegGraph() {
+RegionGraph::RegionGraph( const FactorGraph &fg, const vector<VarSet> &cl ) : FactorGraph(fg), G(), ORs(), IRs(), fac2OR() {
     // Retain only maximal clusters
     ClusterGraph cg( cl );
     cg.eraseNonMaximal();
     
     // Create outer regions, giving them counting number 1.0
-    ORs().reserve( cg.size() );
+    ORs.reserve( cg.size() );
     for( ClusterGraph::const_iterator alpha = cg.begin(); alpha != cg.end(); alpha++ )
-        ORs().push_back( FRegion(Factor(*alpha, 1.0), 1.0) );
-
-    // Incorporate factors into outer regions
-/*  if( !_hasNegatives ) {
-        // For each factor, find the outer regions that subsume that factor.
-        // Then, "divide" the factor over its subsuming outer regions.
-        for( size_t I = 0; I < nrFactors(); I++ ) {
-            vector<size_t> subsuming_alphas;
-
-            for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
-                if( OR(alpha).vars() >> factor(I).vars() )
-                    subsuming_alphas.push_back( alpha );
-            }
-
-            assert( subsuming_alphas.size() >= 1 );
+        ORs.push_back( FRegion(Factor(*alpha, 1.0), 1.0) );
 
-            for( vector<size_t>::const_iterator alpha = subsuming_alphas.begin(); alpha != subsuming_alphas.end(); alpha++ )
-                OR(*alpha) *= factor(I) ^ (1.0 / subsuming_alphas.size());
-        }
-    } else {
-        cout << "Careful composition of outer regions" << endl;
-*/
     // For each factor, find an outer regions that subsumes that factor.
     // Then, multiply the outer region with that factor.
+    fac2OR.reserve( nrFactors() );
     for( size_t I = 0; I < nrFactors(); I++ ) {
         size_t alpha;
-        for( alpha = 0; alpha < nr_ORs(); alpha++ )
+        for( alpha = 0; alpha < nrORs(); alpha++ )
             if( OR(alpha).vars() >> factor(I).vars() ) {
-                _fac2OR[I] = alpha;
+                fac2OR.push_back( alpha );
 //              OR(alpha) *= factor(I);
                 break;
             }
-        assert( alpha != nr_ORs() );
+        assert( alpha != nrORs() );
     }
     RecomputeORs();
-//  }
     
     // Create inner regions - first pass
     set<VarSet> betas;
@@ -160,36 +112,39 @@ RegionGraph::RegionGraph(const FactorGraph & fg, const vector<VarSet> & cl) : Fa
     } while( new_betas.size() );
 
     // Create inner regions - store them in the bipartite graph
-    IRs().reserve( betas.size() );
+    IRs.reserve( betas.size() );
     for( set<VarSet>::const_iterator beta = betas.begin(); beta != betas.end(); beta++ )
-        IRs().push_back( Region(*beta,NAN) );
+        IRs.push_back( Region(*beta,NAN) );
     
     // Create edges
-    for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
-        for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+    vector<pair<size_t,size_t> > edges;
+    for( size_t beta = 0; beta < nrIRs(); beta++ ) {
+        for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
             if( OR(alpha).vars() >> IR(beta) )
-                Redges().push_back(R_edge_t(alpha,beta));
+                edges.push_back( pair<size_t,size_t>(alpha,beta) );
         }
     }
     
-    // Regenerate BipartiteGraph internals
-    BipRegGraph::Regenerate();
+    // create bipartite graph
+    G.create( nrORs(), nrIRs(), edges.begin(), edges.end() );
 
     // Calculate counting numbers
     Calc_Counting_Numbers();
 
     // Check counting numbers
+#ifdef DAI_DEBUG
     Check_Counting_Numbers();
+#endif
 }
 
 
 void RegionGraph::Calc_Counting_Numbers() {
     // Calculates counting numbers of inner regions based upon counting numbers of outer regions
     
-    vector<vector<size_t> > ancestors(nr_IRs());
-    for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
+    vector<vector<size_t> > ancestors(nrIRs());
+    for( size_t beta = 0; beta < nrIRs(); beta++ ) {
         IR(beta).c() = NAN;
-        for( size_t beta2 = 0; beta2 < nr_IRs(); beta2++ )
+        for( size_t beta2 = 0; beta2 < nrIRs(); beta2++ )
             if( (beta2 != beta) && IR(beta2) >> IR(beta) )
                 ancestors[beta].push_back(beta2);
     }
@@ -197,7 +152,7 @@ void RegionGraph::Calc_Counting_Numbers() {
     bool new_counting;
     do {
         new_counting = false;
-        for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
+        for( size_t beta = 0; beta < nrIRs(); beta++ ) {
             if( isnan( IR(beta).c() ) ) {
                 bool has_nan_ancestor = false;
                 for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); (beta2 != ancestors[beta].end()) && !has_nan_ancestor; beta2++ )
@@ -205,8 +160,8 @@ void RegionGraph::Calc_Counting_Numbers() {
                         has_nan_ancestor = true;
                 if( !has_nan_ancestor ) {
                     double c = 1.0;
-                    for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ )
-                        c -= OR(*alpha).c();
+                    foreach( const Neighbor &alpha, nbIR(beta) )
+                        c -= OR(alpha).c();
                     for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); beta2 != ancestors[beta].end(); beta2++ )
                         c -= IR(*beta2).c();
                     IR(beta).c() = c;
@@ -222,12 +177,12 @@ bool RegionGraph::Check_Counting_Numbers() {
     // Checks whether the counting numbers satisfy the fundamental relation
     
     bool all_valid = true;
-    for( vector<Var>::const_iterator n = vars().begin(); n != vars().end(); n++ ) {
+    for( vector<Var>::const_iterator n = vars.begin(); n != vars.end(); n++ ) {
         double c_n = 0.0;
-        for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+        for( size_t alpha = 0; alpha < nrORs(); alpha++ )
             if( OR(alpha).vars() && *n )
                 c_n += OR(alpha).c();
-        for( size_t beta = 0; beta < nr_IRs(); beta++ )
+        for( size_t beta = 0; beta < nrIRs(); beta++ )
             if( IR(beta) && *n )
                 c_n += IR(beta).c();
         if( fabs(c_n - 1.0) > 1e-15 ) {
@@ -241,41 +196,44 @@ bool RegionGraph::Check_Counting_Numbers() {
 
 
 void RegionGraph::RecomputeORs() {
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         OR(alpha).fill( 1.0 );
-    for( fac2OR_cit I = _fac2OR.begin(); I != _fac2OR.end(); I++ )
-        OR( I->second ) *= factor( I->first );
+    for( size_t I = 0; I < nrFactors(); I++ )
+        if( fac2OR[I] != -1U )
+            OR( fac2OR[I] ) *= factor( I );
 }
 
 
 void RegionGraph::RecomputeORs( const VarSet &ns ) {
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         if( OR(alpha).vars() && ns )
             OR(alpha).fill( 1.0 );
-    for( fac2OR_cit I = _fac2OR.begin(); I != _fac2OR.end(); I++ )
-        if( OR( I->second ).vars() && ns )
-            OR( I->second ) *= factor( I->first );
+    for( size_t I = 0; I < nrFactors(); I++ )
+        if( fac2OR[I] != -1U )
+            if( OR( fac2OR[I] ).vars() && ns )
+                OR( fac2OR[I] ) *= factor( I );
 }
 
 
 void RegionGraph::RecomputeOR( size_t I ) {
-    if( _fac2OR.count(I) ) {
-        size_t alpha = _fac2OR[I];
+    assert( I < nrFactors() );
+    if( fac2OR[I] != -1U ) {
+        size_t alpha = fac2OR[I];
         OR(alpha).fill( 1.0 );
-        for( fac2OR_cit I = _fac2OR.begin(); I != _fac2OR.end(); I++ )
-            if( I->second == alpha )
-                OR(alpha) *= factor( I->first );
+        for( size_t J = 0; J < nrFactors(); J++ )
+            if( fac2OR[J] == alpha )
+                OR(alpha) *= factor( J );
     }
 }
 
 
 ostream & operator << (ostream & os, const RegionGraph & rg) {
     os << "Outer regions" << endl;
-    for( size_t alpha = 0; alpha < rg.nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < rg.nrORs(); alpha++ )
         os << rg.OR(alpha).vars() << ": c = " << rg.OR(alpha).c() << endl;
 
     os << "Inner regions" << endl;
-    for( size_t beta = 0; beta < rg.nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < rg.nrIRs(); beta++ )
         os << (VarSet)rg.IR(beta) << ": c = " << rg.IR(beta).c() << endl;
 
     return(os);
index f0b0ba7..a12ff1b 100644 (file)
@@ -190,7 +190,7 @@ double TreeEPSubTree::logZ( const vector<Factor> &Qa, const vector<Factor> &Qb )
 TreeEP::TreeEP( const FactorGraph &fg, const Properties &opts ) : JTree(fg, opts("updates",string("HUGIN")), false) {
     assert( checkProperties() );
 
-    assert( fg.isConnected() );
+    assert( fg.G.isConnected() );
 
     if( opts.hasKey("tree") ) {
         ConstructRG( opts.GetAs<DEdgeVec>("tree") );
@@ -199,24 +199,25 @@ TreeEP::TreeEP( const FactorGraph &fg, const Properties &opts ) : JTree(fg, opts
             // construct weighted graph with as weights a crude estimate of the
             // mutual information between the nodes
             WeightedGraph<double> wg;
-            for( vector<Var>::const_iterator i = vars().begin(); i != vars().end(); i++ ) {
-                VarSet di = delta(*i);
+            for( size_t i = 0; i < nrVars(); ++i ) {
+                Var v_i = var(i);
+                VarSet di = delta(i);
                 for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
-                    if( *i < *j ) {
+                    if( v_i < *j ) {
                         Factor piet;
                         for( size_t I = 0; I < nrFactors(); I++ ) {
                             VarSet Ivars = factor(I).vars();
-                            if( (Ivars == *i) || (Ivars == *j) )
+                            if( (Ivars == v_i) || (Ivars == *j) )
                                 piet *= factor(I);
-                            else if( Ivars >> (*i | *j) )
-                                piet *= factor(I).marginal( *i | *j );
+                            else if( Ivars >> (v_i | *j) )
+                                piet *= factor(I).marginal( v_i | *j );
                         }
-                        if( piet.vars() >> (*i | *j) ) {
-                            piet = piet.marginal( *i | *j );
-                            Factor pietf = piet.marginal(*i) * piet.marginal(*j);
-                            wg[UEdge(findVar(*i),findVar(*j))] = real( KL_dist( piet, pietf ) );
+                        if( piet.vars() >> (v_i | *j) ) {
+                            piet = piet.marginal( v_i | *j );
+                            Factor pietf = piet.marginal(v_i) * piet.marginal(*j);
+                            wg[UEdge(i,findVar(*j))] = real( KL_dist( piet, pietf ) );
                         } else
-                            wg[UEdge(findVar(*i),findVar(*j))] = 0;
+                            wg[UEdge(i,findVar(*j))] = 0;
                     }
             }
 
@@ -233,17 +234,18 @@ TreeEP::TreeEP( const FactorGraph &fg, const Properties &opts ) : JTree(fg, opts
             // construct weighted graph with as weights an upper bound on the
             // effective interaction strength between pairs of nodes
             WeightedGraph<double> wg;
-            for( vector<Var>::const_iterator i = vars().begin(); i != vars().end(); i++ ) {
-                VarSet di = delta(*i);
+            for( size_t i = 0; i < nrVars(); ++i ) {
+                Var v_i = var(i);
+                VarSet di = delta(i);
                 for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
-                    if( *i < *j ) {
+                    if( v_i < *j ) {
                         Factor piet;
                         for( size_t I = 0; I < nrFactors(); I++ ) {
                             VarSet Ivars = factor(I).vars();
-                            if( Ivars >> (*i | *j) )
+                            if( Ivars >> (v_i | *j) )
                                 piet *= factor(I);
                         }
-                        wg[UEdge(findVar(*i),findVar(*j))] = piet.strength(*i, *j);
+                        wg[UEdge(i,findVar(*j))] = piet.strength(v_i, *j);
                     }
             }
 
@@ -277,19 +279,21 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
     // Construct corresponding region graph
 
     // Create outer regions
-    ORs().reserve( Cliques.size() );
+    ORs.reserve( Cliques.size() );
     for( size_t i = 0; i < Cliques.size(); i++ )
-        ORs().push_back( FRegion( Factor(Cliques[i], 1.0), 1.0 ) );
+        ORs.push_back( FRegion( Factor(Cliques[i], 1.0), 1.0 ) );
 
     // For each factor, find an outer region that subsumes that factor.
     // Then, multiply the outer region with that factor.
     // If no outer region can be found subsuming that factor, label the
     // factor as off-tree.
+    fac2OR.clear();
+    fac2OR.resize( nrFactors(), -1U );
     for( size_t I = 0; I < nrFactors(); I++ ) {
         size_t alpha;
-        for( alpha = 0; alpha < nr_ORs(); alpha++ )
+        for( alpha = 0; alpha < nrORs(); alpha++ )
             if( OR(alpha).vars() >> factor(I).vars() ) {
-                _fac2OR[I] = alpha;
+                fac2OR[I] = alpha;
                 break;
             }
     // DIFF WITH JTree::GenerateJT:      assert
@@ -297,30 +301,32 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
     RecomputeORs();
 
     // Create inner regions and edges
-    IRs().reserve( _RTree.size() );
-    Redges().reserve( 2 * _RTree.size() );
+    IRs.reserve( _RTree.size() );
+    typedef pair<size_t,size_t> Edge;
+    vector<Edge> edges;
+    edges.reserve( 2 * _RTree.size() );
     for( size_t i = 0; i < _RTree.size(); i++ ) {
-        Redges().push_back( R_edge_t( _RTree[i].n1, IRs().size() ) );
-        Redges().push_back( R_edge_t( _RTree[i].n2, IRs().size() ) );
+        edges.push_back( Edge( _RTree[i].n1, IRs.size() ) );
+        edges.push_back( Edge( _RTree[i].n2, IRs.size() ) );
         // inner clusters have counting number -1
-        IRs().push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
+        IRs.push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
     }
 
-    // Regenerate BipartiteGraph internals
-    Regenerate();
+    // create bipartite graph
+    G.create( nrORs(), nrIRs(), edges.begin(), edges.end() );
 
     // Check counting numbers
     Check_Counting_Numbers();
     
     // Create messages and beliefs
     _Qa.clear();
-    _Qa.reserve( nr_ORs() );
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    _Qa.reserve( nrORs() );
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         _Qa.push_back( OR(alpha) );
 
     _Qb.clear();
-    _Qb.reserve( nr_IRs() );
-    for( size_t beta = 0; beta < nr_IRs(); beta++ ) 
+    _Qb.reserve( nrIRs() );
+    for( size_t beta = 0; beta < nrIRs(); beta++ ) 
         _Qb.push_back( Factor( IR(beta), 1.0 ) );
 
     // DIFF with JTree::GenerateJT:  no messages
@@ -459,13 +465,13 @@ Complex TreeEP::logZ() const {
     double sum = 0.0;
 
     // entropy of the tree
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         sum -= real(_Qb[beta].entropy());
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         sum += real(_Qa[alpha].entropy());
 
     // energy of the on-tree factors
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
 
     // energy of the off-tree factors
index dc88a21..ff29c84 100644 (file)
@@ -284,7 +284,7 @@ int main( int argc, char *argv[] ) {
                 throw "Please specify all required arguments";
             do {
                 MakeHOIFG( N, M, k, beta, fg );
-            } while( !fg.isConnected() );
+            } while( !fg.G.isConnected() );
 
             cout << "# N = " << N << endl;
             cout << "# M = " << M << endl;
index a505b7e..d50c776 100644 (file)
@@ -55,8 +55,9 @@ int main( int argc, char *argv[] ) {
                 cout << "node[shape=box,style=filled,color=lightgrey,width=0.3,height=0.3,fixedsize=true];" << endl;
                 for( size_t I = 0; I < fg.nrFactors(); I++ )
                     cout << "\tp" << I << ";" << endl;
-                for( size_t iI = 0; iI < fg.nr_edges(); iI++ )
-                    cout << "\tx" << fg.var(fg.edge(iI).first).label() << " -- p" << fg.edge(iI).second << ";" << endl;
+                for( size_t i = 0; i < fg.nrVars(); i++ )
+                    foreach( const FactorGraph::Neighbor &I, fg.nbV(i) )
+                        cout << "\tx" << fg.var(i).label() << " -- p" << I << ";" << endl;
                 cout << "}" << endl;
             }
 
index a335106..e957abd 100644 (file)
@@ -44,15 +44,15 @@ int main( int argc, char *argv[] ) {
         } else {
             cout << "Number of variables: " << fg.nrVars() << endl;
             cout << "Number of factors:   " << fg.nrFactors() << endl;
-            cout << "Connected:           " << fg.isConnected() << endl;
+            cout << "Connected:           " << fg.G.isConnected() << endl;
 //          cout << "Treewidth:           " << endl;
 
             double cavsum_lcbp = 0.0;
             double cavsum_lcbp2 = 0.0;
             size_t max_Delta_size = 0;
             for( size_t i = 0; i < fg.nrVars(); i++ ) {
-                VarSet di = fg.delta(fg.var(i));
-                size_t Ds = fg.Delta(fg.var(i)).stateSpace();
+                VarSet di = fg.delta(i);
+                size_t Ds = fg.Delta(i).stateSpace();
                 if( Ds > max_Delta_size )
                     max_Delta_size = Ds;
                 cavsum_lcbp += di.stateSpace();