Ticket #168: hopcroft_karp-ec5a4ceb4481.patch

lemon/Makefile.am

@@ -97 +97 @@
         lemon/karp_mmc.h \
         lemon/kruskal.h \
         lemon/hao_orlin.h \
+        lemon/hopcroft_karp.h \
         lemon/lgf_reader.h \
         lemon/lgf_writer.h \
         lemon/list_graph.h \

new file lemon/hopcroft_karp.h

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2011
+ * 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_HOPCROFT_KARP_H
+#define LEMON_HOPCROFT_KARP_H
+
+#include <lemon/core.h>
+#include <list>
+
+/// \ingroup matching
+/// \file
+/// \brief Hopcroft-Karp algorithm.
+namespace lemon {
+
+  /// \brief The Hopcroft-Karp bipartite matching algorithm
+  ///
+  /// Finds maximal matching in a given bipartite
+  /// graph using the Hopcroft-Karp algorithm,
+  /// having \f$O(e\sqrt{n})\f$ complexity.
+  template <typename BPG>
+  class HopcroftKarp {
+  public:
+      /// Type of the bipartite graph
+    typedef BPG BpGraph;
+    /// Type of the matching map
+    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
+  private:
+    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
+
+    const BpGraph& _bpg;
+    MatchingMap* _matching;
+    bool _local_matching;
+
+  protected:
+    void createStructures() {
+      if (!_matching) {
+        _matching = new MatchingMap(_bpg, INVALID);
+        _local_matching = true;
+      }
+    }
+
+    void destroyStructures() {
+      if (_local_matching) {
+        delete _matching;
+      }
+    }
+
+    HopcroftKarp() {}
+
+  public:
+    /// \brief Constructor
+    ///
+    /// Constructs the class on the given bipartite graph.
+    HopcroftKarp(const BpGraph& bpg) : _bpg(bpg),
+                                       _matching(0),
+                                       _local_matching(false)
+    {}
+
+
+    /// \brief Destructor
+    ///
+    /// Destructor.
+    ~HopcroftKarp() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the matching map
+    ///
+    /// Sets the matching map.
+    /// 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) {
+      _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 every node (and \c INVALID for the
+    /// unmatched nodes).
+    const MatchingMap& matchingMap() const {
+      return *_matching;
+    }
+
+
+    /// \brief Initializes the algorithm
+    ///
+    /// Allocates the matching map if necessary, and sets
+    /// to empty matching.
+    ///
+    /// \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.
+    ///
+    /// 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 found matching.
+    ///
+    /// \note heuristicInit() replaces init() regarding the preconditions.
+    int heuristicInit() {
+      init();
+      int size = 0;
+
+      for (BlueIt 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;
+            }
+          }
+        }
+      }
+      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.
+    ///
+    /// 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();
+      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);
+          }
+        }
+      }
+      return valid;
+    }
+
+    /// \brief Executes an augmenting phase
+    ///
+    /// Searches a maximal set of vertex disjoint shortest alternating paths.
+    /// Meaning:
+    ///  - alternating: connents an unmatched red- and an unmatched blue node,
+    ///    and exactly every second edge is in the current matching;
+    ///  - shortest: contain a minimal number of edges;
+    ///  - vertex disjoint: every vertex belong to at most one path;
+    ///  - maximal set: a set of path is maximal when it is not expandable
+    ///    further (and not when there does not exist a set with
+    ///    more vertex disjoint shortest alternating paths).
+    ///
+    /// After a maximal set is found, it applies the augmenting paths,
+    /// so edges of the matching are taken out, the others are put in
+    /// the matching.
+    ///
+    /// \return The length of the found augmenting paths, or -1 when none
+    /// found, which occurs if and only if the current matching is maximal.
+    ///
+    /// \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;
+
+      for (RedIt r_it(_bpg); r_it != INVALID; ++r_it) {
+        if ((*_matching)[r_it] == INVALID) {
+          act_rednodes.push_front(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
+      // augmenting paths.
+      std::list<BlueNode> path_ends;
+
+      // Layer counter
+      int cur_level = 0;
+
+      // 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);
+            }
+          }
+        }
+
+        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;
+      }
+
+      // The current_edge structure assures that edges iterated at most once
+      typename BpGraph::template BlueMap<IncEdgeIt*> current_edge(_bpg, 0);
+
+      // Stack for the Dfs
+      std::list<BlueNode> stack;
+
+      // 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();
+
+      while (n_it != path_ends.end()) {
+        stack.push_front(*n_it);
+        path_found = false;
+        while (!stack.empty() && !path_found) {
+          BlueNode b = stack.front();
+          if (current_edge[b] == 0) {
+            current_edge[b] = new IncEdgeIt(_bpg, b);
+          } 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]));
+            }
+          }
+        }
+        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);
+          }
+        }
+
+        stack.clear();
+        ++n_it;
+      }
+      return cur_level - 1;
+    }
+
+    /// \brief Executes the algorithm
+    ///
+    /// It runs augmenting phases until the optimal solution is reached.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    void start() {
+      while (augment() > 0) {}
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// hk.run() is just a shorthand for:
+    ///
+    ///\code
+    /// hk.heuristicInit();
+    /// hk.start();
+    ///\endcode
+    void run() {
+      heuristicInit();
+      start();
+    }
+
+    /// \brief Size of the matching
+    ///
+    /// Returns the size of the current matching.
+    int matchingSize() const {
+      int size = 0;
+      for (RedIt 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;
+    }
+  };
+
+}
+#endif

