Fixed nasty bug in BipartiteGraph::erase1 and erase2
authorJoris Mooij <jorism@marvin.jorismooij.nl>
Mon, 15 Sep 2008 10:20:12 +0000 (12:20 +0200)
committerJoris Mooij <jorism@marvin.jorismooij.nl>
Mon, 15 Sep 2008 10:20:12 +0000 (12:20 +0200)
- Fixed bug in example.cpp (logdomain option was missing)
- Moved some code from the BipartiteGraph header to a
  new source file src/bipgraph.cpp
- Added some consistency checks to BipartiteGraph

Makefile
example.cpp
include/dai/bipgraph.h
src/bipgraph.cpp [new file with mode: 0644]

index f469d67..154debc 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -111,8 +111,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 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)
+$(LIB)/libdai.a : bipgraph.o daialg.o alldai.o clustergraph.o factorgraph.o properties.o regiongraph.o util.o weightedgraph.o x2x.o $(OBJECTS)
+       ar rcs $(LIB)/libdai.a bipgraph.o daialg.o alldai.o clustergraph.o factorgraph.o properties.o regiongraph.o util.o weightedgraph.o x2x.o $(OBJECTS)
 
 tests : tests/test
 
@@ -129,6 +129,8 @@ clean :
        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
 
+bipgraph.o : $(SRC)/bipgraph.cpp $(HEADERS)
+       $(CC) $(CCFLAGS) -c $(SRC)/bipgraph.cpp
 
 daialg.o : $(SRC)/daialg.cpp $(HEADERS)
        $(CC) $(CCFLAGS) -c $(SRC)/daialg.cpp
@@ -216,8 +218,8 @@ matlab/matlab.o : matlab/matlab.cpp matlab/matlab.h $(HEADERS)
 # UTILS
 ########
 
-utils/createfg : utils/createfg.cpp $(HEADERS) factorgraph.o weightedgraph.o util.o
-       $(CC) $(CCFLAGS) -o utils/createfg utils/createfg.cpp factorgraph.o weightedgraph.o util.o $(BOOSTFLAGS)
+utils/createfg : utils/createfg.cpp $(HEADERS) factorgraph.o weightedgraph.o util.o bipgraph.o
+       $(CC) $(CCFLAGS) -o utils/createfg utils/createfg.cpp factorgraph.o weightedgraph.o util.o bipgraph.o $(BOOSTFLAGS)
 
 utils/fg2dot : utils/fg2dot.cpp $(HEADERS) factorgraph.o
        $(CC) $(CCFLAGS) -o utils/fg2dot utils/fg2dot.cpp factorgraph.o
@@ -225,5 +227,5 @@ utils/fg2dot : utils/fg2dot.cpp $(HEADERS) factorgraph.o
 utils/remove_short_loops : utils/remove_short_loops.cpp $(HEADERS) factorgraph.o
        $(CC) $(CCFLAGS) -o utils/remove_short_loops utils/remove_short_loops.cpp factorgraph.o
 
