Context Navigation

Back to Ticket #168

Ticket #168: bpm_upgrade.2.patch

File bpm_upgrade.2.patch, 81.2 KB (added by Poroszkai Daniel, 12 years ago)
Improved Hopcroft-Karp; two preflow based matching

lemon/hopcroft_karp.h

# HG changeset patch
# User Daniel Poroszkai <poroszd@inf.elte.hu>
# Date 1330129198 -3600
# Node ID 6cd8dd8bd86a11f12ed8c617e6af6d5256026603
# Parent  4928e7fa0a4062b96c9fada8abd40d3c7041bcae
Hopcroft-Karp improved

diff --git a/lemon/hopcroft_karp.h b/lemon/hopcroft_karp.h

-                      a
+ *
  * This file is a part of LEMON, a generic C++ optimization library.
+ *
  * Copyright (C) 2003-2011
+ * Copyright (C) 2003-2012
  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
  * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
 …
 #define LEMON_HOPCROFT_KARP_H
 #include <lemon/core.h>
+#include <vector>
+#include <queue>
 #include <list>
 /// \ingroup matching
 …
   template <typename BPG>
   class HopcroftKarp {
   public:
       /// Type of the bipartite graph
+    /// Type of the bipartite graph
     typedef BPG BpGraph;
     /// Type of the matching map
     typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
 …
     /// \brief Sets the matching map
     ///
+    /// Sets the matching map.
+    /// Sets the matching map. It should contain a valid matching,
+    /// for example the empty matching is fine.
     /// If you don't use this function before calling \ref init(),
     /// an instance will be allocated automatically.
     /// The destructor deallocates this automatically allocated map,
     /// of course.
-    /// This member is not to initialize the algorithm with a valid
-    /// matching; use \ref matchingInit() instead.
     ///
     /// \return <tt>(*this)</tt>
     HopcroftKarp& matchingMap(MatchingMap& map) {
 …
     /// \brief Initializes the algorithm
     ///
+    /// Allocates the matching map if necessary, and sets
+    /// to empty matching.
+    /// Allocates the matching map if necessary.
     ///
     /// \pre \ref init() is not called.
     void init() {
       createStructures();
-      if (!_local_matching) {
-        for (NodeIt it(_bpg); it!=INVALID; ++it) {
-          _matching->set(it, INVALID);
+        }
+      }
+    }
     /// \brief Smarter initialization of the matching.
 …
     ///
     /// \note heuristicInit() replaces init() regarding the preconditions.
     int heuristicInit() {
       init();
+      createStructures();
       int size = 0;
+      for (BlueNodeIt b_it(_bpg); b_it!=INVALID; ++b_it) {
+        if ((*_matching)[b_it] == INVALID) {
+          bool matched = false;
+          for (IncEdgeIt inc(_bpg, b_it); inc != INVALID && !matched; ++inc) {
+            if ((*_matching)[_bpg.redNode(inc)] == INVALID) {
+              _matching->set(_bpg.redNode(inc), inc);
+              _matching->set(b_it, inc);
+              matched = true;
+              ++size;
+            }
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        bool matched = false;
+        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
+          if ((*_matching)[_bpg.target(oi)] == INVALID) {
+            _matching->set(_bpg.source(oi), oi);
+            _matching->set(_bpg.target(oi), oi);
+            matched = true;
+            ++size;
+          }
+        }
+      }
       return size;
+    }
 …
     /// in which an edge is true if it is in the matching.
     ///
     /// If the given matching is invalid, some edges will be left out.
     ///
+    ///
     /// \return \c false if the given matching is invalid.
     ///
     /// \note matchingInit() replaces init() regarding the preconditions.
     bool matchingInit(const BoolEdgeMap& matching) {
       init();
+      createStructures();
       bool valid = true;
       for(EdgeIt it(_bpg); it!=INVALID; ++it) {
         if (matching[it]) {
 …
     ///
     /// \pre \ref init() must be called before using this function.
     int augment() {
+      IntNodeMap level(_bpg, -1);
+      std::list<RedNode> act_rednodes;
+      std::list<BlueNode> act_bluenodes;
+      IntBlueNodeMap level(_bpg, -2);
+      std::vector<RedNode> act_red_nodes;
       for (RedNodeIt r_it(_bpg); r_it != INVALID; ++r_it) {
         if ((*_matching)[r_it] == INVALID) {
           act_rednodes.push_front(r_it);
+          act_red_nodes.push_back(r_it);
+        }
+      }
+      // Raise this flag when a shortest augmenting path is found.
+      bool path_found = false;
+      // This nodelist will contain the end nodes of the possible
+      // This vector will contain the end nodes of the possible
       // augmenting paths.
       std::list<BlueNode> path_ends;
+      std::vector<BlueNode> path_ends;
       // Layer counter
       int cur_level = 0;
+      int cur_level = -1;
       // Starting from the unmatched red nodes search for unmatched
       // blue nodes, using Bfs (but only on alternating paths).
+      while (!path_found) {
+        while (!act_rednodes.empty()) {
+          RedNode red = act_rednodes.front();
+          act_rednodes.pop_front();
+          level[red] = cur_level;
+          for (IncEdgeIt it(_bpg, red); it!=INVALID; ++it) {
+            BlueNode blue(_bpg.blueNode(it));
+            if (level[blue] == -1) {
+              act_bluenodes.push_front(blue);
+              level[blue] = cur_level + 1;
+              path_found |= ((*_matching)[blue] == INVALID);
+      // Blue nodes with level==-2 are unreached.
+      while (path_ends.empty()) {
+        cur_level += 2;
+        std::vector<RedNode> next_red_nodes;
+        for (typename std::vector<RedNode>::iterator r=act_red_nodes.begin();
+             r != act_red_nodes.end(); ++r) {
+          for (OutArcIt it(_bpg, *r); it!=INVALID; ++it) {
+            BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+            if (level[blue] != -2) continue;
+            level[blue] = cur_level;
+            if ((*_matching)[blue] == INVALID) {
+              path_ends.push_back(blue);
+            } else {
+              next_red_nodes
+                .push_back(_bpg.redNode((*_matching)[blue]));
+            }
+          }
+        }
+        if (!path_found) {
+          if (act_bluenodes.empty()) return -1;
+          while (!act_bluenodes.empty()) {
+            BlueNode blue = act_bluenodes.front();
+            act_bluenodes.pop_front();
+            RedNode red = _bpg.redNode((*_matching)[blue]);
+            act_rednodes.push_front(red);
+          }
+        } else {
+          for (typename std::list<BlueNode>::iterator it=act_bluenodes.begin();
+               it != act_bluenodes.end();
+               ++it) {
+            if ((*_matching)[*it] == INVALID) path_ends.push_front(*it);
+          }
+          // Now path_ends contains nodes that are ends of
+          // shortest alternating paths.
+        }
+        cur_level += 2;
+        if (path_ends.empty() && next_red_nodes.empty())
+          return -1;
+        act_red_nodes.swap(next_red_nodes);
+        next_red_nodes.clear();
+      }
+      // The current_edge structure assures that edges iterated at most once
+      typename BpGraph::template BlueNodeMap<IncEdgeIt*> current_edge(_bpg, 0);
+      // length of found augmenting paths
+      int max_level = cur_level;
       // Stack for the Dfs
+      std::list<BlueNode> stack;
+      std::vector<OutArcIt> stack;
+      // Raise this flag when a shortest augmenting path is found.
+      bool path_found;
       // Searching backward, starting from the previously found
       // blue nodes, we build vertex disjoint alternating paths.
+      typename std::list<BlueNode>::iterator n_it = path_ends.begin();
+      // Blue nodes with level!=-2 are unreached.
+      typename std::vector<BlueNode>::iterator end = path_ends.begin();
+      while (n_it != path_ends.end()) {
+        stack.push_front(*n_it);
+      while (end != path_ends.end()) {
+        stack.clear();
+        stack.push_back(OutArcIt(_bpg, *end));
+        cur_level = max_level;
         path_found = false;
         while (!stack.empty() && !path_found) {
+          BlueNode b = stack.front();
+          if (current_edge[b] == 0) {
+            current_edge[b] = new IncEdgeIt(_bpg, b);
+          OutArcIt &o = stack.back();
+          // Search an arc that goes one level higher.
+          // If none found, retreat.
+          while (o != INVALID &&
+                 (*_matching)[_bpg.target(o)] != INVALID &&
+                 level[_bpg.blueNode((*_matching)[_bpg.target(o)])] !=
+                 cur_level - 2) {
+            ++o;
+          }
+          if (o==INVALID) {
+            stack.pop_back();
+            cur_level += 2;
+          } else if ((*_matching)[_bpg.target(o)] == INVALID) {
+            path_found = true;
           } else {
+            ++(*current_edge[b]);
+          }
+          while (*current_edge[b] != INVALID &&
+                 level[_bpg.redNode(*current_edge[b])] != level[b] - 1) {
+            ++(*current_edge[b]);
+          }
+          if (*current_edge[b] == INVALID) {
+            stack.pop_front();
+          } else {
+            RedNode r = _bpg.redNode(*current_edge[b]);
+            if ((*_matching)[r] == INVALID) {
+              path_found = true;
+            } else {
+              level[r] = -1; // Do not visit this node again
+              stack.push_front(_bpg.blueNode((*_matching)[r]));
+            }
+            BlueNode b = _bpg.blueNode((*_matching)[_bpg.target(o)]);
+            stack.push_back(OutArcIt(_bpg, b));
+            level[b] = -2;        // Do not visit this node again
+            cur_level -= 2;
+          }
+        }
         if (path_found) {
+          BlueNode next = *n_it;
+          RedNode r;
+          BlueNode b;
+          Edge e;
+          while (next != INVALID) {
+            e = *current_edge[next];
+            b = next;
+            r = _bpg.redNode(e);
+            next = (*_matching)[r] == INVALID ?
+              INVALID :_bpg.blueNode((*_matching)[r]);
+            _matching->set(b, e);
+            _matching->set(r, e);
+          OutArcIt o;
+          while (!stack.empty()) {
+            o = stack.back();
+            _matching->set(_bpg.source(o), o);
+            _matching->set(_bpg.target(o), o);
+            stack.pop_back();
+          }
+        }
+        stack.clear();
+        ++n_it;
+        ++end;
+      }
       return cur_level - 1;
+      return max_level;
+    }
     /// \brief Executes the algorithm
 …
     /// \brief Size of the matching
     ///
     /// Returns the size of the current matching.
+    /// Function runs in \f$O(n)\f$ time.
     int matchingSize() const {
       int size = 0;
       for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
 …
       return (*_matching)[n] != INVALID ?
         _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
+    }
+    /// \brief Query a vertex cover.
+    ///
+    /// This function finds a vertex cover based on the current matching,
+    /// and in \c cover sets node in the cover \c true and \c false
+    /// otherwise. In case of maximal matching, this will be a minimal
+    /// vertex cover.
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return The number of nodes in the cover.
+    int cover(BoolNodeMap &cover) const {
+      int count = 0;
+      std::queue<RedNode> q;
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
+        cover[b] = false;
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          cover[r] = false;
+          q.push(r);
+        } else {
+          cover[r] = true;
+          ++count;
+        }
+      }
+      while (!q.empty()) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (cover[blue]) continue;
+          cover[blue] = true;
+          ++count;
+          if ((*_matching)[blue] != INVALID) {
+            RedNode nr = _bpg.redNode((*_matching)[blue]);
+            q.push(nr);
+            cover[nr] = false;
+            --count;
+          }
+        }
+        q.pop();
+      }
+      return count;
+    }
+    /// \brief Query a barrier of red nodes.
+    ///
+    /// This function tries to create a barrier of red nodes, which should:
+    ///  - contain every unmatched red nodes,
+    ///  - contain exactly that many matched red nodes as much blue nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool redBarrier(BoolRedNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<RedNode> q;
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          barrier[r] = true;
+          q.push(r);
+        } else {
+          barrier[r] = false;
+        }
+      }
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (exist &= ((*_matching)[blue] != INVALID)) {
+            RedNode r = _bpg.redNode((*_matching)[blue]);
+            if (!barrier[r]) {
+              q.push(r);
+              barrier[r] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+      return exist;
+    }
+    /// \brief Query a barrier of blue nodes.
+    ///
+    /// This function tries to create a barrier of blue nodes, which should:
+    ///  - contain every unmatched blue nodes,
+    ///  - contain exactly that many matched blue nodes as much red nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool blueBarrier(BoolBlueNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<BlueNode> q;
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        if ((*_matching)[b] == INVALID) {
+          barrier[b] = true;
+          q.push(b);
+        } else {
+          barrier[b] = false;
+        }
+      }
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
+          if (exist &= ((*_matching)[r] != INVALID)) {
+            BlueNode b = _bpg.blueNode((*_matching)[r]);
+            if (!barrier[b]) {
+              q.push(b);
+              barrier[b] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+      return exist;
+    }
   };
+}

test/hopcroft_karp_test.cc

diff --git a/test/hopcroft_karp_test.cc b/test/hopcroft_karp_test.cc

-                      a
+ *
  * This file is a part of LEMON, a generic C++ optimization library.
+ *
  * Copyright (C) 2003-2011
+ * Copyright (C) 2003-2012
  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
  * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
 …
   BoolEdgeMap init_matching(graph);
   HK hk_test(graph);
   const HK& const_hk_test = hk_test;
+  HK::MatchingMap matching_map(graph);
+  HK::MatchingMap matching_map(graph, INVALID);
+  BoolNodeMap cov(graph);
+  int cover_size;
+  BoolRedNodeMap rb(graph);
+  BoolBlueNodeMap bb(graph);
   hk_test.matchingMap(matching_map);
   hk_test.init();
 …
   const_hk_test.mate(node);
   const HK::MatchingMap& max_matching = const_hk_test.matchingMap();
   edge = max_matching[node];
+  cover_size = const_hk_test.cover(cov);
+  const_hk_test.redBarrier(rb);
+  const_hk_test.blueBarrier(bb);
+}
 typedef SmartBpGraph BpGraph;
 …
+}
+void checkMatchingOptimality(const BpGraph& bpg, const HK& hk) {
+  list<RedNode> reds;
+  list<BlueNode> blues;
+  BoolBlueNodeMap reached(bpg, false);
+void checkCover(const BpGraph& bpg, const HK& hk) {
+  int cov_size;
+  BoolNodeMap cov(bpg);
+  cov_size = hk.cover(cov);
+  BoolEdgeMap covered(bpg, false);
+  for (NodeIt n(bpg); n!=INVALID; ++n) {
+    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
+      covered[e] = true;
+    }
+  }
+  bool all = true;
+  for (EdgeIt e(bpg); e!=INVALID; ++e) {
+    all&=covered[e];
+  }
+  for (RedNodeIt it(bpg); it != INVALID; ++it) {
+    if (hk.matching(it) == INVALID) {
+      reds.push_front(it);
+  check(all, "Invalid cover");
+  check(hk.matchingSize() == cov_size, "Suboptimal matching.");
+}
+void checkBarriers(const BpGraph& bpg, const HK& hk) {
+  BoolRedNodeMap rb(bpg);
+  BoolBlueNodeMap bb(bpg);
+  check(hk.redBarrier(rb), "Suboptimal matching.");
+  check(hk.blueBarrier(bb), "Suboptimal matching.");
+  BoolNodeMap reached(bpg, false);
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (!rb[r]) continue;
+    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (!bb[b]) continue;
+    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  while (!reds.empty()) {
+    while (!reds.empty()) {
+      RedNode red = reds.front();
+      reds.pop_front();
+      for (IncEdgeIt it(bpg, red); it != INVALID; ++it) {
+        BlueNode blue = bpg.blueNode(it);
+        if (!reached[blue]) {
+          blues.push_front(blue);
+          reached.set(blue, true);
+        }
+      }
+    }
+    while (!blues.empty()) {
+      BlueNode blue = blues.front();
+      blues.pop_front();
+      check(hk.matching(blue) != INVALID, "Suboptimal matching");
+      reds.push_front(bpg.asRedNode(hk.mate(blue)));
+    }
+  int ur = 0, ub = 0, rc = 0, bc = 0;
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (rb[r]) ++rc;
+    if (hk.mate(r) == INVALID) ++ur;
+    if (reached[r]) --bc;
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (bb[b]) ++bc;
+    if (hk.mate(b) == INVALID) ++ub;
+    if (reached[b]) --rc;
+  }
+  check(rc == ur, "Red barrier is wrong.");
+  check(bc == ub, "Blue barrier is wrong.");
+}
+void checkMatching(const BpGraph& bpg, const HK& hk) {
+  checkMatchingValidity(bpg, hk);
+  checkBarriers(bpg, hk);
+  checkCover(bpg, hk);
+}
 int main() {
   // Checking hard wired extremal graphs
   for (int i=0; i<lgfn; ++i) {
 …
     bpGraphReader(bpg, lgfs).run();
     HK hk(bpg);
     hk.run();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+    checkMatching(bpg, hk);
+  }
   // Checking some use cases
 …
     hk.init();
     hk.augment();
     hk.start();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+    checkMatching(bpg, hk);
+  }
+  {
     HK hk(bpg);
     hk.heuristicInit();
     hk.augment();
     hk.start();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+    checkMatching(bpg, hk);
+  }
+  {
     HK hk(bpg);
     HK::MatchingMap m(bpg);
+    HK::MatchingMap m(bpg, INVALID);
     hk.matchingMap(m);
     hk.run();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+    checkMatching(bpg, hk);
+  }
+  {
     HK hk(bpg);
 …
     init_matching[e] = true;
     check(hk.matchingInit(init_matching), "Wrong result from matchingInit()");
     hk.start();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+    checkMatching(bpg, hk);
+  }
+  {
     HK hk(bpg);
     HK::MatchingMap m(bpg);
+    HK::MatchingMap m(bpg, INVALID);
     hk.matchingMap(m);
     BpGraph::EdgeMap<bool> init_matching(bpg, true);
     check(!hk.matchingInit(init_matching), "Wrong result from matchingInit()");
     hk.start();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+    checkMatching(bpg, hk);
+  }
   // Checking some random generated graphs
   const int random_test_n = 10;
   const int max_red_n = 100;
+  const int max_red_n = 1000;
   const int min_red_n = 10;
   const int max_blue_n = 100;
+  const int max_blue_n = 1000;
   const int min_blue_n = 10;
   for (int i=0; i<random_test_n; ++i) {
     int red_n = rnd.integer(min_red_n, max_red_n);
 …
+    }
     for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
       for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
         if (rnd() < prob) {
+        if (rnd() < prob/10.) {
           bpg.addEdge(r,b);
+        }
+      }
 …
     HK hk(bpg);
     hk.run();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+    checkMatching(bpg, hk);
+  }

lemon/elevator.h

# HG changeset patch
# User Daniel Poroszkai <poroszd@inf.elte.hu>
# Date 1329865937 -3600
# Node ID b5f7d8e60bc7f599bf0e2bf60bda1e05cf3c6d47
# Parent  dbf9516f49c509e59a9c73da18eb868cfafb751c
Fix initialization bug in elevator

diff --git a/lemon/elevator.h b/lemon/elevator.h

-                      a
+ *
  * This file is a part of LEMON, a generic C++ optimization library.
+ *
  * Copyright (C) 2003-2009
+ * Copyright (C) 2003-2012
  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
  * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
 …
     Elevator(const GR &graph,int max_level) :
       _g(graph),
       _max_level(max_level),
       _item_num(_max_level),
+      _item_num(countItems<GR, Item>(graph)),
       _where(graph),
       _level(graph,0),
       _items(_max_level),
+      _items(_item_num),
       _first(_max_level+2),
       _last_active(_max_level+2),
       _highest_active(-1) {}
 …
     ///\param max_level The maximum allowed level.
     ///Set the range of the possible labels to <tt>[0..max_level]</tt>.
     LinkedElevator(const GR& graph, int max_level)
+      : _graph(graph), _max_level(max_level), _item_num(_max_level),
+      : _graph(graph), _max_level(max_level),
+        _item_num(countItems<GR, Item>(graph)),
         _first(_max_level + 1), _last(_max_level + 1),
         _prev(graph), _next(graph),
         _highest_active(-1), _level(graph), _active(graph) {}

new file lemon/preflow_bp_matching.h

# HG changeset patch
# User Daniel Poroszkai <poroszd@inf.elte.hu>
# Date 1330129742 -3600
# Node ID 29f2b9488e4ea58b107564fa272ae38b402b2c1e
# Parent  b5f7d8e60bc7f599bf0e2bf60bda1e05cf3c6d47
Preflow based bipartite matching

diff --git a/lemon/preflow_bp_matching.h b/lemon/preflow_bp_matching.h
new file mode 100644

-                      -
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2012
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+#ifndef LEMON_PREFLOW_BP_MATCHING_H
+#define LEMON_PREFLOW_BP_MATCHING_H
+#include <lemon/core.h>
+#include <lemon/elevator.h>
+#include <vector>
+#include <list>
+#include <queue>
+/// \ingroup matching
+/// \file
+/// \brief A preflow based bipartite matching algorithm.
+namespace lemon {
+  /// \brief A preflow based bipartite matching algorithm
+  ///
+  /// Finds maximum cardinality matching in a given bipartite
+  /// graph using specialized preflow algorithm.
+  template <typename BPG>
+  class PreflowBpMatching {
+  public:
+    /// Type of the bipartite graph
+    typedef BPG BpGraph;
+    /// Type of the matching map
+    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
+    /// Type of the elevator.
+    typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t;
+  private:
+    TEMPLATE_BPGRAPH_TYPEDEFS(BPG);
+  protected:
+    typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow;
+    typedef std::list<typename BPG::Arc> CurrentArcSet;
+    typedef typename
+    std::list<typename BPG::Arc>::iterator CurrentArcIterator;
+    const BpGraph& _bpg;
+    MatchingMap* _matching;
+    bool _local_matching;
+    Elevator_t* _elevator;
+    bool _local_elevator;
+    typename BPG::template BlueNodeMap<CurrentArcSet> _cas;
+    typename BPG::template BlueNodeMap<CurrentArcIterator> _ca;
+    Preflow _flow;
+    void createStructures() {
+      if (!_matching) {
+        _matching = new MatchingMap(_bpg, INVALID);
+        _local_matching = true;
+      }
+      if (!_elevator) {
+        _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1,
+                                                  countBlueNodes(_bpg)));
+        _local_elevator = true;
+      }
+    }
+    void destroyStructures() {
+      if (_local_matching) {
+        delete _matching;
+      }
+      if (_local_elevator) {
+        delete _elevator;
+      }
+    }
+    /// \brief Initialize the elevator
+    ///
+    /// When this function is called, the elevator is initialized with
+    /// the exact distance labels corresponding the current matching.
+    void initElevator() {
+      // Set the initial flow:
+      // choose arbitrary outgoing arc for every unmatched red node.
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        Arc a = (*_matching)[r] != INVALID ?
+          _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r);
+        if (a != INVALID) {
+          _flow[r] = a;
+          _cas[_bpg.asBlueNodeUnsafe(_bpg.target(a))].push_front(a);
+        }
+      }
+      // Now we start a backward Bfs from the unmatched blue nodes in the
+      // residual graph, to supply the exact distance labels to the elevator.
+      std::vector<BlueNode> act_blue_nodes, next_blue_nodes;
+      BoolBlueNodeMap reached(_bpg, false);
+      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
+        if (_cas[b].empty()) {
+          act_blue_nodes.push_back(b);
+          reached.set(b, true);
+        }
+      }
+      _elevator->initStart();
+      while (!act_blue_nodes.empty()) {
+        for (typename std::vector<BlueNode>::iterator b_it =
+               act_blue_nodes.begin();
+             b_it != act_blue_nodes.end(); ++b_it) {
+          _elevator->initAddItem(*b_it);
+          for (OutArcIt oi(_bpg, *b_it); oi!=INVALID; ++oi) {
+            RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(oi));
+            BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r]));
+            if (b != *b_it && !reached[b]) {
+              next_blue_nodes.push_back(b);
+              reached.set(b, true);
+            }
+          }
+        }
+        _elevator->initNewLevel();
+        next_blue_nodes.swap(act_blue_nodes);
+        next_blue_nodes.clear();
+      }
+      _elevator->initFinish();
+      // Activate blue nodes with excess, i. e. those
+      // which get flow from more than one red node.
+      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
+        if (_cas[b].size() > 1 && (*_elevator)[b] < _elevator->maxLevel()) {
+          _elevator->activate(b);
+        }
+        _ca[b] = _cas[b].begin();
+      }
+    }
+    PreflowBpMatching() {}
+  public:
+    /// \brief Constructor
+    ///
+    /// Constructs the class on the given bipartite graph.
+    PreflowBpMatching(const BpGraph& bpg) : _bpg(bpg),
+                                            _matching(0),
+                                            _local_matching(false),
+                                            _elevator(0),
+                                            _local_elevator(false),
+                                            _cas(bpg),
+                                            _ca(bpg),
+                                            _flow(bpg, INVALID)
+    {}
+    /// \brief Destructor
+    ///
+    /// Destructor.
+    ~PreflowBpMatching() {
+      destroyStructures();
+    }
+    /// \brief Sets the matching map
+    ///
+    /// Sets the matching map. It should contain a valid matching,
+    /// for example the empty matching is fine.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated map,
+    /// of course.
+    ///
+    /// \return <tt>(*this)</tt>
+    PreflowBpMatching& matchingMap(MatchingMap& map) {
+      _matching = &map;
+      _local_matching = false;
+      return *this;
+    }
+    /// \brief Returns a const reference to the matching map.
+    ///
+    /// Returns a const reference to the matching map, which contains
+    /// the matching edge for each node (or INVALID for unmatched nodes).
+    const MatchingMap& matchingMap() const {
+      return *_matching;
+    }
+    /// \brief Sets the elevator used by algorithm.
+    ///
+    /// Sets the elevator.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated elevator,
+    /// of course.
+    /// \return <tt>(*this)</tt>
+    PreflowBpMatching& elevator(Elevator_t& elevator) {
+      _elevator = &elevator;
+      _local_elevator = false;
+      return *this;
+    }
+    /// \brief Returns a const reference to the elevator.
+    ///
+    /// Returns a const reference to the elevator, which keeps track
+    /// of the distance labels of nodes.
+    const Elevator_t& elevator() const {
+      return *_elevator;
+    }
+    /// \brief Initializes the algorithm
+    ///
+    /// Allocates the elevator and the matching map if necessary,
+    /// and initializes the elevator.
+    ///
+    /// \pre \ref init() is not called.
+    void init() {
+      createStructures();
+      initElevator();
+    }
+    /// \brief Smarter initialization of the matching.
+    ///
+    /// Allocates the matching map if necessary, and initializes
+    /// the algorithm with a trivially found matching: iterating
+    /// on every edge, take it in if both endnodes are unmatched.
+    ///
+    /// \return Size of the so found matching.
+    ///
+    /// \note heuristicInit() replaces init() regarding the preconditions.
+    int heuristicInit() {
+      createStructures();
+      int size = 0;
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        bool matched = false;
+        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
+          if ((*_matching)[_bpg.target(oi)] == INVALID) {
+            _matching->set(_bpg.source(oi), oi);
+            _matching->set(_bpg.target(oi), oi);
+            matched = true;
+            ++size;
+          }
+        }
+      }
+      initElevator();
+      return size;
+    }
+    /// \brief Initialize the matching from a map.
+    ///
+    /// Allocates the matching map if necessary, and
+    /// initializes the algorithm with a matching given as a bool edgemap,
+    /// in which an edge is true if it is in the matching.
+    /// Also initializes the elevator.
+    ///
+    /// \return \c true if the matching is valid.
+    ///
+    /// \note matchingInit() replaces init() regarding the preconditions.
+    bool matchingInit(const BoolEdgeMap& matching) {
+      createStructures();
+      bool valid = true;
+      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
+        if (matching[it]) {
+          Node red = _bpg.redNode(it);
+          Node blue = _bpg.blueNode(it);
+          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
+            valid = false;
+          } else {
+            _matching->set(red, it);
+            _matching->set(blue, it);
+          }
+        }
+      }
+      initElevator();
+      return valid;
+    }
+    /// \brief Perform a push/relabel operation on a node
+    ///
+    /// Perform a push/relabel operation on a node.
+    ///
+    /// \pre \ref init() is called, and the node is active.
+    void push_relabel(const BlueNode& chosen) {
+      Arc new_flow = INVALID;
+      while (_ca[chosen] != _cas[chosen].end() && new_flow == INVALID) {
+        OutArcIt r_oi(_bpg, _bpg.source(*_ca[chosen]));
+        while (r_oi != INVALID && new_flow == INVALID) {
+          if ((*_elevator)[_bpg.asBlueNodeUnsafe(_bpg.target(r_oi))] ==
+              (*_elevator)[chosen] - 1) {
+            new_flow = r_oi;
+          }
+          ++r_oi;
+        }
+        if (new_flow == INVALID) ++_ca[chosen];
+      }
+      // If an incident blue node is found on one level lower,
+      // then perform a push, i. e. flip the flow from the
+      // current arc to 'new_flow'.
+      if (new_flow != INVALID) {
+        _ca[chosen] = _cas[chosen].erase(_ca[chosen]);
+        if (_cas[chosen].size() == 1) {
+          _elevator->deactivate(chosen);
+        }
+        _flow[_bpg.asRedNodeUnsafe(_bpg.source(new_flow))] = new_flow;
+        BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(new_flow));
+        _cas[b].push_front(new_flow);
+        if (_cas[b].size() == 2) {
+          _elevator->activate(b);
+        }
+      }
+      // If none found, then relabel.
+      else {
+        int min = _elevator->maxLevel() + 1;
+        for (CurrentArcIterator it1 = _cas[chosen].begin();
+             it1 != _cas[chosen].end();
+             ++it1) {
+          for (OutArcIt it2(_bpg, _bpg.asRedNodeUnsafe(_bpg.source(*it1)));
+               it2 != INVALID;
+               ++it2) {
+            BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(it2));
+            if (b != chosen && min > (*_elevator)[b]) {
+              min = (*_elevator)[b];
+            }
+          }
+        }
+        int chosen_level = (*_elevator)[chosen];
+        if (min >= _elevator->maxLevel() - 1) {
+          _elevator->lift(chosen, _elevator->maxLevel());
+          _elevator->deactivate(chosen);
+        } else {
+          _elevator->lift(chosen, min+1);
+        }
+        if (_elevator->emptyLevel(chosen_level)) {
+          _elevator->liftToTop(chosen_level);
+        }
+        // Restore the current arc iterator.
+        _ca[chosen] = _cas[chosen].begin();
+      }
+    }
+    /// \brief Executes the algorithm
+    ///
+    /// Executes push/relabel operation on the highest active node
+    /// until there is. After start(), writeMatching() should be called to
+    /// remove unnecessary flow and choose a valid maximum matching.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    void start() {
+      while (_elevator->highestActive() != INVALID)
+        push_relabel(_elevator->highestActive());
+    }
+    /// \brief Deduces and writes a matching from the preflow
+    ///
+    /// The underlying preflow algorithm does not use the matching map,
+    /// so before using the matching query functions ( \ref mate(), \ref
+    /// matchingSize(), etc.), you should explicitly call this function
+    /// to write a matching based on the current preflow to the matching
+    /// map. When the preflow results an ambiguous matching (a node has
+    /// incoming flow on more than one arc), it chooses arbitrary one.
+    void writeMatching() {
+      for (NodeIt it(_bpg); it != INVALID; ++it) {
+        _matching->set(it, INVALID);
+      }
+      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
+        if (!_cas[b].empty()) {
+          Arc a = _cas[b].front();
+          _matching->set(_bpg.source(a), a);
+          _matching->set(_bpg.target(a), a);
+        }
+      }
+    }
+    /// \brief Runs the algorithm.
+    ///
+    /// bp1.run() is just a shorthand for:
+    ///
+    ///\code
+    /// bp1.heuristicInit();
+    /// bp1.start();
+    /// bp1.writeMatching();
+    ///\endcode
+    void run() {
+      heuristicInit();
+      start();
+      writeMatching();
+    }
+    /// \brief Size of the matching
+    ///
+    /// Returns the size of the current matching.
+    /// Function runs in \f$O(n)\f$ time.
+    int matchingSize() const {
+      int size = 0;
+      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
+        if ((*_matching)[it] != INVALID) ++size;
+      }
+      return size;
+    }
+    /// \brief Return \c true if the given edge is in the matching.
+    ///
+    /// This function returns \c true if the given edge is in the current
+    /// matching.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpg.redNode(edge)];
+    }
+    /// \brief Return the matching edge incident to the given node.
+    ///
+    /// This function returns the matching edge incident to the
+    /// given node in the current matching or \c INVALID if the node is
+    /// not covered by the matching.
+    Edge matching(const Node& n) const {
+      return (*_matching)[n];
+    }
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the current
+    /// matching or \c INVALID if the node is not covered by the matching.
+    Node mate(const Node& n) const {
+      return (*_matching)[n] != INVALID ?
+        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
+    }
+    /// \brief Query a vertex cover.
+    ///
+    /// This function finds a vertex cover based on the current matching,
+    /// and in \c cover sets node in the cover \c true and \c false
+    /// otherwise. In case of maximal matching, this will be a minimal
+    /// vertex cover.
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return The number of nodes in the cover.
+    int cover(BoolNodeMap &cover) const {
+      int count = 0;
+      std::queue<RedNode> q;
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
+        cover[b] = false;
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          cover[r] = false;
+          q.push(r);
+        } else {
+          cover[r] = true;
+          ++count;
+        }
+      }
+      while (!q.empty()) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (cover[blue]) continue;
+          cover[blue] = true;
+          ++count;
+          if ((*_matching)[blue] != INVALID) {
+            RedNode nr = _bpg.redNode((*_matching)[blue]);
+            q.push(nr);
+            cover[nr] = false;
+            --count;
+          }
+        }
+        q.pop();
+      }
+      return count;
+    }
+    /// \brief Query a barrier of red nodes.
+    ///
+    /// This function tries to create a barrier of red nodes, which should:
+    ///  - contain every unmatched red nodes,
+    ///  - contain exactly that many matched red nodes as much blue nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool redBarrier(BoolRedNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<RedNode> q;
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          barrier[r] = true;
+          q.push(r);
+        } else {
+          barrier[r] = false;
+        }
+      }
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (exist &= ((*_matching)[blue] != INVALID)) {
+            RedNode r = _bpg.redNode((*_matching)[blue]);
+            if (!barrier[r]) {
+              q.push(r);
+              barrier[r] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+      return exist;
+    }
+    /// \brief Query a barrier of blue nodes.
+    ///
+    /// This function tries to create a barrier of blue nodes, which should:
+    ///  - contain every unmatched blue nodes,
+    ///  - contain exactly that many matched blue nodes as much red nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool blueBarrier(BoolBlueNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<BlueNode> q;
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        if ((*_matching)[b] == INVALID) {
+          barrier[b] = true;
+          q.push(b);
+        } else {
+          barrier[b] = false;
+        }
+      }
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
+          if (exist &= ((*_matching)[r] != INVALID)) {
+            BlueNode b = _bpg.blueNode((*_matching)[r]);
+            if (!barrier[b]) {
+              q.push(b);
+              barrier[b] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+      return exist;
+    }
+  };
+}
+#endif

