Ticket #168: h_k_45d2e43f52ea.patch

lemon/Makefile.am

@@ -92 +92 @@
         lemon/grid_graph.h \
         lemon/grosso_locatelli_pullan_mc.h \
         lemon/hartmann_orlin_mmc.h \
+        lemon/hopcroft_karp.h \
         lemon/howard_mmc.h \
         lemon/hypercube_graph.h \
         lemon/karp_mmc.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:
+    /// The bipartite graph type of the algorithm.
+    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 run(),
+    /// 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.
+    void init() {
+      createStructures();
+      for (NodeIt it(_bpg); it!=INVALID; ++it) {
+        _matching->set(it, INVALID);
+      }
+    }
+
+    /// \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.
+    /// \return \c true if the matching is valid.
+    bool matchingInit(const BoolEdgeMap& matching) {
+      createStructures();
+
+      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
+        if (matching[it]) {
+          Node red = _bpg.redNode(it);
+          if ((*_matching)[red] != INVALID) return false;
+          _matching->set(red, it);
+
+          Node blue = _bpg.blueNode(it);
+          if ((*_matching)[blue] != INVALID) return false;
+          _matching->set(blue, it);
+        }
+      }
+      return true;
+    }
+
+    /// \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 %True when augmenting paths are found, and %False when none,
+    /// which occurs if and only if the current matching is maximal.
+    ///
+    /// \pre \ref init() or \ref matchingInit() must be called before using
+    /// this function.
+    bool augment() {
+      BoolNodeMap free_node(_bpg, false);
+      BoolBlueMap reached_node(_bpg, false);
+      IntRedMap level(_bpg, -1);
+      std::list<RedNode> act_rednodes;
+      std::list<BlueNode> act_bluenodes;
+
+      for (NodeIt it(_bpg); it!=INVALID; ++it) {
+        if ((*_matching)[it] == INVALID) {
+          free_node.set(it, true);
+        }
+      }
+      for (RedIt it(_bpg); it!=INVALID; ++it) {
+        if (free_node[it]) {
+          act_rednodes.push_front(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<Node> path_ends;
+
+      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 (!reached_node[blue]) {
+              act_bluenodes.push_front(blue);
+              reached_node[blue] = true;
+              path_found |= free_node[blue];
+            }
+          }
+        }
+
+        if (!path_found) {
+          if (act_bluenodes.empty()) return false;
+          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 (free_node[*it]) path_ends.push_front(*it);
+          }
+          // Now path_ends contains nodes that are possible ends of
+          // shortest alternating paths.
+        }
+        ++cur_level;
+      }
+      // List of possible augmenting paths (paths are list of edges).
+      std::list<std::list<Edge> > paths;
+      paths.resize(path_ends.size());
+
+      // Searching backward, starting from the previously found
+      // blue nodes, we build vertex disjoint alternating paths.
+      // Paths that cannot be built further are erased from the list.
+      do {
+        typename std::list<std::list<Edge> >::iterator p_it = paths.begin();
+        typename std::list<Node>::iterator n_it = path_ends.begin();
+        --cur_level;
+        while (p_it != paths.end()) {
+          IncEdgeIt it(_bpg, *n_it);
+          while (it!=INVALID && level[_bpg.redNode(it)] != cur_level) ++it;
+          if (it != INVALID) {
+            Node red = _bpg.redNode(it);
+            level[red] = -1; // Don't visit this node again.
+            *n_it = red;
+            p_it->push_front(it);
+            ++n_it;
+            ++p_it;
+          } else {
+            p_it = paths.erase(p_it);
+            n_it = path_ends.erase(n_it);
+          }
+        }
+
+        if (cur_level > 0) {
+          p_it = paths.begin();
+          n_it = path_ends.begin();
+          while (p_it != paths.end()) {
+            Node blue = _bpg.blueNode((*_matching)[*n_it]);
+            *n_it = blue;
+            p_it->push_front((*_matching)[*n_it]);
+            ++n_it;
+            ++p_it;
+          }
+        }
+      } while (cur_level > 0);
+
+      // Now 'paths' contains a maximal set of shortest augmenting paths.
+      for (typename std::list<std::list<Edge> >::iterator p_it = paths.begin();
+           p_it != paths.end();
+           ++p_it) {
+        for (typename std::list<Edge>::iterator e_it = p_it->begin();
+             e_it != p_it->end();) {
+          _matching->set(_bpg.redNode(*e_it), *e_it);
+          _matching->set(_bpg.blueNode(*e_it), *e_it);
+          ++e_it;
+          if (e_it != p_it->end()) ++e_it;
+        }
+      }
+      return true;
+    }
+
+    /// \brief Executes the algorithm
+    ///
+    /// It runs augmenting phases until the optimal solution is reached.
+    ///
+    /// \pre \ref init() or \ref matchingInit() must be called before using
+    /// this function.
+    void start() {
+      while (augment()) {}
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// hk.run() is just a shorthand for:
+    ///
+    ///\code
+    /// hk.init();
+    /// hk.start();
+    ///\endcode
+    void run() {
+      init();
+      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

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

test/Makefile.am

@@ -30 +30 @@
         test/graph_utils_test \
         test/hao_orlin_test \
         test/heap_test \
+        test/hopcroft_karp_test \
         test/kruskal_test \
         test/maps_test \
         test/matching_test \
@@ -83 +84 @@
 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",
+};
+
+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.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 generateBpGraph(BpGraph& bpg, int rn, int bn, double p) {
+  for (int i=0; i<rn; ++i) bpg.addRedNode();
+  for (int i=0; i<bn; ++i) bpg.addBlueNode();
+  for (RedIt r_it(bpg); r_it != INVALID; ++r_it) {
+    for (BlueIt b_it(bpg); b_it != INVALID; ++b_it) {
+      if (rnd.boolean(p)) bpg.addEdge(r_it, b_it);
+    }
+  }
+}
+
+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 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 random generated graphs
+  const int random_test_n = 5;
+  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();
+    BpGraph bpg;
+    generateBpGraph(bpg, red_n, blue_n, prob);
+    HK hk(bpg);
+    hk.run();
+    checkMatchingValidity(bpg, hk);
+    checkMatchingOptimality(bpg, hk);
+  }
+  return 0;
+}
+
+