COIN-OR::LEMON - Graph Library

Ticket #168: h_k_1a2daf1f6bf7.patch

File h_k_1a2daf1f6bf7.patch, 17.5 KB (added by Poroszkai Daniel, 10 years ago)

Hopcroft-Karp bipartite matching

  • lemon/Makefile.am

    # HG changeset patch
    # User Poroszkai Daniel <poroszd@gmail.com>
    # Date 1295810211 -3600
    # Node ID 1a2daf1f6bf7ce6b53e7317986fcd91eb814f306
    # Parent  6c08e82b576ec66d1fafdef16971c1103fe719b6
    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 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, true);
     167      BoolNodeMap reached_node(_bpg, false);
     168      std::list<RedNode> act_rednodes;
     169      std::list<BlueNode> act_bluenodes;
     170
     171      for (NodeIt it(_bpg); it!=INVALID; ++it) {
     172        if ((*_matching)[it] != INVALID) {
     173          free_node.set(it, false);
     174        }
     175      }
     176      for (RedIt it(_bpg); it!=INVALID; ++it) {
     177        if (free_node[it]) {
     178          act_rednodes.push_front(it);
     179          reached_node[it] = true;
     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      // Starting from the unmatched red nodes search for unmatched
     190      // blue nodes, using Bfs (but only on alternating paths).
     191      while (!path_found) {
     192        while (!act_rednodes.empty()) {
     193          RedNode red = act_rednodes.front();
     194          act_rednodes.pop_front();
     195          for (IncEdgeIt it(_bpg, red); it!=INVALID; ++it) {
     196            BlueNode blue(_bpg.blueNode(it));
     197            if (!reached_node[blue]) {
     198              act_bluenodes.push_front(blue);
     199              reached_node[blue] = true;
     200              path_found |= free_node[blue];
     201            }
     202          }
     203        }
     204
     205        if (!path_found) {
     206          if (act_bluenodes.empty()) return false;
     207          while (!act_bluenodes.empty()) {
     208            BlueNode blue = act_bluenodes.front();
     209            act_bluenodes.pop_front();
     210            RedNode red = _bpg.redNode((*_matching)[blue]);
     211            reached_node[red] = true;
     212            act_rednodes.push_front(red);
     213          }
     214        } else {
     215          for(typename std::list<BlueNode>::iterator it = act_bluenodes.begin();
     216              it != act_bluenodes.end();
     217              ++it) {
     218            if (free_node[*it]) path_ends.push_front(*it);
     219          }
     220          // Now path_ends contains nodes that are possible ends of
     221          // shortest alternating paths.
     222        }
     223      }
     224      // List of possible augmenting paths (paths are list of edges).
     225      std::list<std::list<Edge> > paths;
     226      paths.resize(path_ends.size());
     227
     228      // Flag will be raised when a shortest augmenting path is found
     229      // (but now in reverse direction)
     230      path_found = false;
     231
     232      // Searching backward, starting from the previously found
     233      // blue nodes, we build vertex disjoint alternating paths.
     234      // Paths that cannot be built further are erased from the list.
     235      while (!path_found) {
     236        typename std::list<std::list<Edge> >::iterator p_it = paths.begin();
     237        typename std::list<Node>::iterator n_it = path_ends.begin();
     238
     239        while (p_it != paths.end()) {
     240          IncEdgeIt it(_bpg, *n_it);
     241          while (it!=INVALID && !reached_node[_bpg.redNode(it)]) ++it;
     242          if (it != INVALID) {
     243            Node red = _bpg.redNode(it);
     244            reached_node[red] = false;
     245            *n_it = red;
     246            p_it->push_front(it);
     247            path_found |= free_node[red];
     248            ++n_it;
     249            ++p_it;
     250          } else {
     251            typename std::list<std::list<Edge> >::iterator p_del(p_it);
     252            typename std::list<Node>::iterator n_del(n_it);
     253            ++p_it;
     254            ++n_it;
     255            paths.erase(p_del);
     256            path_ends.erase(n_del);
     257          }
     258        }
     259
     260        if (!path_found) {
     261          p_it = paths.begin();
     262          n_it = path_ends.begin();
     263          while (p_it != paths.end()) {
     264            Node blue = _bpg.blueNode((*_matching)[*n_it]);
     265            reached_node[blue] = false;
     266            *n_it = blue;
     267            p_it->push_front((*_matching)[*n_it]);
     268            ++n_it;
     269            ++p_it;
     270          }
     271        } else {
     272          // Now 'paths' contains a maximal set of shortest augmenting paths.
     273          for (p_it = paths.begin(); p_it != paths.end(); ++p_it) {
     274            for (typename std::list<Edge>::iterator it = p_it->begin();
     275                 it != p_it->end();) {
     276              _matching->set(_bpg.redNode(*it), *it);
     277              _matching->set(_bpg.blueNode(*it), *it);
     278              ++it;
     279              if (it != p_it->end()) ++it;
     280            }
     281          }
     282        }
     283      }
     284
     285      return true;
     286    }
     287
     288    /// \brief Executes the algorithm
     289    ///
     290    /// It runs augmenting phases until the optimal solution is reached.
     291    ///
     292    /// \pre \ref init() or \ref matchingInit() must be called before using
     293    /// this function.
     294    void start() {
     295      while (augment()) {}
     296    }
     297
     298    /// \brief Runs the algorithm.
     299    ///
     300    /// hk.run() is just a shorthand for:
     301    ///
     302    ///\code
     303    /// hk.init();
     304    /// hk.start();
     305    ///\endcode
     306    void run() {
     307      init();
     308      start();
     309    }
     310
     311    /// \brief Size of the matching
     312    ///
     313    /// Returns the size of the current matching.
     314    int matchingSize() const {
     315      int size = 0;
     316      for (RedIt it(_bpg); it!=INVALID; ++it) {
     317        if ((*_matching)[it] != INVALID) ++size;
     318      }
     319      return size;
     320    }
     321
     322    /// \brief Return \c true if the given edge is in the matching.
     323    ///
     324    /// This function returns \c true if the given edge is in the current
     325    /// matching.
     326    bool matching(const Edge& edge) const {
     327      return edge == (*_matching)[_bpg.redNode(edge)];
     328    }
     329
     330    /// \brief Return the matching edge incident to the given node.
     331    ///
     332    /// This function returns the matching edge incident to the
     333    /// given node in the current matching or \c INVALID if the node is
     334    /// not covered by the matching.
     335    Edge matching(const Node& n) const {
     336      return (*_matching)[n];
     337    }
     338
     339    /// \brief Return the mate of the given node.
     340    ///
     341    /// This function returns the mate of the given node in the current
     342    /// matching or \c INVALID if the node is not covered by the matching.
     343    Node mate(const Node& n) const {
     344      return (*_matching)[n] != INVALID ?
     345        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
     346    }
     347  };
     348
     349}
     350#endif
  • test/CMakeLists.txt

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

    diff --git a/test/Makefile.am b/test/Makefile.am
    a b  
    3131        test/graph_utils_test \
    3232        test/hao_orlin_test \
    3333        test/heap_test \
     34        test/hopcroft_karp_test \
    3435        test/kruskal_test \
    3536        test/maps_test \
    3637        test/matching_test \
     
    8586test_heap_test_SOURCES = test/heap_test.cc
    8687test_kruskal_test_SOURCES = test/kruskal_test.cc
    8788test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
     89test_hopcroft_karp_test_SOURCES = test/hopcroft_karp_test.cc
    8890test_lp_test_SOURCES = test/lp_test.cc
    8991test_maps_test_SOURCES = test/maps_test.cc
    9092test_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