test/CMakeLists.txt

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt

-                      a
   path_test
   planarity_test
   preflow_test
+  preflow_bp_matching_test
   radix_sort_test
   random_test
   suurballe_test

new file test/preflow_bp_matching_test.cc

diff --git a/test/preflow_bp_matching_test.cc b/test/preflow_bp_matching_test.cc
new file mode 100644

-                      -
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2012
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+#include <iostream>
+#include <sstream>
+#include <list>
+#include <lemon/random.h>
+#include <lemon/preflow_bp_matching.h>
+#include <lemon/smart_graph.h>
+#include <lemon/concepts/bpgraph.h>
+#include <lemon/concepts/maps.h>
+#include <lemon/lgf_reader.h>
+#include "test_tools.h"
+using namespace std;
+using namespace lemon;
+const int lgfn = 3;
+const string lgf[lgfn] = {
+  // Empty graph
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+  // No red nodes
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "1\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+  // No blue nodes
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+};
+const string u_lgf =
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "@blue_nodes\n"
+  "label\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "@edges\n"
+  "label\n"
+  "1 4 1\n"
+  "1 5 2\n"
+  "2 5 3\n"
+  "2 6 4\n"
+  "2 7 5\n"
+  "3 7 6\n"
+;
+void checkPreflowBpMatchingCompile() {
+  typedef concepts::BpGraph BpGraph;
+  BPGRAPH_TYPEDEFS(BpGraph);
+  typedef PreflowBpMatching<BpGraph> PBM;
+  BpGraph graph;
+  RedNode red;
+  BlueNode blue;
+  Node node;
+  Edge edge;
+  BoolEdgeMap init_matching(graph);
+  PBM pbm_test(graph);
+  const PBM& const_pbm_test = pbm_test;
+  PBM::MatchingMap matching_map(graph, INVALID);
+  PBM::Elevator_t elevator(graph);
+  BpGraph::NodeMap<bool> cov(graph);
+  int cover_size;
+  BoolRedNodeMap rb(graph);
+  BoolBlueNodeMap bb(graph);
+  pbm_test.matchingMap(matching_map);
+  pbm_test.elevator(elevator);
+  pbm_test.init();
+  pbm_test.heuristicInit();
+  pbm_test.matchingInit(init_matching);
+  pbm_test.start();
+  pbm_test.push_relabel(blue);
+  pbm_test.run();
+  pbm_test.writeMatching();
+  const_pbm_test.matchingSize();
+  const_pbm_test.matching(node);
+  const_pbm_test.matching(edge);
+  const_pbm_test.mate(red);
+  const_pbm_test.mate(blue);
+  const_pbm_test.mate(node);
+  const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap();
+  edge = max_matching[node];
+  const PBM::Elevator_t& elev = const_pbm_test.elevator();
+  blue = elev.highestActive();
+  cover_size = const_pbm_test.cover(cov);
+  const_pbm_test.redBarrier(rb);
+  const_pbm_test.blueBarrier(bb);
+}
+typedef SmartBpGraph BpGraph;
+typedef PreflowBpMatching<BpGraph> PBM;
+BPGRAPH_TYPEDEFS(BpGraph);
+void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) {
+  BoolBlueNodeMap matched(bpg, false);
+  for (RedNodeIt it(bpg); it != INVALID; ++it) {
+    if (pbm.matching(it) != INVALID) {
+      check(pbm.mate(pbm.mate(it)) == it, "Wrong matching");
+      matched.set(bpg.asBlueNode(pbm.mate(it)), true);
+    }
+  }
+  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
+    // != is used as logical xor here
+    check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching");
+  }
+}
+void checkCover(const BpGraph& bpg, const PBM& pbm) {
+  int cov_size;
+  BoolNodeMap cov(bpg);
+  cov_size = pbm.cover(cov);
+  BoolEdgeMap covered(bpg, false);
+  for (NodeIt n(bpg); n!=INVALID; ++n) {
+    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
+      covered[e] = true;
+    }
+  }
+  bool all = true;
+  for (EdgeIt e(bpg); e!=INVALID; ++e) {
+    all&=covered[e];
+  }
+  check(all, "Invalid cover");
+  check(pbm.matchingSize() == cov_size, "Suboptimal matching.");
+}
+void checkBarriers(const BpGraph& bpg, const PBM& pbm) {
+  BoolRedNodeMap rb(bpg);
+  BoolBlueNodeMap bb(bpg);
+  check(pbm.redBarrier(rb), "Suboptimal matching.");
+  check(pbm.blueBarrier(bb), "Suboptimal matching.");
+  BoolNodeMap reached(bpg, false);
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (!rb[r]) continue;
+    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (!bb[b]) continue;
+    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  int ur = 0, ub = 0, rc = 0, bc = 0;
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (rb[r]) ++rc;
+    if (pbm.mate(r) == INVALID) ++ur;
+    if (reached[r]) --bc;
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (bb[b]) ++bc;
+    if (pbm.mate(b) == INVALID) ++ub;
+    if (reached[b]) --rc;
+  }
+  check(rc == ur, "Red barrier is wrong.");
+  check(bc == ub, "Blue barrier is wrong.");
+}
+void checkMatching(const BpGraph& bpg, const PBM& pbm) {
+  checkMatchingValidity(bpg, pbm);
+  checkBarriers(bpg, pbm);
+  checkCover(bpg, pbm);
+}
+int main() {
+  // Checking hard wired extremal graphs
+  for (int i=0; i<lgfn; ++i) {
+    BpGraph bpg;
+    istringstream lgfs(lgf[i]);
+    bpGraphReader(bpg, lgfs).run();
+    PBM pbm(bpg);
+    pbm.run();
+    checkMatching(bpg, pbm);
+  }
+  // Checking some use cases
+  BpGraph bpg;
+  istringstream u_lgfs(u_lgf);
+  bpGraphReader(bpg, u_lgfs).run();
+  {
+    PBM pbm(bpg);
+    pbm.init();
+    BpGraph::BlueNode b = pbm.elevator().highestActive();
+    if (b != INVALID) {
+      pbm.push_relabel(b);
+    }
+    pbm.start();
+    pbm.writeMatching();
+    check(pbm.elevator().highestActive() == INVALID, "Error in elevator");
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    pbm.heuristicInit();
+    BpGraph::BlueNode b = pbm.elevator().highestActive();
+    if (b != INVALID) {
+      pbm.push_relabel(b);
+    }
+    pbm.start();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    PBM::MatchingMap m(bpg, INVALID);
+    pbm.matchingMap(m);
+    pbm.run();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    BpGraph::EdgeMap<bool> init_matching(bpg, false);
+    BpGraph::Edge e;
+    bpg.first(e);
+    init_matching[e] = true;
+    check(pbm.matchingInit(init_matching),
+          "Wrong result from matchingInit()");
+    pbm.start();
+    pbm.writeMatching();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    PBM::MatchingMap m(bpg, INVALID);
+    pbm.matchingMap(m);
+    BpGraph::EdgeMap<bool> init_matching(bpg, true);
+    check(!pbm.matchingInit(init_matching),
+          "Wrong result from matchingInit()");
+    pbm.start();
+    pbm.writeMatching();
+    checkMatching(bpg, pbm);
+  }
+  // Checking some random generated graphs
+  const int random_test_n = 10;
+  const int max_red_n = 1000;
+  const int min_red_n = 10;
+  const int max_blue_n = 1000;
+  const int min_blue_n = 10;
+  for (int i=0; i<random_test_n; ++i) {
+    int red_n = rnd.integer(min_red_n, max_red_n);
+    int blue_n = rnd.integer(min_blue_n, max_blue_n);
+    double prob = rnd();
+    SmartBpGraph bpg;
+    for (int j=0; j<red_n; ++j) {
+      bpg.addRedNode();
+    }
+    for (int j=0; j<blue_n; ++j) {
+      bpg.addBlueNode();
+    }
+    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
+      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (rnd() < prob/10.) {
+          bpg.addEdge(r,b);
+        }
+      }
+    }
+    PBM pbm(bpg);
+    pbm.run();
+    checkMatching(bpg, pbm);
+    // do not always use the highest active node
+    PBM pbm2(bpg);
+    PBM::Elevator_t elev(bpg);
+    pbm2.elevator(elev);
+    pbm2.init();
+    while (elev.highestActive() != INVALID) {
+      for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (elev.active(b)) pbm2.push_relabel(b);
+      }
+    }
+    pbm2.writeMatching();
+    checkMatching(bpg, pbm2);
+  }
+  return 0;
+}