-utils/fginfo : utils/fginfo.cpp $(HEADERS) factorgraph.o
-       $(CC) $(CCFLAGS) -o utils/fginfo utils/fginfo.cpp factorgraph.o
+utils/fginfo : utils/fginfo.cpp $(HEADERS) factorgraph.o bipgraph.o
+       $(CC) $(CCFLAGS) -o utils/fginfo utils/fginfo.cpp factorgraph.o bipgraph.o
index 53adcfc..527bf70 100644 (file)
@@ -57,7 +57,7 @@ int main( int argc, char *argv[] ) {
             for( size_t i = 0; i < fg.nrVars(); i++ )
                 cout << jt.belief(fg.var(i)) << endl;
 
-            BP bp(fg, opts("updates",string("SEQMAX")));
+            BP bp(fg, opts("updates",string("SEQMAX"))("logdomain",false));
             bp.init();
             bp.run();
 
index d6e720e..c9701ef 100644 (file)
@@ -119,24 +119,68 @@ class BipartiteGraph {
         }
 
         /// Returns constant reference to the _i2'th neighbor of node i1 of type 1
-        const Neighbor & nb1( size_t i1, size_t _i2 ) const { return _nb1[i1][_i2]; }
+        const Neighbor & nb1( size_t i1, size_t _i2 ) const { 
+#ifdef DAI_DEBUG
+            assert( i1 < _nb1.size() );
+            assert( _i2 < _nb1[i1].size() );
+#endif
+            return _nb1[i1][_i2]; 
+        }
         /// Returns reference to the _i2'th neighbor of node i1 of type 1
-        Neighbor & nb1( size_t i1, size_t _i2 ) { return _nb1[i1][_i2]; }
+        Neighbor & nb1( size_t i1, size_t _i2 ) {
+#ifdef DAI_DEBUG
+            assert( i1 < _nb1.size() );
+            assert( _i2 < _nb1[i1].size() );
+#endif
+            return _nb1[i1][_i2]; 
+        }
 
         /// Returns constant reference to the _i1'th neighbor of node i2 of type 2
-        const Neighbor & nb2( size_t i2, size_t _i1 ) const { return _nb2[i2][_i1]; }
+        const Neighbor & nb2( size_t i2, size_t _i1 ) const { 
+#ifdef DAI_DEBUG
+            assert( i2 < _nb2.size() );
+            assert( _i1 < _nb2[i2].size() );
+#endif
+            return _nb2[i2][_i1]; 
+        }
         /// Returns reference to the _i1'th neighbor of node i2 of type 2
-        Neighbor & nb2( size_t i2, size_t _i1 ) { return _nb2[i2][_i1]; }
+        Neighbor & nb2( size_t i2, size_t _i1 ) { 
+#ifdef DAI_DEBUG
+            assert( i2 < _nb2.size() );
+            assert( _i1 < _nb2[i2].size() );
+#endif
+            return _nb2[i2][_i1]; 
+        }
 
         /// Returns constant reference to all neighbors of node i1 of type 1
-        const Neighbors & nb1( size_t i1 ) const { return _nb1[i1]; }
+        const Neighbors & nb1( size_t i1 ) const { 
+#ifdef DAI_DEBUG
+            assert( i1 < _nb1.size() );
+#endif
+            return _nb1[i1]; 
+        }
         /// Returns reference to all neighbors of node of i1 type 1
-        Neighbors & nb1( size_t i1 ) { return _nb1[i1]; }
+        Neighbors & nb1( size_t i1 ) { 
+#ifdef DAI_DEBUG
+            assert( i1 < _nb1.size() );
+#endif
+            return _nb1[i1]; 
+        }
 
         /// Returns constant reference to all neighbors of node i2 of type 2
-        const Neighbors & nb2( size_t i2 ) const { return _nb2[i2]; }
+        const Neighbors & nb2( size_t i2 ) const { 
+#ifdef DAI_DEBUG
+            assert( i2 < _nb2.size() );
+#endif
+            return _nb2[i2]; 
+        }
         /// Returns reference to all neighbors of node i2 of type 2
-        Neighbors & nb2( size_t i2 ) { return _nb2[i2]; }
+        Neighbors & nb2( size_t i2 ) { 
+#ifdef DAI_DEBUG
+            assert( i2 < _nb2.size() );
+#endif
+            return _nb2[i2]; 
+        }
 
         /// Returns number of nodes of type 1
         size_t nr1() const { return _nb1.size(); }
@@ -152,14 +196,10 @@ class BipartiteGraph {
         }
         
         /// Add node of type 1 without neighbors.
-        void add1() {
-            _nb1.push_back( Neighbors() );
-        }
+        void add1() { _nb1.push_back( Neighbors() ); }
         
         /// Add node of type 2 without neighbors.
-        void add2() {
-            _nb2.push_back( Neighbors() );
-        }
+        void add2() { _nb2.push_back( Neighbors() ); }
 
         /// Add node of type 1 with neighbors specified by a range.
         /** The value_type of an NodeInputIterator should be a size_t, corresponding to
@@ -202,222 +242,35 @@ class BipartiteGraph {
         }
 
         /// Remove node of type 1 and all incident edges.
-        void erase1( size_t n1 ) {
-            assert( n1 < nr1() );
-            // Erase neighbor entry of node n1
-            _nb1.erase( _nb1.begin() + n1 );
-            // Adjust neighbor entries of nodes of type 2
-            for( size_t n2 = 0; n2 < nr2(); n2++ )
-                for( size_t iter = 0; iter < nb2(n2).size(); ) {
-                    if( nb2(n2, iter).node == n1 ) {
-                        // delete this entry, because it points to the deleted node
-                        nb2(n2).erase( nb2(n2).begin() + iter );
-                        // adjust all subsequent entries:
-                        // update their iter and the corresponding dual of the neighboring node of type 1
-                        for( size_t newiter = iter; newiter < nb2(n2).size(); newiter++ ) {
-                            nb2( n2, newiter ).iter = newiter;
-                            nb1( nb2(n2, newiter).node, nb2(n2, newiter).dual ).dual = newiter;
-                        }
-                    } else if( nb2(n2, iter).node > n1 ) {
-                        nb2(n2, iter).node--;
-                        iter++;
-                    } else
-                        iter++;
-                }
-        }
+        void erase1( size_t n1 );
 
         /// Remove node of type 2 and all incident edges.
-        void erase2( size_t n2 ) {
-            assert( n2 < nr2() );
-            // Erase neighbor entry of node n2
-            _nb2.erase( _nb2.begin() + n2 );
-            // Adjust neighbor entries of nodes of type 1
-            for( size_t n1 = 0; n1 < nr1(); n1++ )
-                for( size_t iter = 0; iter < nb1(n1).size(); ) {
-                    if( nb1(n1, iter).node == n2 ) {
-                        // delete this entry, because it points to the deleted node
-                        nb1(n1).erase( nb1(n1).begin() + iter );
-                        // adjust all subsequent entries:
-                        // update their iter and the corresponding dual of the neighboring node of type 2
-                        for( size_t newiter = iter; newiter < nb1(n1).size(); newiter++ ) {
-                            nb1( n1, newiter ).iter = newiter;
-                            nb2( nb1(n1, newiter).node, nb1(n1, newiter).dual ).dual = newiter;
-                        }
-                    } else if( nb1(n1, iter).node > n2 ) {
-                        nb1(n1, iter).node--;
-                        iter++;
-                    } else
-                        iter++;
-                }
-        }
+        void erase2( size_t n2 );
 
         /// Calculate second-order neighbors (i.e., neighbors of neighbors) of node n1 of type 1.
         /** If include == true, include n1 itself, otherwise exclude n1.
          */
-        std::vector<size_t> delta1( size_t n1, bool include = false ) const {
-            std::vector<size_t> result;
-            foreach( const Neighbor &n2, nb1(n1) )
-                foreach( const Neighbor &m1, nb2(n2) )
-                    if( include || (m1 != n1) )
-                        result.push_back( m1 );
-            // remove duplicates
-            std::vector<size_t>::iterator it = std::unique( result.begin(), result.end() );
-            result.erase( it, result.end() );
-            return result;
-        }
+        std::vector<size_t> delta1( size_t n1, bool include = false ) const;
 
         /// Calculate second-order neighbors (i.e., neighbors of neighbors) of node n2 of type 2.
         /** If include == true, include n2 itself, otherwise exclude n2.
          */
-        std::vector<size_t> delta2( size_t n2, bool include = false ) const {
-            std::vector<size_t> result;
-            foreach( const Neighbor &n1, nb2(n2) )
-                foreach( const Neighbor &m2, nb1(n1) )
-                    if( include || (m2 != n2) )
-                        result.push_back( m2 );
-            // remove duplicates
-            std::vector<size_t>::iterator it = std::unique( result.begin(), result.end() );
-            result.erase( it, result.end() );
-            return result;
-        }
+        std::vector<size_t> delta2( size_t n2, bool include = false ) const;
 
         /// Returns true if the graph is connected
-        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 type 2, 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 type 1, 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;
-            }
-        }
+        bool isConnected() const;
 
         /// Returns true if the graph is a tree, i.e., if it is singly connected and connected.
         /** This is equivalent to whether for each pair of vertices in the graph, there exists
          *  a unique path in the graph that starts at the first and ends at the second vertex.
          */
