COIN-OR::LEMON - Graph Library

Ticket #168: h_k_45d2e43f52ea.patch

File h_k_45d2e43f52ea.patch, 17.3 KB (added by Poroszkai Daniel, 10 years ago)

In the previous patch the algorithm was erroneous, it is corrected.

  • lemon/Makefile.am

    # HG changeset patch
    # User Poroszkai Daniel <poroszd@gmail.com>
    # Date 1301950286 -7200
    # Node ID 45d2e43f52ea2bbcc2f41146fac78d59b5ab7503
    # Parent  9312d6c89d02e8a5dc3c8c97332a4b347ca00e9c
    A working Hopcroft-Karp bipartite matching.
    
    diff --git a/lemon/Makefile.am b/lemon/Makefile.am
    a b  
    9292        lemon/grid_graph.h \
    9393        lemon/grosso_locatelli_pullan_mc.h \
    9494        lemon/hartmann_orlin_mmc.h \
     95        lemon/hopcroft_karp.h \
    9596        lemon/howard_mmc.h \
    9697        lemon/hypercube_graph.h \
    9798        lemon/karp_mmc.h \
  • new file lemon/hopcroft_karp.h

    diff --git a/lemon/hopcroft_karp.h b/lemon/hopcroft_karp.h
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2011
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_HOPCROFT_KARP_H
     20#define LEMON_HOPCROFT_KARP_H
     21
     22#include <lemon/core.h>
     23#include <list>
     24
     25/// \ingroup matching
     26/// \file
     27/// \brief Hopcroft-Karp algorithm.
     28
     29namespace lemon {
     30
     31  /// \brief The Hopcroft-Karp bipartite matching algorithm
     32  ///
     33  /// Finds maximal matching in a given bipartite
     34  /// graph using the Hopcroft-Karp algorithm,
     35  /// having \f$O(e\sqrt{n})\f$ complexity.
     36  template <typename BPG>
     37  class HopcroftKarp {
     38  public:
     39    /// The bipartite graph type of the algorithm.
     40    typedef BPG BpGraph;
     41    /// Type of the matching map.
     42    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
     43
     44  private:
     45    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     46
     47    const BpGraph& _bpg;
     48    MatchingMap* _matching;
     49    bool _local_matching;
     50
     51  protected:
     52
     53    void createStructures() {
     54      if (!_matching) {
     55        _matching = new MatchingMap(_bpg, INVALID);
     56        _local_matching = true;
     57      }
     58    }
     59
     60    void destroyStructures() {
     61      if (_local_matching) {
     62        delete _matching;
     63      }
     64    }
     65
     66    HopcroftKarp() {}
     67
     68  public:
     69
     70    /// \brief Constructor
     71    ///
     72    /// Constructs the class on the given bipartite graph.
     73    HopcroftKarp(const BpGraph& bpg) : _bpg(bpg),
     74                                       _matching(0),
     75                                       _local_matching(false)
     76    {}
     77
     78    /// \brief Destructor
     79    ///
     80    /// Destructor.
     81    ~HopcroftKarp() {
     82      destroyStructures();
     83    }
     84
     85    /// \brief Sets the matching map
     86    ///
     87    /// Sets the matching map.
     88    /// If you don't use this function before calling \ref run(),
     89    /// an instance will be allocated automatically.
     90    /// The destructor deallocates this automatically allocated map,
     91    /// of course.
     92    /// This member is not to initialize the algorithm with a valid
     93    /// matching; use \ref matchingInit() instead.
     94    /// \return <tt>(*this)</tt>
     95    HopcroftKarp& matchingMap(MatchingMap& map) {
     96      _matching = &map;
     97      _local_matching = false;
     98      return *this;
     99    }
     100
     101    /// \brief Returns a const reference to the matching map
     102    ///
     103    /// Returns a const reference to the matching map, which contains
     104    /// the matching edge for every node (and \c INVALID for the
     105    /// unmatched nodes).
     106    const MatchingMap& matchingMap() const {
     107      return *_matching;
     108    }
     109
     110    /// \brief Initializes the algorithm
     111    ///
     112    /// Allocates the matching map if necessary, and sets
     113    /// to empty matching.
     114    void init() {
     115      createStructures();
     116      for (NodeIt it(_bpg); it!=INVALID; ++it) {
     117        _matching->set(it, INVALID);
     118      }
     119    }
     120
     121    /// \brief Initialize the matching from a map.
     122    ///
     123    /// Allocates the matching map if necessary, and
     124    /// initializes the algorithm with a matching given as a bool edgemap,
     125    /// in which an edge is true if it is in the matching.
     126    /// \return \c true if the matching is valid.
     127    bool matchingInit(const BoolEdgeMap& matching) {
     128      createStructures();
     129
     130      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
     131        if (matching[it]) {
     132          Node red = _bpg.redNode(it);
     133          if ((*_matching)[red] != INVALID) return false;
     134          _matching->set(red, it);
     135
     136          Node blue = _bpg.blueNode(it);
     137          if ((*_matching)[blue] != INVALID) return false;
     138          _matching->set(blue, it);
     139        }
     140      }
     141      return true;
     142    }
     143
     144    /// \brief Executes an augmenting phase
     145    ///
     146    /// Searches a maximal set of vertex disjoint shortest alternating paths.
     147    /// Meaning:
     148    ///  - alternating: connents an unmatched red- and an unmatched blue node,
     149    ///    and exactly every second edge is in the current matching;
     150    ///  - shortest: contain a minimal number of edges;
     151    ///  - vertex disjoint: every vertex belong to at most one path;
     152    ///  - maximal set: a set of path is maximal when it is not expandable
     153    ///    further (and not when there does not exist a set with
     154    ///    more vertex disjoint shortest alternating paths).
     155    ///
     156    /// After a maximal set is found, it applies the augmenting paths,
     157    /// so edges of the matching are taken out, the others are put in
     158    /// the matching.
     159    ///
     160    /// \return %True when augmenting paths are found, and %False when none,
     161    /// which occurs if and only if the current matching is maximal.
     162    ///
     163    /// \pre \ref init() or \ref matchingInit() must be called before using
     164    /// this function.
     165    bool augment() {
     166      BoolNodeMap free_node(_bpg, false);
     167      BoolBlueMap reached_node(_bpg, false);
     168      IntRedMap level(_bpg, -1);
     169      std::list<RedNode> act_rednodes;
     170      std::list<BlueNode> act_bluenodes;
     171
     172      for (NodeIt it(_bpg); it!=INVALID; ++it) {
     173        if ((*_matching)[it] == INVALID) {
     174          free_node.set(it, true);
     175        }
     176      }
     177      for (RedIt it(_bpg); it!=INVALID; ++it) {
     178        if (free_node[it]) {
     179          act_rednodes.push_front(it);
     180        }
     181      }
     182
     183      // Raise this flag when a shortest augmenting path is found.
     184      bool path_found = false;
     185      // This nodelist will contain the end nodes of the possible
     186      // augmenting paths.
     187      std::list<Node> path_ends;
     188
     189      int cur_level = 0;
     190
     191      // Starting from the unmatched red nodes search for unmatched
     192      // blue nodes, using Bfs (but only on alternating paths).
     193      while (!path_found) {
     194        while (!act_rednodes.empty()) {
     195          RedNode red = act_rednodes.front();
     196          act_rednodes.pop_front();
     197          level[red] = cur_level;
     198          for (IncEdgeIt it(_bpg, red); it!=INVALID; ++it) {
     199            BlueNode blue(_bpg.blueNode(it));
     200            if (!reached_node[blue]) {
     201              act_bluenodes.push_front(blue);
     202              reached_node[blue] = true;
     203              path_found |= free_node[blue];
     204            }
     205          }
     206        }
     207
     208        if (!path_found) {
     209          if (act_bluenodes.empty()) return false;
     210          while (!act_bluenodes.empty()) {
     211            BlueNode blue = act_bluenodes.front();
     212            act_bluenodes.pop_front();
     213            RedNode red = _bpg.redNode((*_matching)[blue]);
     214            act_rednodes.push_front(red);
     215          }
     216        } else {
     217          for (typename std::list<BlueNode>::iterator it=act_bluenodes.begin();
     218               it != act_bluenodes.end();
     219               ++it) {
     220            if (free_node[*it]) path_ends.push_front(*it);
     221          }
     222          // Now path_ends contains nodes that are possible ends of
     223          // shortest alternating paths.
     224        }
     225        ++cur_level;
     226      }
     227      // List of possible augmenting paths (paths are list of edges).
     228      std::list<std::list<Edge> > paths;
     229      paths.resize(path_ends.size());
     230
     231      // Searching backward, starting from the previously found
     232      // blue nodes, we build vertex disjoint alternating paths.
     233      // Paths that cannot be built further are erased from the list.
     234      do {
     235        typename std::list<std::list<Edge> >::iterator p_it = paths.begin();
     236        typename std::list<Node>::iterator n_it = path_ends.begin();
     237        --cur_level;
     238        while (p_it != paths.end()) {
     239          IncEdgeIt it(_bpg, *n_it);
     240          while (it!=INVALID && level[_bpg.redNode(it)] != cur_level) ++it;
     241          if (it != INVALID) {
     242            Node red = _bpg.redNode(it);
     243            level[red] = -1; // Don't visit this node again.
     244            *n_it = red;
     245            p_it->push_front(it);
     246            ++n_it;
     247            ++p_it;
     248          } else {
     249            p_it = paths.erase(p_it);
     250            n_it = path_ends.erase(n_it);
     251          }
     252        }
     253
     254        if (cur_level > 0) {
     255          p_it = paths.begin();
     256          n_it = path_ends.begin();
     257          while (p_it != paths.end()) {
     258            Node blue = _bpg.blueNode((*_matching)[*n_it]);
     259            *n_it = blue;
     260            p_it->push_front((*_matching)[*n_it]);
     261            ++n_it;
     262            ++p_it;
     263          }
     264        }
     265      } while (cur_level > 0);
     266
     267      // Now 'paths' contains a maximal set of shortest augmenting paths.
     268      for (typename std::list<std::list<Edge> >::iterator p_it = paths.begin();
     269           p_it != paths.end();
     270           ++p_it) {
     271        for (typename std::list<Edge>::iterator e_it = p_it->begin();
     272             e_it != p_it->end();) {
     273          _matching->set(_bpg.redNode(*e_it), *e_it);
     274          _matching->set(_bpg.blueNode(*e_it), *e_it);
     275          ++e_it;
     276          if (e_it != p_it->end()) ++e_it;
     277        }
     278      }
     279      return true;
     280    }
     281
     282    /// \brief Executes the algorithm
     283    ///
     284    /// It runs augmenting phases until the optimal solution is reached.
     285    ///
     286    /// \pre \ref init() or \ref matchingInit() must be called before using
     287    /// this function.
     288    void start() {
     289      while (augment()) {}
     290    }
     291
     292    /// \brief Runs the algorithm.
     293    ///
     294    /// hk.run() is just a shorthand for:
     295    ///
     296    ///\code
     297    /// hk.init();
     298    /// hk.start();
     299    ///\endcode
     300    void run() {
     301      init();
     302      start();
     303    }
     304
     305    /// \brief Size of the matching
     306    ///
     307    /// Returns the size of the current matching.
     308    int matchingSize() const {
     309      int size = 0;
     310      for (RedIt it(_bpg); it!=INVALID; ++it) {
     311        if ((*_matching)[it] != INVALID) ++size;
     312      }
     313      return size;
     314    }
     315
     316    /// \brief Return \c true if the given edge is in the matching.
     317    ///
     318    /// This function returns \c true if the given edge is in the current
     319    /// matching.
     320    bool matching(const Edge& edge) const {
     321      return edge == (*_matching)[_bpg.redNode(edge)];
     322    }
     323
     324    /// \brief Return the matching edge incident to the given node.
     325    ///
     326    /// This function returns the matching edge incident to the
     327    /// given node in the current matching or \c INVALID if the node is
     328    /// not covered by the matching.
     329    Edge matching(const Node& n) const {
     330      return (*_matching)[n];
     331    }
     332
     333    /// \brief Return the mate of the given node.
     334    ///
     335    /// This function returns the mate of the given node in the current
     336    /// matching or \c INVALID if the node is not covered by the matching.
     337    Node mate(const Node& n) const {
     338      return (*_matching)[n] != INVALID ?
     339        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
     340    }
     341  };
     342
     343}
     344#endif
  • test/CMakeLists.txt

    diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
    a b  
    2828  graph_utils_test
    2929  hao_orlin_test
    3030  heap_test
     31  hopcroft_karp_test
    3132  kruskal_test
    3233  maps_test
    3334  matching_test
  • test/Makefile.am

    diff --git a/test/Makefile.am b/test/Makefile.am
    a b  
    3030        test/graph_utils_test \
    3131        test/hao_orlin_test \
    3232        test/heap_test \
     33        test/hopcroft_karp_test \
    3334        test/kruskal_test \
    3435        test/maps_test \
    3536        test/matching_test \
     
    8384test_heap_test_SOURCES = test/heap_test.cc
    8485test_kruskal_test_SOURCES = test/kruskal_test.cc
    8586test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
     87test_hopcroft_karp_test_SOURCES = test/hopcroft_karp_test.cc
    8688test_lp_test_SOURCES = test/lp_test.cc
    8789test_maps_test_SOURCES = test/maps_test.cc
    8890test_mip_test_SOURCES = test/mip_test.cc
  • new file test/hopcroft_karp_test.cc

    diff --git a/test/hopcroft_karp_test.cc b/test/hopcroft_karp_test.cc
    new file mode 100644
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2011
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <sstream>
     21#include <list>
     22
     23#include <lemon/random.h>
     24#include <lemon/hopcroft_karp.h>
     25#include <lemon/smart_graph.h>
     26#include <lemon/concepts/bpgraph.h>
     27#include <lemon/concepts/maps.h>
     28#include <lemon/lgf_reader.h>
     29
     30#include "test_tools.h"
     31
     32using namespace std;
     33using namespace lemon;
     34
     35const int lgfn = 3;
     36const string lgf[lgfn] = {
     37  // Empty graph
     38  "@red_nodes\n"
     39  "label\n"
     40  "\n"
     41  "@blue_nodes\n"
     42  "label\n"
     43  "\n"
     44  "@edges\n"
     45  "  label\n"
     46  "\n",
     47
     48  // No red nodes
     49  "@red_nodes\n"
     50  "label\n"
     51  "\n"
     52  "@blue_nodes\n"
     53  "label\n"
     54  "1\n"
     55  "@edges\n"
     56  "  label\n"
     57  "\n",
     58
     59  // No blue nodes
     60  "@red_nodes\n"
     61  "label\n"
     62  "1\n"
     63  "@blue_nodes\n"
     64  "label\n"
     65  "\n"
     66  "@edges\n"
     67  "  label\n"
     68  "\n",
     69};
     70
     71void checkHopcroftKarpCompile() {
     72  typedef concepts::BpGraph BpGraph;
     73  BPGRAPH_TYPEDEFS(BpGraph);
     74  typedef HopcroftKarp<BpGraph> HK;
     75
     76  BpGraph graph;
     77  RedNode red;
     78  BlueNode blue;
     79  Node node;
     80  Edge edge;
     81  BoolEdgeMap init_matching(graph);
     82  HK hk_test(graph);
     83  const HK& const_hk_test = hk_test;
     84  HK::MatchingMap matching_map(graph);
     85
     86  hk_test.matchingMap(matching_map);
     87  hk_test.init();
     88  hk_test.matchingInit(init_matching);
     89  hk_test.start();
     90  hk_test.augment();
     91  hk_test.run();
     92
     93  const_hk_test.matchingSize();
     94  const_hk_test.matching(node);
     95  const_hk_test.matching(edge);
     96  const_hk_test.mate(red);
     97  const_hk_test.mate(blue);
     98  const_hk_test.mate(node);
     99  const HK::MatchingMap& max_matching = const_hk_test.matchingMap();
     100  edge = max_matching[node];
     101}
     102
     103typedef SmartBpGraph BpGraph;
     104typedef HopcroftKarp<BpGraph> HK;
     105BPGRAPH_TYPEDEFS(BpGraph);
     106
     107void generateBpGraph(BpGraph& bpg, int rn, int bn, double p) {
     108  for (int i=0; i<rn; ++i) bpg.addRedNode();
     109  for (int i=0; i<bn; ++i) bpg.addBlueNode();
     110  for (RedIt r_it(bpg); r_it != INVALID; ++r_it) {
     111    for (BlueIt b_it(bpg); b_it != INVALID; ++b_it) {
     112      if (rnd.boolean(p)) bpg.addEdge(r_it, b_it);
     113    }
     114  }
     115}
     116
     117void checkMatchingValidity(const BpGraph& bpg, const HK& hk) {
     118  BoolBlueMap matched(bpg, false);
     119  for (RedIt it(bpg); it != INVALID; ++it) {
     120    if (hk.matching(it) != INVALID) {
     121      check(hk.mate(hk.mate(it)) == it, "Wrong matching");
     122      matched.set(hk.mate(it), true);
     123    }
     124  }
     125  for (BlueIt it(bpg); it != INVALID; ++it) {
     126    // != is used as logical xor here
     127    check(matched[it] != (hk.matching(it) == INVALID), "Wrong matching");
     128  }
     129
     130}
     131
     132void checkMatchingOptimality(const BpGraph& bpg, const HK& hk) {
     133  list<RedNode> reds;
     134  list<BlueNode> blues;
     135  BoolBlueMap reached(bpg, false);
     136
     137  for (RedIt it(bpg); it != INVALID; ++it) {
     138    if (hk.matching(it) == INVALID) {
     139      reds.push_front(it);
     140    }
     141  }
     142
     143  while (!reds.empty()) {
     144    while (!reds.empty()) {
     145      RedNode red = reds.front();
     146      reds.pop_front();
     147      for (IncEdgeIt it(bpg, red); it != INVALID; ++it) {
     148        BlueNode blue = bpg.blueNode(it);
     149        if (!reached[blue]) {
     150          blues.push_front(blue);
     151          reached.set(blue, true);
     152        }
     153      }
     154    }
     155
     156    while (!blues.empty()) {
     157      BlueNode blue = blues.front();
     158      blues.pop_front();
     159      check(hk.matching(blue) != INVALID, "Suboptimal matching");
     160      reds.push_front(hk.mate(blue));
     161    }
     162  }
     163
     164}
     165
     166int main() {
     167  // Checking hard wired graphs
     168  for (int i=0; i<lgfn; ++i) {
     169    BpGraph bpg;
     170    istringstream lgfs(lgf[i]);
     171    bpGraphReader(bpg, lgfs).run();
     172    HK hk(bpg);
     173    hk.run();
     174    checkMatchingValidity(bpg, hk);
     175    checkMatchingOptimality(bpg, hk);
     176  }
     177
     178  // Checking some random generated graphs
     179  const int random_test_n = 5;
     180  const int max_red_n = 100;
     181  const int min_red_n = 10;
     182  const int max_blue_n = 100;
     183  const int min_blue_n = 10;
     184  for (int i=0; i<random_test_n; ++i) {
     185    int red_n = rnd.integer(min_red_n, max_red_n);
     186    int blue_n = rnd.integer(min_blue_n, max_blue_n);
     187    double prob = rnd();
     188    BpGraph bpg;
     189    generateBpGraph(bpg, red_n, blue_n, prob);
     190    HK hk(bpg);
     191    hk.run();
     192    checkMatchingValidity(bpg, hk);
     193    checkMatchingOptimality(bpg, hk);
     194  }
     195  return 0;
     196}
     197
     198