new file lemon/preflow_bp_matching_2.h

# HG changeset patch
# User Daniel Poroszkai <poroszd@inf.elte.hu>
# Date 1330129984 -3600
# Node ID ba89bdcf59f238ea9917f1d5d9515760c4847df7
# Parent  b348df30553a63d6ecda56c1c2f1e1f18401f87a
Other variant of preflow based bipartite matching

diff --git a/lemon/preflow_bp_matching_2.h b/lemon/preflow_bp_matching_2.h
new file mode 100644

-                      -
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2012
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+#ifndef LEMON_PREFLOW_BP_MATCHING_2_H
+#define LEMON_PREFLOW_BP_MATCHING_2_H
+#include <lemon/core.h>
+#include <lemon/elevator.h>
+#include <vector>
+#include <queue>
+/// \ingroup matching
+/// \file
+/// \brief A preflow based bipartite matching algorithm.
+namespace lemon {
+  /// \brief A preflow based bipartite matching algorithm
+  ///
+  /// Finds maximum cardinality matching in a given bipartite
+  /// graph using specialized preflow algorithm.
+  template <typename BPG>
+  class PreflowBpMatching2 {
+  public:
+    /// Type of the bipartite graph
+    typedef BPG BpGraph;
+    /// Type of the matching map
+    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
+    /// Type of the elevator.
+    typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t;
+  private:
+    TEMPLATE_BPGRAPH_TYPEDEFS(BPG);
+  protected:
+    typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow;
+    typedef typename BPG::template BlueNodeMap<int> Excess;
+    typedef std::vector<typename BPG::BlueNode> BlueNodes;
+    typedef std::vector<typename BPG::Arc> BlueFlow;
+    typedef typename BPG::template BlueNodeMap<BlueFlow> IncomingFlows;
+    const BpGraph& _bpg;
+    MatchingMap* _matching;
+    bool _local_matching;
+    Elevator_t* _elevator;
+    bool _local_elevator;
+    Preflow _flow;
+    Excess _excess;
+    void createStructures() {
+      if (!_matching) {
+        _matching = new MatchingMap(_bpg, INVALID);
+        _local_matching = true;
+      }
+      if (!_elevator) {
+        _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1,
+                                                  countBlueNodes(_bpg)));
+        _local_elevator = true;
+      }
+    }
+    void destroyStructures() {
+      if (_local_matching) {
+        delete _matching;
+      }
+      if (_local_elevator) {
+        delete _elevator;
+      }
+    }
+    /// \brief Initialize the elevator
+    ///
+    /// When this function is called, the elevator is initialized with
+    /// the exact distance labels corresponding the current matching.
+    void initElevator() {
+      BlueNodes act_blue_nodes, next_blue_nodes;
+      BoolBlueNodeMap reached(_bpg, false);
+      IncomingFlows incoming(_bpg);
+      // Set the initial flow:
+      // choose arbitrary outgoing arc for every unmatched red node.
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        Arc a = (*_matching)[r] != INVALID ?
+          _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r);
+        if (a != INVALID) {
+          _flow[r] = a;
+          BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(a));
+          incoming[b].push_back(a);
+          if (++_excess[b] == 2) {
+            act_blue_nodes.push_back(b);
+            reached.set(b, true);
+          }
+        }
+      }
+      // Now we start a Bfs from the blue nodes with excess
+      // (those which get flow from more than one red node) in the
+      // residual graph, to supply the exact distance labels to the elevator.
+      _elevator->initStart();
+      while (!act_blue_nodes.empty()) {
+        for (typename BlueNodes::iterator b_it = act_blue_nodes.begin();
+             b_it != act_blue_nodes.end(); ++b_it) {
+          _elevator->initAddItem(*b_it);
+          for (typename BlueFlow::iterator b_in = incoming[*b_it].begin();
+               b_in != incoming[*b_it].end(); ++b_in) {
+            RedNode r = _bpg.asRedNodeUnsafe(_bpg.source(*b_in));
+            for (OutArcIt r_oi(_bpg, r); r_oi!=INVALID; ++r_oi) {
+              BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(r_oi));
+              if (!reached[b]) {
+                next_blue_nodes.push_back(b);
+                reached.set(b, true);
+              }
+            }
+          }
+        }
+        _elevator->initNewLevel();
+        next_blue_nodes.swap(act_blue_nodes);
+        next_blue_nodes.clear();
+      }
+      _elevator->initFinish();
+      // Activate blue nodes with 0 excess, i. e. those
+      // which does not get any flow.
+      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
+        if (_excess[b] == 0  && (*_elevator)[b] < _elevator->maxLevel()) {
+          _elevator->activate(b);
+        }
+      }
+    }
+    PreflowBpMatching2() {}
+  public:
+    /// \brief Constructor
+    ///
+    /// Constructs the class on the given bipartite graph.
+    PreflowBpMatching2(const BpGraph& bpg) : _bpg(bpg),
+                                             _matching(0),
+                                             _local_matching(false),
+                                             _elevator(0),
+                                             _local_elevator(false),
+                                             _flow(bpg, INVALID),
+                                             _excess(_bpg, 0)
+    {}
+    /// \brief Destructor
+    ///
+    /// Destructor.
+    ~PreflowBpMatching2() {
+      destroyStructures();
+    }
+    /// \brief Sets the matching map
+    ///
+    /// Sets the matching map. It should contain a valid matching,
+    /// for example the empty matching is fine.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated map,
+    /// of course.
+    ///
+    /// \return <tt>(*this)</tt>
+    PreflowBpMatching2& matchingMap(MatchingMap& map) {
+      _matching = &map;
+      _local_matching = false;
+      return *this;
+    }
+    /// \brief Returns a const reference to the matching map.
+    ///
+    /// Returns a const reference to the matching map, which contains
+    /// the matching edge for each node (or INVALID for unmatched nodes).
+    const MatchingMap& matchingMap() const {
+      return *_matching;
+    }
+    /// \brief Sets the elevator used by algorithm.
+    ///
+    /// Sets the elevator.
+    /// If you don't use this function before calling \ref init(),
+    /// an instance will be allocated automatically.
+    /// The destructor deallocates this automatically allocated elevator,
+    /// of course.
+    /// \return <tt>(*this)</tt>
+    PreflowBpMatching2& elevator(Elevator_t& elevator) {
+      _elevator = &elevator;
+      _local_elevator = false;
+      return *this;
+    }
+    /// \brief Returns a const reference to the elevator.
+    ///
+    /// Returns a const reference to the elevator, which keeps track
+    /// of the distance labels of nodes.
+    const Elevator_t& elevator() const {
+      return *_elevator;
+    }
+    /// \brief Initializes the algorithm
+    ///
+    /// Allocates the elevator and the matching map if necessary,
+    /// and initializes the elevator.
+    ///
+    /// \pre \ref init() is not called.
+    void init() {
+      createStructures();
+      initElevator();
+    }
+    /// \brief Smarter initialization of the matching.
+    ///
+    /// Allocates the matching map if necessary, and initializes
+    /// the algorithm with a trivially found matching: iterating
+    /// on every edge, take it in if both endnodes are unmatched.
+    ///
+    /// \return Size of the so found matching.
+    ///
+    /// \note heuristicInit() replaces init() regarding the preconditions.
+    int heuristicInit() {
+      createStructures();
+      int size = 0;
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        bool matched = false;
+        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
+          if ((*_matching)[_bpg.target(oi)] == INVALID) {
+            _matching->set(_bpg.source(oi), oi);
+            _matching->set(_bpg.target(oi), oi);
+            matched = true;
+            ++size;
+          }
+        }
+      }
+      initElevator();
+      return size;
+    }
+    /// \brief Initialize the matching from a map.
+    ///
+    /// Allocates the matching map if necessary, and
+    /// initializes the algorithm with a matching given as a bool edgemap,
+    /// in which an edge is true if it is in the matching.
+    /// Also initializes the elevator.
+    ///
+    /// \return \c true if the matching is valid.
+    ///
+    /// \note matchingInit() replaces init() regarding the preconditions.
+    bool matchingInit(const BoolEdgeMap& matching) {
+      createStructures();
+      bool valid = true;
+      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
+        if (matching[it]) {
+          Node red = _bpg.redNode(it);
+          Node blue = _bpg.blueNode(it);
+          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
+            valid = false;
+          } else {
+            _matching->set(red, it);
+            _matching->set(blue, it);
+          }
+        }
+      }
+      initElevator();
+      return valid;
+    }
+    /// \brief Perform a pull/relabel operation on a node
+    ///
+    /// Perform a pull/relabel operation on a node:
+    /// search an incident blue node on one level lower,
+    /// and flip its flow to the 'chosen', or relabel
+    /// if none found.
+    ///
+    /// \pre \ref init() is called, and the node is active.
+    void pull_relabel(const BlueNode& chosen) {
+      int min_inc_level = _elevator->maxLevel() + 1;
+      for (OutArcIt b_oi(_bpg, chosen); b_oi!=INVALID; ++b_oi) {
+        RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(b_oi));
+        BlueNode lower = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r]));
+        if ((*_elevator)[(lower)] < (*_elevator)[chosen]) {
+          if (--_excess[lower] == 0) {
+            _elevator->activate(lower);
+          }
+          _flow[r] = _bpg.direct(b_oi, true);
+          ++_excess[chosen];
+          _elevator->deactivate(chosen);
+          break;
+        } else {
+          if ((*_elevator)[lower] < min_inc_level &&
+              lower != chosen) {
+            min_inc_level = (*_elevator)[lower];
+          }
+        }
+      }
+      // If none found, then relabel.
+      if (_excess[chosen] == 0) {
+        int chosen_level = (*_elevator)[chosen];
+        if (min_inc_level >= _elevator->maxLevel() - 1) {
+          _elevator->lift(chosen, _elevator->maxLevel());
+          _elevator->deactivate(chosen);
+        } else {
+          _elevator->lift(chosen, min_inc_level+1);
+        }
+        if (_elevator->emptyLevel(chosen_level)) {
+          _elevator->liftToTop(chosen_level);
+        }
+      }
+    }
+    /// \brief Executes the algorithm
+    ///
+    /// Executes pull/relabel operation on the highest active node
+    /// until there is. After start(), writeMatching() should be called to
+    /// remove unnecessary flow and choose a valid maximum matching.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    void start() {
+      while (_elevator->highestActive() != INVALID)
+        pull_relabel(_elevator->highestActive());
+    }
+    /// \brief Deduces and writes a matching from the preflow
+    ///
+    /// The underlying preflow algorithm does not use the matching map,
+    /// so before using the matching query functions ( \ref mate(), \ref
+    /// matchingSize(), etc.), you should explicitly call this function
+    /// to write a matching based on the current preflow to the matching
+    /// map. When the preflow results an ambiguous matching (a node has
+    /// incoming flow on more than one arc), it chooses arbitrary one.
+    void writeMatching() {
+      for (NodeIt it(_bpg); it != INVALID; ++it) {
+        _matching->set(it, INVALID);
+      }
+      for (RedNodeIt r(_bpg); r != INVALID; ++r) {
+        if (_flow[r] != INVALID) {
+          BlueNode mate = _bpg.blueNode(_flow[r]);
+          if (_excess[mate] > 1) {
+            --_excess[mate];
+          } else {
+            _matching->set(mate, _flow[r]);
+            _matching->set(r, _flow[r]);
+          }
+        }
+      }
+    }
+    /// \brief Runs the algorithm.
+    ///
+    /// bp1.run() is just a shorthand for:
+    ///
+    ///\code
+    /// bp1.heuristicInit();
+    /// bp1.start();
+    /// bp1.writeMatching();
+    ///\endcode
+    void run() {
+      heuristicInit();
+      start();
+      writeMatching();
+    }
+    /// \brief Size of the matching
+    ///
+    /// Returns the size of the current matching.
+    /// Function runs in \f$O(n)\f$ time.
+    int matchingSize() const {
+      int size = 0;
+      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
+        if ((*_matching)[it] != INVALID) ++size;
+      }
+      return size;
+    }
+    /// \brief Return \c true if the given edge is in the matching.
+    ///
+    /// This function returns \c true if the given edge is in the current
+    /// matching.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_bpg.redNode(edge)];
+    }
+    /// \brief Return the matching edge incident to the given node.
+    ///
+    /// This function returns the matching edge incident to the
+    /// given node in the current matching or \c INVALID if the node is
+    /// not covered by the matching.
+    Edge matching(const Node& n) const {
+      return (*_matching)[n];
+    }
+    /// \brief Return the mate of the given node.
+    ///
+    /// This function returns the mate of the given node in the current
+    /// matching or \c INVALID if the node is not covered by the matching.
+    Node mate(const Node& n) const {
+      return (*_matching)[n] != INVALID ?
+        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
+    }
+    /// \brief Query a vertex cover.
+    ///
+    /// This function finds a vertex cover based on the current matching,
+    /// and in \c cover sets node in the cover \c true and \c false
+    /// otherwise. In case of maximal matching, this will be a minimal
+    /// vertex cover.
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return The number of nodes in the cover.
+    int cover(BoolNodeMap &cover) const {
+      int count = 0;
+      std::queue<RedNode> q;
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
+        cover[b] = false;
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          cover[r] = false;
+          q.push(r);
+        } else {
+          cover[r] = true;
+          ++count;
+        }
+      }
+      while (!q.empty()) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (cover[blue]) continue;
+          cover[blue] = true;
+          ++count;
+          if ((*_matching)[blue] != INVALID) {
+            RedNode nr = _bpg.redNode((*_matching)[blue]);
+            q.push(nr);
+            cover[nr] = false;
+            --count;
+          }
+        }
+        q.pop();
+      }
+      return count;
+    }
+    /// \brief Query a barrier of red nodes.
+    ///
+    /// This function tries to create a barrier of red nodes, which should:
+    ///  - contain every unmatched red nodes,
+    ///  - contain exactly that many matched red nodes as much blue nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool redBarrier(BoolRedNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<RedNode> q;
+      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
+        if ((*_matching)[r] == INVALID) {
+          barrier[r] = true;
+          q.push(r);
+        } else {
+          barrier[r] = false;
+        }
+      }
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
+          if (exist &= ((*_matching)[blue] != INVALID)) {
+            RedNode r = _bpg.redNode((*_matching)[blue]);
+            if (!barrier[r]) {
+              q.push(r);
+              barrier[r] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+      return exist;
+    }
+    /// \brief Query a barrier of blue nodes.
+    ///
+    /// This function tries to create a barrier of blue nodes, which should:
+    ///  - contain every unmatched blue nodes,
+    ///  - contain exactly that many matched blue nodes as much red nodes
+    ///    there are in its neighbor set.
+    ///
+    /// Barrier exist if and only if the matching is maximal. If exist, the
+    /// map is set true for nodes of the barrier, and false for the others.
+    ///
+    /// Function runs in \f$O(n + e)\f$ time.
+    ///
+    /// \return \c true if a barrier found.
+    bool blueBarrier(BoolBlueNodeMap &barrier) const {
+      bool exist = true;
+      std::queue<BlueNode> q;
+      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
+        if ((*_matching)[b] == INVALID) {
+          barrier[b] = true;
+          q.push(b);
+        } else {
+          barrier[b] = false;
+        }
+      }
+      while (!q.empty() && exist) {
+        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
+          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
+          if (exist &= ((*_matching)[r] != INVALID)) {
+            BlueNode b = _bpg.blueNode((*_matching)[r]);
+            if (!barrier[b]) {
+              q.push(b);
+              barrier[b] = true;
+            }
+          }
+        }
+        q.pop();
+      }
+      return exist;
+    }
+  };
+}
+#endif