-        bool isTree() const {
-            using namespace std;
-            vector<levelType> levels;
-
-            bool foundCycle = false;
-            size_t nr_1 = 0;
-            size_t nr_2 = 0;
-
-            if( nr1() == 0 || nr2() == 0 )
-                return true;
-            else {
-                levelType newLevel;
-                do {
-                    newLevel.ind1.clear();
-                    newLevel.ind2.clear();
-                    if( levels.size() == 0 ) {
-                        size_t n1 = 0;
-                        // add n1 to ind1
-                        newLevel.ind1 = vector<size_t>( 1, n1 );
-                        // add all neighbors of n1 to ind2
-                        newLevel.ind2.reserve( nb1(n1).size() );
-                        foreach( const Neighbor &n2, nb1(n1) )
-                            newLevel.ind2.push_back( n2 );
-                    } else {
-                        const levelType &prevLevel = levels.back();
-                        // build newLevel.ind1
-                        foreach( size_t n2, prevLevel.ind2 ) { // for all n2 in the previous level
-                            foreach( const Neighbor &n1, nb2(n2) ) { // for all neighbors n1 of n2
-                                if( find( prevLevel.ind1.begin(), prevLevel.ind1.end(), n1 ) == prevLevel.ind1.end() ) { // n1 not in previous level
-                                    if( find( newLevel.ind1.begin(), newLevel.ind1.end(), n1 ) != newLevel.ind1.end() )
-                                        foundCycle = true; // n1 already in new level: we found a cycle
-                                    else
-                                        newLevel.ind1.push_back( n1 ); // add n1 to new level
-                                }
-                                if( foundCycle )
-                                    break;
-                            }
-                            if( foundCycle )
-                                break;
-                        }
-                        // build newLevel.ind2
-                        foreach( size_t n1, newLevel.ind1 ) { // for all n1 in this level
-                            foreach( const Neighbor &n2, nb1(n1) ) { // for all neighbors n2 of n1
-                                if( find( prevLevel.ind2.begin(), prevLevel.ind2.end(), n2 ) == prevLevel.ind2.end() ) { // n2 not in previous level
-                                    if( find( newLevel.ind2.begin(), newLevel.ind2.end(), n2 ) != newLevel.ind2.end() )
-                                        foundCycle = true; // n2 already in new level: we found a cycle
-                                    else
-                                        newLevel.ind2.push_back( n2 ); // add n2 to new level
-                                }
-                                if( foundCycle )
-                                    break;
-                            }
-                            if( foundCycle )
-                                break;
-                        } 
-                    }
-                    levels.push_back( newLevel );
-                    nr_1 += newLevel.ind1.size();
-                    nr_2 += newLevel.ind2.size();
-                } while( ((newLevel.ind1.size() != 0) || (newLevel.ind2.size() != 0)) && !foundCycle );
-                if( nr_1 == nr1() && nr_2 == nr2() && !foundCycle )
-                    return true;
-                else
-                    return false;
-            }
-        }
+        bool isTree() const;
 
         /// Stream to output stream os in graphviz .dot syntax