test/CMakeLists.txt

@@ -33 +33 @@
   graph_utils_test
   hao_orlin_test
   heap_test
+  hopcroft_karp_test
   kruskal_test
   maps_test
   matching_test

test/Makefile.am

@@ -31 +31 @@
         test/graph_utils_test \
         test/hao_orlin_test \
         test/heap_test \
+        test/hopcroft_karp_test \
         test/kruskal_test \
         test/maps_test \
         test/matching_test \
@@ -85 +86 @@
 test_heap_test_SOURCES = test/heap_test.cc
 test_kruskal_test_SOURCES = test/kruskal_test.cc
 test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
+test_hopcroft_karp_test_SOURCES = test/hopcroft_karp_test.cc
 test_lp_test_SOURCES = test/lp_test.cc
 test_maps_test_SOURCES = test/maps_test.cc
 test_mip_test_SOURCES = test/mip_test.cc

new file test/hopcroft_karp_test.cc

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2011
+ * 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/hopcroft_karp.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 checkHopcroftKarpCompile() {
+  typedef concepts::BpGraph BpGraph;
+  BPGRAPH_TYPEDEFS(BpGraph);
+  typedef HopcroftKarp<BpGraph> HK;
+
+  BpGraph graph;
+  RedNode red;
+  BlueNode blue;
+  Node node;
+  Edge edge;
+  BoolEdgeMap init_matching(graph);
+  HK hk_test(graph);
+  const HK& const_hk_test = hk_test;
+  HK::MatchingMap matching_map(graph);
+
+  hk_test.matchingMap(matching_map);
+  hk_test.init();
+  hk_test.heuristicInit();
+  hk_test.matchingInit(init_matching);
+
+  hk_test.start();
+  hk_test.augment();
+  hk_test.run();
+
+  const_hk_test.matchingSize();
+  const_hk_test.matching(node);
+  const_hk_test.matching(edge);
+  const_hk_test.mate(red);
+  const_hk_test.mate(blue);
+  const_hk_test.mate(node);
+  const HK::MatchingMap& max_matching = const_hk_test.matchingMap();
+  edge = max_matching[node];
+}
+
+typedef SmartBpGraph BpGraph;
+typedef HopcroftKarp<BpGraph> HK;
+BPGRAPH_TYPEDEFS(BpGraph);
+
+void checkMatchingValidity(const BpGraph& bpg, const HK& hk) {
+  BoolBlueMap matched(bpg, false);
+  for (RedIt it(bpg); it != INVALID; ++it) {
+    if (hk.matching(it) != INVALID) {
+      check(hk.mate(hk.mate(it)) == it, "Wrong matching");
+      matched.set(hk.mate(it), true);
+    }
+  }
+  for (BlueIt it(bpg); it != INVALID; ++it) {
+    // != is used as logical xor here
+    check(matched[it] != (hk.matching(it) == INVALID), "Wrong matching");
+  }
+
+}
+
+void checkMatchingOptimality(const BpGraph& bpg, const HK& hk) {
+  list<RedNode> reds;
+  list<BlueNode> blues;
+  BoolBlueMap reached(bpg, false);
+
+  for (RedIt it(bpg); it != INVALID; ++it) {
+    if (hk.matching(it) == INVALID) {
+      reds.push_front(it);
+    }
+  }
+
+  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(hk.mate(blue));
+    }
+  }
+
+}
+
+int main() {
+  // Checking hard wired extremal graphs
+  for (int i=0; i<lgfn; ++i) {
+    BpGraph bpg;
+    istringstream lgfs(lgf[i]);
+    bpGraphReader(bpg, lgfs).run();
+    HK hk(bpg);
+    hk.run();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+  }
+
+  // Checking some use cases
+  BpGraph bpg;
+  istringstream u_lgfs(u_lgf);
+  bpGraphReader(bpg, u_lgfs).run();
+  {
+    HK hk(bpg);
+    hk.init();
+    hk.augment();
+    hk.start();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+  }
+  {
+    HK hk(bpg);
+    hk.heuristicInit();
+    hk.augment();
+    hk.start();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+  }
+  {
+    HK hk(bpg);
+    HK::MatchingMap m(bpg);
+    hk.matchingMap(m);
+    hk.run();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+  }
+  {
+    HK hk(bpg);
+    BpGraph::EdgeMap<bool> init_matching(bpg, false);
+    BpGraph::Edge e;
+    bpg.first(e);
+    init_matching[e] = true;
+    check(hk.matchingInit(init_matching), "Wrong result from matchingInit()");
+    hk.start();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+  }
+  {
+    HK hk(bpg);
+    HK::MatchingMap m(bpg);
+    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);
+  }
+
+  // Checking some random generated graphs
+  const int random_test_n = 10;
+  const int max_red_n = 100;
+  const int min_red_n = 10;
+  const int max_blue_n = 100;
+  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 i=0; i<red_n; ++i) {
+      bpg.addRedNode();
+    }
+    for (int i=0; i<blue_n; ++i) {
+      bpg.addBlueNode();
+    }
+    for (SmartBpGraph::RedIt r(bpg); r!=INVALID; ++r) {
+      for (SmartBpGraph::BlueIt b(bpg); b!=INVALID; ++b) {
+        if (rnd() < prob) {
+          bpg.addEdge(r,b);
+        }
+      }
+    }
+
+    HK hk(bpg);
+    hk.run();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+  }
+
+
+  return 0;
+}
+
+