test/CMakeLists.txt

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt

-                      a
   path_test
   planarity_test
   preflow_test
+  preflow_bp_matching_2_test
   preflow_bp_matching_test
   radix_sort_test
   random_test

new file test/preflow_bp_matching_2_test.cc

diff --git a/test/preflow_bp_matching_2_test.cc b/test/preflow_bp_matching_2_test.cc
new file mode 100644

-                      -
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2012
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+#include <iostream>
+#include <sstream>
+#include <list>
+#include <lemon/random.h>
+#include <lemon/preflow_bp_matching_2.h>
+#include <lemon/smart_graph.h>
+#include <lemon/concepts/bpgraph.h>
+#include <lemon/concepts/maps.h>
+#include <lemon/lgf_reader.h>
+#include "test_tools.h"
+using namespace std;
+using namespace lemon;
+const int lgfn = 3;
+const string lgf[lgfn] = {
+  // Empty graph
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+  // No red nodes
+  "@red_nodes\n"
+  "label\n"
+  "\n"
+  "@blue_nodes\n"
+  "label\n"
+  "1\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+  // No blue nodes
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "@blue_nodes\n"
+  "label\n"
+  "\n"
+  "@edges\n"
+  "  label\n"
+  "\n",
+};
+const string u_lgf =
+  "@red_nodes\n"
+  "label\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "@blue_nodes\n"
+  "label\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "@edges\n"
+  "label\n"
+  "1 4 1\n"
+  "1 5 2\n"
+  "2 5 3\n"
+  "2 6 4\n"
+  "2 7 5\n"
+  "3 7 6\n"
+;
+void checkPreflowBpMatching2Compile() {
+  typedef concepts::BpGraph BpGraph;
+  BPGRAPH_TYPEDEFS(BpGraph);
+  typedef PreflowBpMatching2<BpGraph> PBM;
+  BpGraph graph;
+  RedNode red;
+  BlueNode blue;
+  Node node;
+  Edge edge;
+  BoolEdgeMap init_matching(graph);
+  PBM pbm_test(graph);
+  const PBM& const_pbm_test = pbm_test;
+  PBM::MatchingMap matching_map(graph, INVALID);
+  PBM::Elevator_t elevator(graph);
+  BpGraph::NodeMap<bool> cov(graph);
+  int cover_size;
+  BoolRedNodeMap rb(graph);
+  BoolBlueNodeMap bb(graph);
+  pbm_test.matchingMap(matching_map);
+  pbm_test.elevator(elevator);
+  pbm_test.init();
+  pbm_test.heuristicInit();
+  pbm_test.matchingInit(init_matching);
+  pbm_test.start();
+  pbm_test.pull_relabel(blue);
+  pbm_test.run();
+  pbm_test.writeMatching();
+  const_pbm_test.matchingSize();
+  const_pbm_test.matching(node);
+  const_pbm_test.matching(edge);
+  const_pbm_test.mate(red);
+  const_pbm_test.mate(blue);
+  const_pbm_test.mate(node);
+  const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap();
+  edge = max_matching[node];
+  const PBM::Elevator_t& elev = const_pbm_test.elevator();
+  blue = elev.highestActive();
+  cover_size = const_pbm_test.cover(cov);
+  const_pbm_test.redBarrier(rb);
+  const_pbm_test.blueBarrier(bb);
+}
+typedef SmartBpGraph BpGraph;
+typedef PreflowBpMatching2<BpGraph> PBM;
+BPGRAPH_TYPEDEFS(BpGraph);
+void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) {
+  BoolBlueNodeMap matched(bpg, false);
+  for (RedNodeIt it(bpg); it != INVALID; ++it) {
+    if (pbm.matching(it) != INVALID) {
+      check(pbm.mate(pbm.mate(it)) == it, "Wrong matching");
+      matched.set(bpg.asBlueNode(pbm.mate(it)), true);
+    }
+  }
+  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
+    // != is used as logical xor here
+    check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching");
+  }
+}
+void checkCover(const BpGraph& bpg, const PBM& pbm) {
+  int cov_size;
+  BoolNodeMap cov(bpg);
+  cov_size = pbm.cover(cov);
+  BoolEdgeMap covered(bpg, false);
+  for (NodeIt n(bpg); n!=INVALID; ++n) {
+    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
+      covered[e] = true;
+    }
+  }
+  bool all = true;
+  for (EdgeIt e(bpg); e!=INVALID; ++e) {
+    all&=covered[e];
+  }
+  check(all, "Invalid cover");
+  check(pbm.matchingSize() == cov_size, "Suboptimal matching.");
+}
+void checkBarriers(const BpGraph& bpg, const PBM& pbm) {
+  BoolRedNodeMap rb(bpg);
+  BoolBlueNodeMap bb(bpg);
+  check(pbm.redBarrier(rb), "Suboptimal matching.");
+  check(pbm.blueBarrier(bb), "Suboptimal matching.");
+  BoolNodeMap reached(bpg, false);
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (!rb[r]) continue;
+    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (!bb[b]) continue;
+    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
+      reached.set(bpg.target(a), true);
+    }
+  }
+  int ur = 0, ub = 0, rc = 0, bc = 0;
+  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
+    if (rb[r]) ++rc;
+    if (pbm.mate(r) == INVALID) ++ur;
+    if (reached[r]) --bc;
+  }
+  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+    if (bb[b]) ++bc;
+    if (pbm.mate(b) == INVALID) ++ub;
+    if (reached[b]) --rc;
+  }
+  check(rc == ur, "Red barrier is wrong.");
+  check(bc == ub, "Blue barrier is wrong.");
+}
+void checkMatching(const BpGraph& bpg, const PBM& pbm) {
+  checkMatchingValidity(bpg, pbm);
+  checkBarriers(bpg, pbm);
+  checkCover(bpg, pbm);
+}
+int main() {
+  // Checking hard wired extremal graphs
+  for (int i=0; i<lgfn; ++i) {
+    BpGraph bpg;
+    istringstream lgfs(lgf[i]);
+    bpGraphReader(bpg, lgfs).run();
+    PBM pbm(bpg);
+    pbm.run();
+    checkMatching(bpg, pbm);
+  }
+  // Checking some use cases
+  BpGraph bpg;
+  istringstream u_lgfs(u_lgf);
+  bpGraphReader(bpg, u_lgfs).run();
+  {
+    PBM pbm(bpg);
+    pbm.init();
+    BpGraph::BlueNode b = pbm.elevator().highestActive();
+    if (b != INVALID) {
+      pbm.pull_relabel(b);
+    }
+    pbm.start();
+    pbm.writeMatching();
+    check(pbm.elevator().highestActive() == INVALID, "Error in elevator");
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    pbm.heuristicInit();
+    BpGraph::BlueNode b = pbm.elevator().highestActive();
+    if (b != INVALID) {
+      pbm.pull_relabel(b);
+    }
+    pbm.start();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    PBM::MatchingMap m(bpg, INVALID);
+    pbm.matchingMap(m);
+    pbm.run();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    BpGraph::EdgeMap<bool> init_matching(bpg, false);
+    BpGraph::Edge e;
+    bpg.first(e);
+    init_matching[e] = true;
+    check(pbm.matchingInit(init_matching),
+          "Wrong result from matchingInit()");
+    pbm.start();
+    pbm.writeMatching();
+    checkMatching(bpg, pbm);
+  }
+  {
+    PBM pbm(bpg);
+    PBM::MatchingMap m(bpg, INVALID);
+    pbm.matchingMap(m);
+    BpGraph::EdgeMap<bool> init_matching(bpg, true);
+    check(!pbm.matchingInit(init_matching),
+          "Wrong result from matchingInit()");
+    pbm.start();
+    pbm.writeMatching();
+    checkMatching(bpg, pbm);
+  }
+  // Checking some random generated graphs
+  const int random_test_n = 10;
+  const int max_red_n = 1000;
+  const int min_red_n = 10;
+  const int max_blue_n = 1000;
+  const int min_blue_n = 10;
+  for (int i=0; i<random_test_n; ++i) {
+    int red_n = rnd.integer(min_red_n, max_red_n);
+    int blue_n = rnd.integer(min_blue_n, max_blue_n);
+    double prob = rnd();
+    SmartBpGraph bpg;
+    for (int j=0; j<red_n; ++j) {
+      bpg.addRedNode();
+    }
+    for (int j=0; j<blue_n; ++j) {
+      bpg.addBlueNode();
+    }
+    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
+      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (rnd() < prob/10.) {
+          bpg.addEdge(r,b);
+        }
+      }
+    }
+    PBM pbm(bpg);
+    pbm.run();
+    checkMatching(bpg, pbm);
+    // do not always use the highest active node
+    PBM pbm2(bpg);
+    PBM::Elevator_t elev(bpg);
+    pbm2.elevator(elev);
+    pbm2.init();
+    while (elev.highestActive() != INVALID) {
+      for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
+        if (elev.active(b)) pbm2.pull_relabel(b);
+      }
+    }
+    pbm2.writeMatching();
+    checkMatching(bpg, pbm2);
+  }
+  return 0;
+}