-        void display( std::ostream& os ) const {
-            using namespace std;
-            os << "graph G {" << endl;
-            os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
-            for( size_t n1 = 0; n1 < nr1(); n1++ )
-                os << "\tx" << n1 << ";" << endl;
-            os << "node[shape=box,width=0.3,height=0.3,fixedsize=true];" << endl;
-            for( size_t n2 = 0; n2 < nr2(); n2++ )
-                os << "\ty" << n2 << ";" << endl;
-            for( size_t n1 = 0; n1 < nr1(); n1++ )
-                foreach( const Neighbor &n2, nb1(n1) )
-                    os << "\tx" << n1 << " -- y" << n2 << ";" << endl;
-            os << "}" << endl;
-        }
+        void display( std::ostream& os ) const;
+
+        /// Checks internal consistency
+        void check() const;
 };
 
 
diff --git a/src/bipgraph.cpp b/src/bipgraph.cpp
new file mode 100644 (file)
index 0000000..1e87de1
--- /dev/null
@@ -0,0 +1,285 @@
+/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
+    Radboud University Nijmegen, The Netherlands
+    
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+#include <dai/bipgraph.h>
+
+
+namespace dai {
+
+
+using namespace std;
+
+
+/// Remove node of type 1 and all incident edges.
+void BipartiteGraph::erase1( size_t n1 ) {
+    assert( n1 < nr1() );
+    // Erase neighbor entry of node n1
+    _nb1.erase( _nb1.begin() + n1 );
+    // Adjust neighbor entries of nodes of type 2
+    for( size_t n2 = 0; n2 < nr2(); n2++ ) {
+        for( size_t iter = 0; iter < nb2(n2).size(); ) {
+            Neighbor &m1 = nb2(n2, iter);
+            if( m1.node == n1 ) {
+                // delete this entry, because it points to the deleted node
+                nb2(n2).erase( nb2(n2).begin() + iter );
+            } else if( m1.node > n1 ) {
+                // update this entry and the corresponding dual of the neighboring node of type 1
+                m1.iter = iter;
+                m1.node--;
+                nb1( m1.node, m1.dual ).dual = iter;
+                iter++;
+            } else {
+                // skip
+                iter++;
+            }
+        }
+    }
+}
+
+
+/// Remove node of type 2 and all incident edges.
+void BipartiteGraph::erase2( size_t n2 ) {
+    assert( n2 < nr2() );
+    // Erase neighbor entry of node n2
+    _nb2.erase( _nb2.begin() + n2 );
+    // Adjust neighbor entries of nodes of type 1
+    for( size_t n1 = 0; n1 < nr1(); n1++ ) {
+        for( size_t iter = 0; iter < nb1(n1).size(); ) {
+            Neighbor &m2 = nb1(n1, iter);
+            if( m2.node == n2 ) {
+                // delete this entry, because it points to the deleted node
+                nb1(n1).erase( nb1(n1).begin() + iter );
+            } else if( m2.node > n2 ) {
+                // update this entry and the corresponding dual of the neighboring node of type 2
+                m2.iter = iter;
+                m2.node--;
+                nb2( m2.node, m2.dual ).dual = iter;
+                iter++;
+            } else {
+                // skip
+                iter++;
+            }
+        }
+    }
+}
+
+
+/// Calculate second-order neighbors (i.e., neighbors of neighbors) of node n1 of type 1.
+/** If include == true, include n1 itself, otherwise exclude n1.
+ */
+std::vector<size_t> BipartiteGraph::delta1( size_t n1, bool include ) const {
+    std::vector<size_t> result;
+    foreach( const Neighbor &n2, nb1(n1) )
+        foreach( const Neighbor &m1, nb2(n2) )
+            if( include || (m1 != n1) )
+                result.push_back( m1 );
+    // remove duplicates
+    std::vector<size_t>::iterator it = std::unique( result.begin(), result.end() );
+    result.erase( it, result.end() );
+    return result;
+}
+
+
+/// Calculate second-order neighbors (i.e., neighbors of neighbors) of node n2 of type 2.
+/** If include == true, include n2 itself, otherwise exclude n2.
+ */
+std::vector<size_t> BipartiteGraph::delta2( size_t n2, bool include ) const {
+    std::vector<size_t> result;
+    foreach( const Neighbor &n1, nb2(n2) )
+        foreach( const Neighbor &m2, nb1(n1) )
+            if( include || (m2 != n2) )
+                result.push_back( m2 );
+    // remove duplicates
+    std::vector<size_t>::iterator it = std::unique( result.begin(), result.end() );
+    result.erase( it, result.end() );
+    return result;
+}
+
+
+/// Returns true if the graph is connected
+bool BipartiteGraph::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 type 2, 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 type 1, 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;
+    }
+}
+
+
+/// Returns true if the graph is a tree, i.e., if it is singly connected and connected.
+/** This is equivalent to whether for each pair of vertices in the graph, there exists
+ *  a unique path in the graph that starts at the first and ends at the second vertex.
+ */
+bool BipartiteGraph::isTree() const {
+    using namespace std;
+    vector<levelType> levels;
+
+    bool foundCycle = false;
+    size_t nr_1 = 0;
+    size_t nr_2 = 0;
+
+    if( nr1() == 0 || nr2() == 0 )
+        return true;
+    else {
+        levelType newLevel;
+        do {
+            newLevel.ind1.clear();
+            newLevel.ind2.clear();
+            if( levels.size() == 0 ) {
+                size_t n1 = 0;
+                // add n1 to ind1
+                newLevel.ind1 = vector<size_t>( 1, n1 );
+                // add all neighbors of n1 to ind2
+                newLevel.ind2.reserve( nb1(n1).size() );
+                foreach( const Neighbor &n2, nb1(n1) )
+                    newLevel.ind2.push_back( n2 );
+            } else {
+                const levelType &prevLevel = levels.back();
+                // build newLevel.ind1
+                foreach( size_t n2, prevLevel.ind2 ) { // for all n2 in the previous level
+                    foreach( const Neighbor &n1, nb2(n2) ) { // for all neighbors n1 of n2
+                        if( find( prevLevel.ind1.begin(), prevLevel.ind1.end(), n1 ) == prevLevel.ind1.end() ) { // n1 not in previous level
+                            if( find( newLevel.ind1.begin(), newLevel.ind1.end(), n1 ) != newLevel.ind1.end() )
+                                foundCycle = true; // n1 already in new level: we found a cycle
+                            else
+                                newLevel.ind1.push_back( n1 ); // add n1 to new level
+                        }
+                        if( foundCycle )
+                            break;
+                    }
+                    if( foundCycle )
+                        break;
+                }
+                // build newLevel.ind2
+                foreach( size_t n1, newLevel.ind1 ) { // for all n1 in this level
+                    foreach( const Neighbor &n2, nb1(n1) ) { // for all neighbors n2 of n1
+                        if( find( prevLevel.ind2.begin(), prevLevel.ind2.end(), n2 ) == prevLevel.ind2.end() ) { // n2 not in previous level
+                            if( find( newLevel.ind2.begin(), newLevel.ind2.end(), n2 ) != newLevel.ind2.end() )
+                                foundCycle = true; // n2 already in new level: we found a cycle
+                            else
+                                newLevel.ind2.push_back( n2 ); // add n2 to new level
+                        }
+                        if( foundCycle )
+                            break;
+                    }
+                    if( foundCycle )
+                        break;
+                } 
+            }
+            levels.push_back( newLevel );
+            nr_1 += newLevel.ind1.size();
+            nr_2 += newLevel.ind2.size();
+        } while( ((newLevel.ind1.size() != 0) || (newLevel.ind2.size() != 0)) && !foundCycle );
+        if( nr_1 == nr1() && nr_2 == nr2() && !foundCycle )
+            return true;
+        else
+            return false;
+    }
+}
+
+
+/// Stream to output stream os in graphviz .dot syntax
+void BipartiteGraph::display( std::ostream& os ) const {
+    using namespace std;
+    os << "graph G {" << endl;
+    os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
+    for( size_t n1 = 0; n1 < nr1(); n1++ )
+        os << "\tx" << n1 << ";" << endl;
+    os << "node[shape=box,width=0.3,height=0.3,fixedsize=true];" << endl;
+    for( size_t n2 = 0; n2 < nr2(); n2++ )
+        os << "\ty" << n2 << ";" << endl;
+    for( size_t n1 = 0; n1 < nr1(); n1++ )
+        foreach( const Neighbor &n2, nb1(n1) )
+            os << "\tx" << n1 << " -- y" << n2 << ";" << endl;
+    os << "}" << endl;
+}
+
+
+/// Checks internal consistency
+void BipartiteGraph::check() const {
+    size_t N1 = nr1();
+    size_t N2 = nr2();
+    for( size_t n1 = 0; n1 < N1; n1++ ) {
+        size_t iter = 0;
+        foreach( const Neighbor &n2, nb1(n1) ) {
+            assert( n2.iter == iter );
+            assert( n2.node < N2 );
+            assert( n2.dual < nb2(n2).size() );
+            assert( nb2(n2, n2.dual) == n1 );
+            iter++;
+        }
+    }
+    for( size_t n2 = 0; n2 < N2; n2++ ) {
+        size_t iter = 0;
+        foreach( const Neighbor &n1, nb2(n2) ) {
+            assert( n1.iter == iter );
+            assert( n1.node < N1 );
+            assert( n1.dual < nb1(n1).size() );
+            assert( nb1(n1, n1.dual) == n2 );
+            iter++;
+        }
+    }
+}
+
+
+} // end of namespace dai