tools/bp-matching-benchmark.cc

# HG changeset patch
# User Daniel Poroszkai <poroszd@inf.elte.hu>
# Date 1329867694 -3600
# Node ID 93fa3bb078dd0fb7e198d73c42999ceccd1bbe27
# Parent  add2bdd7f7a635e1f7d080dd433d342d9271e779
Benchmark tool extended

diff --git a/tools/bp-matching-benchmark.cc b/tools/bp-matching-benchmark.cc

-                      a
 #include <lemon/hopcroft_karp.h>
 #include <lemon/matching.h>
 #include <lemon/preflow.h>
+#include <lemon/preflow_bp_matching.h>
+#include <lemon/preflow_bp_matching_2.h>
 #include <lemon/bp_matching.h>
 #include <lemon/cost_scaling.h>
 #include <lemon/capacity_scaling.h>
 …
 using namespace std;
 using namespace lemon;
+void check(bool c, const string &msg) {
+  if (!c) {
+    std::cerr << "Error: " << msg << "\n";
+  }
+}
 template <typename BpGraph, typename Digraph>
 void create_network(const BpGraph &bpg,
 …
 ///
 /// Depending on the input, it runs the \ref lemon::HopcroftKarp
 /// "Hopcroft-Karp", the \ref lemon::MaxMatching "general matching"
+/// and the \ref lemon::Preflow "Preflow" algorithms for finding
+/// and three algorithms based on preflow (\ref lemon::Preflow "preflow
+/// on extended graph", \ref lemon::PreflowBpMatching,
+/// \ref lemon::PreflowBpMatching2") for finding
 /// a maximum cardinality bipartite matching, or the
 /// \ref lemon::MaxWeightedBpMatching "max. weighted matching
 /// for sparse bipartite graphs", the \ref lemon::MaxWeightedMatching
 …
 /// simplex" algorithms to find the maximum weighted matching.
 int main() {
   DimacsDescriptor desc = dimacsType(cin);
-  int red_n, blue_n, edge_n = desc.edgeNum;
   if (desc.type == DimacsDescriptor::UBM) {
     // unweighted bipartite matching
+    Timer hk_t(false), mm_t(false), pf_t(false);
+    int hk_r, mm_r, pf_r;
+    Timer hk_t(false),
+      mm_t(false),
+      pf_t(false),
+      pbm_t(false),
+      pbm2_t(false);
     BpGraph *bpg = new BpGraph();
     readDimacsUbm(cin, *bpg, desc);
-    red_n = countRedNodes(*bpg);
-    blue_n = countBlueNodes(*bpg);
+    {
       hk_t.start();
       HopcroftKarp<BpGraph> hk(*bpg);
       hk.run();
-      hk_r = hk.matchingSize();
       hk_t.stop();
+      BpGraph::NodeMap<bool> m(*bpg);
+      check(hk.matchingSize() == hk.cover(m),
+            "Suboptimal matching in Hopcroft-Karp.");
+    }
+    {
       mm_t.start();
       MaxMatching<BpGraph> mm(*bpg);
       mm.run();
-      mm_r = mm.matchingSize();
       mm_t.stop();
+    }
+    {
+      pbm_t.start();
+      PreflowBpMatching<BpGraph> pbm(*bpg);
+      pbm.run();
+      pbm_t.stop();
+      BpGraph::NodeMap<bool> m(*bpg);
+      check(pbm.matchingSize() == pbm.cover(m),
+            "Suboptimal matching in push-relabel.");
+    }
+    {
+      pbm2_t.start();
+      PreflowBpMatching2<BpGraph> pbm2(*bpg);
+      pbm2.run();
+      pbm2_t.stop();
+      BpGraph::NodeMap<bool> m(*bpg);
+      check(pbm2.matchingSize() == pbm2.cover(m),
+            "Suboptimal matching in pull-relabel.");
+    }
     Digraph *network = new Digraph();
     Digraph::Node source, target;
     create_network(*bpg, *network, source, target);
-    delete bpg;
+    {
       pf_t.start();
       Preflow<Digraph, ConstMap<Digraph::Arc, long> >
         pf(*network, constMap<Digraph::Arc, long>(1), source, target);
       pf.run();
-      pf_r = pf.flowValue();
       pf_t.stop();
+    }
     delete network;
     cout << "Benchmarking in unweighted case, using a bipartite graph\n"
+         << "with " << red_n << " red, " << blue_n << " blue nodes, and "
+         << edge_n << " edges.\n"
+         << "with " << countRedNodes(*bpg) << " red, "
+         << countBlueNodes(*bpg) << " blue nodes, and "
+         << desc.edgeNum << " edges.\n"
          << "--------------------------------------------------------\n"
+         << "Algorithm used       Matching size       Time\n"
+         << "Hopcroft-Karp            " << hk_r
+         << "         " << hk_t.realTime() << "\n"
+         << "General matching         " << mm_r
+         << "         " << mm_t.realTime() << "\n"
+         << "Preflow                  " << pf_r
+         << "         " << pf_t.realTime() << endl;
+         << "Algorithm used               Time\n"
+         << "Hopcroft-Karp            " << hk_t.realTime() << "\n"
+         << "General matching         " << mm_t.realTime() << "\n"
+         << "Push-relabel matching    " << pbm_t.realTime() << "\n"
+         << "Pull-relabel matching    " << pbm2_t.realTime() << "\n"
+         << "Preflow                  " << pf_t.realTime() << endl;
+    delete bpg;
   } else if (desc.type == DimacsDescriptor::WBM) {
     // weighted bipartite matching
     Timer bm_t(false),
 …
     BpGraph *bpg = new BpGraph();
     BpGraph::EdgeMap<long> *weight = new BpGraph::EdgeMap<long>(*bpg);
     readDimacsWbm(cin, *bpg, *weight, desc);
-    red_n = countRedNodes(*bpg);
-    blue_n = countBlueNodes(*bpg);
+    {
       bm_t.start();
 …
     delete supply;
     cout << "Benchmarking on weighted bipartite matching problem on a "
+         << "graph\nwith " << red_n << " red, " << blue_n
+         << " blue nodes and " << edge_n << " edges\n"
+         << "graph\nwith " << countRedNodes(*bpg) << " red, "
+         << countBlueNodes(*bpg) << " blue nodes and "
+         << desc.edgeNum << " edges\n"
          << "--------------------------------------------------------\n"
          << "Algorithm used       Maximum weight       Time\n"
          << "Bipartite matching       " << bm_r

Download in other formats:

Original Format