COIN-OR::LEMON - Graph Library

Ticket #168: bpm_upgrade.2.patch

File bpm_upgrade.2.patch, 81.2 KB (added by Poroszkai Daniel, 9 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 b  
    22 *
    33 * This file is a part of LEMON, a generic C++ optimization library.
    44 *
    5  * Copyright (C) 2003-2011
     5 * Copyright (C) 2003-2012
    66 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    77 * (Egervary Research Group on Combinatorial Optimization, EGRES).
    88 *
     
    2020#define LEMON_HOPCROFT_KARP_H
    2121
    2222#include <lemon/core.h>
     23#include <vector>
     24#include <queue>
    2325#include <list>
    2426
    2527/// \ingroup matching
     
    3537  template <typename BPG>
    3638  class HopcroftKarp {
    3739  public:
    38       /// Type of the bipartite graph
     40    /// Type of the bipartite graph
    3941    typedef BPG BpGraph;
    4042    /// Type of the matching map
    4143    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
     
    8183
    8284    /// \brief Sets the matching map
    8385    ///
    84     /// Sets the matching map.
     86    /// Sets the matching map. It should contain a valid matching,
     87    /// for example the empty matching is fine.
    8588    /// If you don't use this function before calling \ref init(),
    8689    /// an instance will be allocated automatically.
    8790    /// The destructor deallocates this automatically allocated map,
    8891    /// of course.
    89     /// This member is not to initialize the algorithm with a valid
    90     /// matching; use \ref matchingInit() instead.
    9192    ///
    9293    /// \return <tt>(*this)</tt>
    9394    HopcroftKarp& matchingMap(MatchingMap& map) {
     
    108109
    109110    /// \brief Initializes the algorithm
    110111    ///
    111     /// Allocates the matching map if necessary, and sets
    112     /// to empty matching.
     112    /// Allocates the matching map if necessary.
    113113    ///
    114114    /// \pre \ref init() is not called.
    115115    void init() {
    116116      createStructures();
    117       if (!_local_matching) {
    118         for (NodeIt it(_bpg); it!=INVALID; ++it) {
    119           _matching->set(it, INVALID);
    120         }
    121       }
    122117    }
    123118
    124119    /// \brief Smarter initialization of the matching.
     
    131126    ///
    132127    /// \note heuristicInit() replaces init() regarding the preconditions.
    133128    int heuristicInit() {
    134       init();
     129      createStructures();
    135130      int size = 0;
    136131
    137       for (BlueNodeIt b_it(_bpg); b_it!=INVALID; ++b_it) {
    138         if ((*_matching)[b_it] == INVALID) {
    139           bool matched = false;
    140           for (IncEdgeIt inc(_bpg, b_it); inc != INVALID && !matched; ++inc) {
    141             if ((*_matching)[_bpg.redNode(inc)] == INVALID) {
    142               _matching->set(_bpg.redNode(inc), inc);
    143               _matching->set(b_it, inc);
    144               matched = true;
    145               ++size;
    146             }
     132      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     133        bool matched = false;
     134        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
     135          if ((*_matching)[_bpg.target(oi)] == INVALID) {
     136            _matching->set(_bpg.source(oi), oi);
     137            _matching->set(_bpg.target(oi), oi);
     138            matched = true;
     139            ++size;
    147140          }
    148141        }
    149142      }
     143
    150144      return size;
    151145    }
    152146
     
    157151    /// in which an edge is true if it is in the matching.
    158152    ///
    159153    /// If the given matching is invalid, some edges will be left out.
    160     /// 
     154    ///
    161155    /// \return \c false if the given matching is invalid.
    162156    ///
    163157    /// \note matchingInit() replaces init() regarding the preconditions.
    164158    bool matchingInit(const BoolEdgeMap& matching) {
    165       init();
     159      createStructures();
    166160      bool valid = true;
    167161      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
    168162        if (matching[it]) {
     
    200194    ///
    201195    /// \pre \ref init() must be called before using this function.
    202196    int augment() {
    203       IntNodeMap level(_bpg, -1);
    204       std::list<RedNode> act_rednodes;
    205       std::list<BlueNode> act_bluenodes;
     197      IntBlueNodeMap level(_bpg, -2);
     198      std::vector<RedNode> act_red_nodes;
    206199
    207200      for (RedNodeIt r_it(_bpg); r_it != INVALID; ++r_it) {
    208201        if ((*_matching)[r_it] == INVALID) {
    209           act_rednodes.push_front(r_it);
     202          act_red_nodes.push_back(r_it);
    210203        }
    211204      }
    212 
    213       // Raise this flag when a shortest augmenting path is found.
    214       bool path_found = false;
    215       // This nodelist will contain the end nodes of the possible
     205      // This vector will contain the end nodes of the possible
    216206      // augmenting paths.
    217       std::list<BlueNode> path_ends;
     207      std::vector<BlueNode> path_ends;
    218208
    219209      // Layer counter
    220       int cur_level = 0;
     210      int cur_level = -1;
    221211
    222212      // Starting from the unmatched red nodes search for unmatched
    223213      // blue nodes, using Bfs (but only on alternating paths).
    224       while (!path_found) {
    225         while (!act_rednodes.empty()) {
    226           RedNode red = act_rednodes.front();
    227           act_rednodes.pop_front();
    228           level[red] = cur_level;
    229           for (IncEdgeIt it(_bpg, red); it!=INVALID; ++it) {
    230             BlueNode blue(_bpg.blueNode(it));
    231             if (level[blue] == -1) {
    232               act_bluenodes.push_front(blue);
    233               level[blue] = cur_level + 1;
    234               path_found |= ((*_matching)[blue] == INVALID);
     214      // Blue nodes with level==-2 are unreached.
     215      while (path_ends.empty()) {
     216        cur_level += 2;
     217        std::vector<RedNode> next_red_nodes;
     218        for (typename std::vector<RedNode>::iterator r=act_red_nodes.begin();
     219             r != act_red_nodes.end(); ++r) {
     220          for (OutArcIt it(_bpg, *r); it!=INVALID; ++it) {
     221            BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     222            if (level[blue] != -2) continue;
     223            level[blue] = cur_level;
     224            if ((*_matching)[blue] == INVALID) {
     225              path_ends.push_back(blue);
     226            } else {
     227              next_red_nodes
     228                .push_back(_bpg.redNode((*_matching)[blue]));
    235229            }
    236230          }
    237231        }
    238232
    239         if (!path_found) {
    240           if (act_bluenodes.empty()) return -1;
    241           while (!act_bluenodes.empty()) {
    242             BlueNode blue = act_bluenodes.front();
    243             act_bluenodes.pop_front();
    244             RedNode red = _bpg.redNode((*_matching)[blue]);
    245             act_rednodes.push_front(red);
    246           }
    247         } else {
    248           for (typename std::list<BlueNode>::iterator it=act_bluenodes.begin();
    249                it != act_bluenodes.end();
    250                ++it) {
    251             if ((*_matching)[*it] == INVALID) path_ends.push_front(*it);
    252           }
    253           // Now path_ends contains nodes that are ends of
    254           // shortest alternating paths.
    255         }
    256         cur_level += 2;
     233        if (path_ends.empty() && next_red_nodes.empty())
     234          return -1;
     235
     236        act_red_nodes.swap(next_red_nodes);
     237        next_red_nodes.clear();
    257238      }
    258 
    259       // The current_edge structure assures that edges iterated at most once
    260       typename BpGraph::template BlueNodeMap<IncEdgeIt*> current_edge(_bpg, 0);
     239      // length of found augmenting paths
     240      int max_level = cur_level;
    261241
    262242      // Stack for the Dfs
    263       std::list<BlueNode> stack;
    264 
     243      std::vector<OutArcIt> stack;
     244      // Raise this flag when a shortest augmenting path is found.
     245      bool path_found;
    265246      // Searching backward, starting from the previously found
    266247      // blue nodes, we build vertex disjoint alternating paths.
    267       typename std::list<BlueNode>::iterator n_it = path_ends.begin();
     248      // Blue nodes with level!=-2 are unreached.
     249      typename std::vector<BlueNode>::iterator end = path_ends.begin();
    268250
    269       while (n_it != path_ends.end()) {
    270         stack.push_front(*n_it);
     251      while (end != path_ends.end()) {
     252        stack.clear();
     253        stack.push_back(OutArcIt(_bpg, *end));
     254        cur_level = max_level;
    271255        path_found = false;
     256
    272257        while (!stack.empty() && !path_found) {
    273           BlueNode b = stack.front();
    274           if (current_edge[b] == 0) {
    275             current_edge[b] = new IncEdgeIt(_bpg, b);
     258          OutArcIt &o = stack.back();
     259          // Search an arc that goes one level higher.
     260          // If none found, retreat.
     261          while (o != INVALID &&
     262                 (*_matching)[_bpg.target(o)] != INVALID &&
     263                 level[_bpg.blueNode((*_matching)[_bpg.target(o)])] !=
     264                 cur_level - 2) {
     265            ++o;
     266          }
     267          if (o==INVALID) {
     268            stack.pop_back();
     269            cur_level += 2;
     270          } else if ((*_matching)[_bpg.target(o)] == INVALID) {
     271            path_found = true;
    276272          } else {
    277             ++(*current_edge[b]);
    278           }
    279           while (*current_edge[b] != INVALID &&
    280                  level[_bpg.redNode(*current_edge[b])] != level[b] - 1) {
    281             ++(*current_edge[b]);
    282           }
    283           if (*current_edge[b] == INVALID) {
    284             stack.pop_front();
    285           } else {
    286             RedNode r = _bpg.redNode(*current_edge[b]);
    287             if ((*_matching)[r] == INVALID) {
    288               path_found = true;
    289             } else {
    290               level[r] = -1; // Do not visit this node again
    291               stack.push_front(_bpg.blueNode((*_matching)[r]));
    292             }
     273            BlueNode b = _bpg.blueNode((*_matching)[_bpg.target(o)]);
     274            stack.push_back(OutArcIt(_bpg, b));
     275            level[b] = -2;        // Do not visit this node again
     276            cur_level -= 2;
    293277          }
    294278        }
    295279        if (path_found) {
    296           BlueNode next = *n_it;
    297           RedNode r;
    298           BlueNode b;
    299           Edge e;
    300           while (next != INVALID) {
    301             e = *current_edge[next];
    302             b = next;
    303             r = _bpg.redNode(e);
    304             next = (*_matching)[r] == INVALID ?
    305               INVALID :_bpg.blueNode((*_matching)[r]);
    306             _matching->set(b, e);
    307             _matching->set(r, e);
     280          OutArcIt o;
     281          while (!stack.empty()) {
     282            o = stack.back();
     283            _matching->set(_bpg.source(o), o);
     284            _matching->set(_bpg.target(o), o);
     285            stack.pop_back();
    308286          }
    309287        }
    310 
    311         stack.clear();
    312         ++n_it;
     288        ++end;
    313289      }
    314       return cur_level - 1;
     290      return max_level;
    315291    }
    316292
    317293    /// \brief Executes the algorithm
     
    339315    /// \brief Size of the matching
    340316    ///
    341317    /// Returns the size of the current matching.
     318    /// Function runs in \f$O(n)\f$ time.
    342319    int matchingSize() const {
    343320      int size = 0;
    344321      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
     
    372349      return (*_matching)[n] != INVALID ?
    373350        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
    374351    }
     352
     353    /// \brief Query a vertex cover.
     354    ///
     355    /// This function finds a vertex cover based on the current matching,
     356    /// and in \c cover sets node in the cover \c true and \c false
     357    /// otherwise. In case of maximal matching, this will be a minimal
     358    /// vertex cover.
     359    /// Function runs in \f$O(n + e)\f$ time.
     360    ///
     361    /// \return The number of nodes in the cover.
     362    int cover(BoolNodeMap &cover) const {
     363      int count = 0;
     364      std::queue<RedNode> q;
     365
     366      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
     367        cover[b] = false;
     368
     369      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     370        if ((*_matching)[r] == INVALID) {
     371          cover[r] = false;
     372          q.push(r);
     373        } else {
     374          cover[r] = true;
     375          ++count;
     376        }
     377      }
     378
     379      while (!q.empty()) {
     380        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     381          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     382          if (cover[blue]) continue;
     383
     384          cover[blue] = true;
     385          ++count;
     386
     387          if ((*_matching)[blue] != INVALID) {
     388            RedNode nr = _bpg.redNode((*_matching)[blue]);
     389            q.push(nr);
     390            cover[nr] = false;
     391            --count;
     392          }
     393        }
     394        q.pop();
     395      }
     396
     397      return count;
     398    }
     399
     400    /// \brief Query a barrier of red nodes.
     401    ///
     402    /// This function tries to create a barrier of red nodes, which should:
     403    ///  - contain every unmatched red nodes,
     404    ///  - contain exactly that many matched red nodes as much blue nodes
     405    ///    there are in its neighbor set.
     406    ///
     407    /// Barrier exist if and only if the matching is maximal. If exist, the
     408    /// map is set true for nodes of the barrier, and false for the others.
     409    ///
     410    /// Function runs in \f$O(n + e)\f$ time.
     411    ///
     412    /// \return \c true if a barrier found.
     413    bool redBarrier(BoolRedNodeMap &barrier) const {
     414      bool exist = true;
     415      std::queue<RedNode> q;
     416
     417      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     418        if ((*_matching)[r] == INVALID) {
     419          barrier[r] = true;
     420          q.push(r);
     421        } else {
     422          barrier[r] = false;
     423        }
     424      }
     425
     426      while (!q.empty() && exist) {
     427        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     428          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     429
     430          if (exist &= ((*_matching)[blue] != INVALID)) {
     431            RedNode r = _bpg.redNode((*_matching)[blue]);
     432            if (!barrier[r]) {
     433              q.push(r);
     434              barrier[r] = true;
     435            }
     436          }
     437        }
     438        q.pop();
     439      }
     440
     441      return exist;
     442    }
     443
     444    /// \brief Query a barrier of blue nodes.
     445    ///
     446    /// This function tries to create a barrier of blue nodes, which should:
     447    ///  - contain every unmatched blue nodes,
     448    ///  - contain exactly that many matched blue nodes as much red nodes
     449    ///    there are in its neighbor set.
     450    ///
     451    /// Barrier exist if and only if the matching is maximal. If exist, the
     452    /// map is set true for nodes of the barrier, and false for the others.
     453    ///
     454    /// Function runs in \f$O(n + e)\f$ time.
     455    ///
     456    /// \return \c true if a barrier found.
     457    bool blueBarrier(BoolBlueNodeMap &barrier) const {
     458      bool exist = true;
     459      std::queue<BlueNode> q;
     460
     461      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     462        if ((*_matching)[b] == INVALID) {
     463          barrier[b] = true;
     464          q.push(b);
     465        } else {
     466          barrier[b] = false;
     467        }
     468      }
     469
     470      while (!q.empty() && exist) {
     471        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     472          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
     473
     474          if (exist &= ((*_matching)[r] != INVALID)) {
     475            BlueNode b = _bpg.blueNode((*_matching)[r]);
     476            if (!barrier[b]) {
     477              q.push(b);
     478              barrier[b] = true;
     479            }
     480          }
     481        }
     482        q.pop();
     483      }
     484
     485      return exist;
     486    }
    375487  };
    376488
    377489}
  • test/hopcroft_karp_test.cc

    diff --git a/test/hopcroft_karp_test.cc b/test/hopcroft_karp_test.cc
    a b  
    22 *
    33 * This file is a part of LEMON, a generic C++ optimization library.
    44 *
    5  * Copyright (C) 2003-2011
     5 * Copyright (C) 2003-2012
    66 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    77 * (Egervary Research Group on Combinatorial Optimization, EGRES).
    88 *
     
    104104  BoolEdgeMap init_matching(graph);
    105105  HK hk_test(graph);
    106106  const HK& const_hk_test = hk_test;
    107   HK::MatchingMap matching_map(graph);
     107  HK::MatchingMap matching_map(graph, INVALID);
     108  BoolNodeMap cov(graph);
     109  int cover_size;
     110  BoolRedNodeMap rb(graph);
     111  BoolBlueNodeMap bb(graph);
    108112
    109113  hk_test.matchingMap(matching_map);
    110114  hk_test.init();
     
    123127  const_hk_test.mate(node);
    124128  const HK::MatchingMap& max_matching = const_hk_test.matchingMap();
    125129  edge = max_matching[node];
     130  cover_size = const_hk_test.cover(cov);
     131  const_hk_test.redBarrier(rb);
     132  const_hk_test.blueBarrier(bb);
    126133}
    127134
    128135typedef SmartBpGraph BpGraph;
     
    144151
    145152}
    146153
    147 void checkMatchingOptimality(const BpGraph& bpg, const HK& hk) {
    148   list<RedNode> reds;
    149   list<BlueNode> blues;
    150   BoolBlueNodeMap reached(bpg, false);
     154void checkCover(const BpGraph& bpg, const HK& hk) {
     155  int cov_size;
     156  BoolNodeMap cov(bpg);
     157  cov_size = hk.cover(cov);
     158  BoolEdgeMap covered(bpg, false);
     159  for (NodeIt n(bpg); n!=INVALID; ++n) {
     160    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
     161      covered[e] = true;
     162    }
     163  }
     164  bool all = true;
     165  for (EdgeIt e(bpg); e!=INVALID; ++e) {
     166    all&=covered[e];
     167  }
    151168
    152   for (RedNodeIt it(bpg); it != INVALID; ++it) {
    153     if (hk.matching(it) == INVALID) {
    154       reds.push_front(it);
     169  check(all, "Invalid cover");
     170  check(hk.matchingSize() == cov_size, "Suboptimal matching.");
     171}
     172
     173void checkBarriers(const BpGraph& bpg, const HK& hk) {
     174  BoolRedNodeMap rb(bpg);
     175  BoolBlueNodeMap bb(bpg);
     176
     177  check(hk.redBarrier(rb), "Suboptimal matching.");
     178  check(hk.blueBarrier(bb), "Suboptimal matching.");
     179
     180  BoolNodeMap reached(bpg, false);
     181
     182  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     183    if (!rb[r]) continue;
     184    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
     185      reached.set(bpg.target(a), true);
     186    }
     187  }
     188  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     189    if (!bb[b]) continue;
     190    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
     191      reached.set(bpg.target(a), true);
    155192    }
    156193  }
    157194
    158   while (!reds.empty()) {
    159     while (!reds.empty()) {
    160       RedNode red = reds.front();
    161       reds.pop_front();
    162       for (IncEdgeIt it(bpg, red); it != INVALID; ++it) {
    163         BlueNode blue = bpg.blueNode(it);
    164         if (!reached[blue]) {
    165           blues.push_front(blue);
    166           reached.set(blue, true);
    167         }
    168       }
    169     }
    170 
    171     while (!blues.empty()) {
    172       BlueNode blue = blues.front();
    173       blues.pop_front();
    174       check(hk.matching(blue) != INVALID, "Suboptimal matching");
    175       reds.push_front(bpg.asRedNode(hk.mate(blue)));
    176     }
     195  int ur = 0, ub = 0, rc = 0, bc = 0;
     196  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     197    if (rb[r]) ++rc;
     198    if (hk.mate(r) == INVALID) ++ur;
     199    if (reached[r]) --bc;
    177200  }
    178201
     202  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     203    if (bb[b]) ++bc;
     204    if (hk.mate(b) == INVALID) ++ub;
     205    if (reached[b]) --rc;
     206  }
     207
     208  check(rc == ur, "Red barrier is wrong.");
     209  check(bc == ub, "Blue barrier is wrong.");
    179210}
    180211
     212void checkMatching(const BpGraph& bpg, const HK& hk) {
     213  checkMatchingValidity(bpg, hk);
     214  checkBarriers(bpg, hk);
     215  checkCover(bpg, hk);
     216}
     217
     218
    181219int main() {
    182220  // Checking hard wired extremal graphs
    183221  for (int i=0; i<lgfn; ++i) {
     
    186224    bpGraphReader(bpg, lgfs).run();
    187225    HK hk(bpg);
    188226    hk.run();
    189     checkMatchingValidity(bpg, hk);
    190     checkMatchingOptimality(bpg, hk);
     227    checkMatching(bpg, hk);
    191228  }
    192229
    193230  // Checking some use cases
     
    199236    hk.init();
    200237    hk.augment();
    201238    hk.start();
    202     checkMatchingValidity(bpg, hk);
    203     checkMatchingOptimality(bpg, hk);
     239    checkMatching(bpg, hk);
    204240  }
    205241  {
    206242    HK hk(bpg);
    207243    hk.heuristicInit();
    208244    hk.augment();
    209245    hk.start();
    210     checkMatchingValidity(bpg, hk);
    211     checkMatchingOptimality(bpg, hk);
     246    checkMatching(bpg, hk);
    212247  }
    213248  {
    214249    HK hk(bpg);
    215     HK::MatchingMap m(bpg);
     250    HK::MatchingMap m(bpg, INVALID);
    216251    hk.matchingMap(m);
    217252    hk.run();
    218     checkMatchingValidity(bpg, hk);
    219     checkMatchingOptimality(bpg, hk);
     253    checkMatching(bpg, hk);
    220254  }
    221255  {
    222256    HK hk(bpg);
     
    226260    init_matching[e] = true;
    227261    check(hk.matchingInit(init_matching), "Wrong result from matchingInit()");
    228262    hk.start();
    229     checkMatchingValidity(bpg, hk);
    230     checkMatchingOptimality(bpg, hk);
     263    checkMatching(bpg, hk);
    231264  }
    232265  {
    233266    HK hk(bpg);
    234     HK::MatchingMap m(bpg);
     267    HK::MatchingMap m(bpg, INVALID);
    235268    hk.matchingMap(m);
    236269    BpGraph::EdgeMap<bool> init_matching(bpg, true);
    237270    check(!hk.matchingInit(init_matching), "Wrong result from matchingInit()");
    238271    hk.start();
    239     checkMatchingValidity(bpg, hk);
    240     checkMatchingOptimality(bpg, hk);
     272    checkMatching(bpg, hk);
    241273  }
    242 
    243274  // Checking some random generated graphs
    244275  const int random_test_n = 10;
    245   const int max_red_n = 100;
     276  const int max_red_n = 1000;
    246277  const int min_red_n = 10;
    247   const int max_blue_n = 100;
     278  const int max_blue_n = 1000;
    248279  const int min_blue_n = 10;
    249280  for (int i=0; i<random_test_n; ++i) {
    250281    int red_n = rnd.integer(min_red_n, max_red_n);
     
    259290    }
    260291    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
    261292      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
    262         if (rnd() < prob) {
     293        if (rnd() < prob/10.) {
    263294          bpg.addEdge(r,b);
    264295        }
    265296      }
     
    267298
    268299    HK hk(bpg);
    269300    hk.run();
    270     checkMatchingValidity(bpg, hk);
    271     checkMatchingOptimality(bpg, hk);
     301    checkMatching(bpg, hk);
    272302  }
    273303
    274304
  • 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 b  
    22 *
    33 * This file is a part of LEMON, a generic C++ optimization library.
    44 *
    5  * Copyright (C) 2003-2009
     5 * Copyright (C) 2003-2012
    66 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    77 * (Egervary Research Group on Combinatorial Optimization, EGRES).
    88 *
     
    108108    Elevator(const GR &graph,int max_level) :
    109109      _g(graph),
    110110      _max_level(max_level),
    111       _item_num(_max_level),
     111      _item_num(countItems<GR, Item>(graph)),
    112112      _where(graph),
    113113      _level(graph,0),
    114       _items(_max_level),
     114      _items(_item_num),
    115115      _first(_max_level+2),
    116116      _last_active(_max_level+2),
    117117      _highest_active(-1) {}
     
    526526    ///\param max_level The maximum allowed level.
    527527    ///Set the range of the possible labels to <tt>[0..max_level]</tt>.
    528528    LinkedElevator(const GR& graph, int max_level)
    529       : _graph(graph), _max_level(max_level), _item_num(_max_level),
     529      : _graph(graph), _max_level(max_level),
     530        _item_num(countItems<GR, Item>(graph)),
    530531        _first(_max_level + 1), _last(_max_level + 1),
    531532        _prev(graph), _next(graph),
    532533        _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
    - +  
     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-2012
     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_PREFLOW_BP_MATCHING_H
     20#define LEMON_PREFLOW_BP_MATCHING_H
     21
     22#include <lemon/core.h>
     23#include <lemon/elevator.h>
     24#include <vector>
     25#include <list>
     26#include <queue>
     27
     28/// \ingroup matching
     29/// \file
     30/// \brief A preflow based bipartite matching algorithm.
     31namespace lemon {
     32
     33  /// \brief A preflow based bipartite matching algorithm
     34  ///
     35  /// Finds maximum cardinality matching in a given bipartite
     36  /// graph using specialized preflow algorithm.
     37  template <typename BPG>
     38  class PreflowBpMatching {
     39  public:
     40    /// Type of the bipartite graph
     41    typedef BPG BpGraph;
     42    /// Type of the matching map
     43    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
     44    /// Type of the elevator.
     45    typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t;
     46
     47  private:
     48    TEMPLATE_BPGRAPH_TYPEDEFS(BPG);
     49
     50  protected:
     51    typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow;
     52    typedef std::list<typename BPG::Arc> CurrentArcSet;
     53    typedef typename
     54    std::list<typename BPG::Arc>::iterator CurrentArcIterator;
     55
     56    const BpGraph& _bpg;
     57    MatchingMap* _matching;
     58    bool _local_matching;
     59    Elevator_t* _elevator;
     60    bool _local_elevator;
     61    typename BPG::template BlueNodeMap<CurrentArcSet> _cas;
     62    typename BPG::template BlueNodeMap<CurrentArcIterator> _ca;
     63    Preflow _flow;
     64
     65    void createStructures() {
     66      if (!_matching) {
     67        _matching = new MatchingMap(_bpg, INVALID);
     68        _local_matching = true;
     69      }
     70      if (!_elevator) {
     71        _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1,
     72                                                  countBlueNodes(_bpg)));
     73        _local_elevator = true;
     74      }
     75    }
     76
     77    void destroyStructures() {
     78      if (_local_matching) {
     79        delete _matching;
     80      }
     81      if (_local_elevator) {
     82        delete _elevator;
     83      }
     84    }
     85
     86    /// \brief Initialize the elevator
     87    ///
     88    /// When this function is called, the elevator is initialized with
     89    /// the exact distance labels corresponding the current matching.
     90    void initElevator() {
     91
     92      // Set the initial flow:
     93      // choose arbitrary outgoing arc for every unmatched red node.
     94      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     95        Arc a = (*_matching)[r] != INVALID ?
     96          _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r);
     97        if (a != INVALID) {
     98          _flow[r] = a;
     99          _cas[_bpg.asBlueNodeUnsafe(_bpg.target(a))].push_front(a);
     100        }
     101      }
     102      // Now we start a backward Bfs from the unmatched blue nodes in the
     103      // residual graph, to supply the exact distance labels to the elevator.
     104      std::vector<BlueNode> act_blue_nodes, next_blue_nodes;
     105      BoolBlueNodeMap reached(_bpg, false);
     106
     107      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
     108        if (_cas[b].empty()) {
     109          act_blue_nodes.push_back(b);
     110          reached.set(b, true);
     111        }
     112      }
     113
     114      _elevator->initStart();
     115      while (!act_blue_nodes.empty()) {
     116        for (typename std::vector<BlueNode>::iterator b_it =
     117               act_blue_nodes.begin();
     118             b_it != act_blue_nodes.end(); ++b_it) {
     119
     120          _elevator->initAddItem(*b_it);
     121          for (OutArcIt oi(_bpg, *b_it); oi!=INVALID; ++oi) {
     122            RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(oi));
     123            BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r]));
     124            if (b != *b_it && !reached[b]) {
     125              next_blue_nodes.push_back(b);
     126              reached.set(b, true);
     127            }
     128          }
     129        }
     130        _elevator->initNewLevel();
     131        next_blue_nodes.swap(act_blue_nodes);
     132        next_blue_nodes.clear();
     133      }
     134      _elevator->initFinish();
     135
     136      // Activate blue nodes with excess, i. e. those
     137      // which get flow from more than one red node.
     138      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
     139        if (_cas[b].size() > 1 && (*_elevator)[b] < _elevator->maxLevel()) {
     140          _elevator->activate(b);
     141        }
     142
     143        _ca[b] = _cas[b].begin();
     144      }
     145    }
     146
     147    PreflowBpMatching() {}
     148
     149  public:
     150
     151    /// \brief Constructor
     152    ///
     153    /// Constructs the class on the given bipartite graph.
     154    PreflowBpMatching(const BpGraph& bpg) : _bpg(bpg),
     155                                            _matching(0),
     156                                            _local_matching(false),
     157                                            _elevator(0),
     158                                            _local_elevator(false),
     159                                            _cas(bpg),
     160                                            _ca(bpg),
     161                                            _flow(bpg, INVALID)
     162    {}
     163
     164    /// \brief Destructor
     165    ///
     166    /// Destructor.
     167    ~PreflowBpMatching() {
     168      destroyStructures();
     169    }
     170
     171    /// \brief Sets the matching map
     172    ///
     173    /// Sets the matching map. It should contain a valid matching,
     174    /// for example the empty matching is fine.
     175    /// If you don't use this function before calling \ref init(),
     176    /// an instance will be allocated automatically.
     177    /// The destructor deallocates this automatically allocated map,
     178    /// of course.
     179    ///
     180    /// \return <tt>(*this)</tt>
     181    PreflowBpMatching& matchingMap(MatchingMap& map) {
     182      _matching = &map;
     183      _local_matching = false;
     184      return *this;
     185    }
     186
     187    /// \brief Returns a const reference to the matching map.
     188    ///
     189    /// Returns a const reference to the matching map, which contains
     190    /// the matching edge for each node (or INVALID for unmatched nodes).
     191    const MatchingMap& matchingMap() const {
     192      return *_matching;
     193    }
     194
     195    /// \brief Sets the elevator used by algorithm.
     196    ///
     197    /// Sets the elevator.
     198    /// If you don't use this function before calling \ref init(),
     199    /// an instance will be allocated automatically.
     200    /// The destructor deallocates this automatically allocated elevator,
     201    /// of course.
     202    /// \return <tt>(*this)</tt>
     203    PreflowBpMatching& elevator(Elevator_t& elevator) {
     204      _elevator = &elevator;
     205      _local_elevator = false;
     206      return *this;
     207    }
     208
     209    /// \brief Returns a const reference to the elevator.
     210    ///
     211    /// Returns a const reference to the elevator, which keeps track
     212    /// of the distance labels of nodes.
     213    const Elevator_t& elevator() const {
     214      return *_elevator;
     215    }
     216
     217    /// \brief Initializes the algorithm
     218    ///
     219    /// Allocates the elevator and the matching map if necessary,
     220    /// and initializes the elevator.
     221    ///
     222    /// \pre \ref init() is not called.
     223    void init() {
     224      createStructures();
     225      initElevator();
     226    }
     227
     228    /// \brief Smarter initialization of the matching.
     229    ///
     230    /// Allocates the matching map if necessary, and initializes
     231    /// the algorithm with a trivially found matching: iterating
     232    /// on every edge, take it in if both endnodes are unmatched.
     233    ///
     234    /// \return Size of the so found matching.
     235    ///
     236    /// \note heuristicInit() replaces init() regarding the preconditions.
     237    int heuristicInit() {
     238      createStructures();
     239      int size = 0;
     240
     241      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     242        bool matched = false;
     243        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
     244          if ((*_matching)[_bpg.target(oi)] == INVALID) {
     245            _matching->set(_bpg.source(oi), oi);
     246            _matching->set(_bpg.target(oi), oi);
     247            matched = true;
     248            ++size;
     249          }
     250        }
     251      }
     252
     253      initElevator();
     254      return size;
     255    }
     256
     257    /// \brief Initialize the matching from a map.
     258    ///
     259    /// Allocates the matching map if necessary, and
     260    /// initializes the algorithm with a matching given as a bool edgemap,
     261    /// in which an edge is true if it is in the matching.
     262    /// Also initializes the elevator.
     263    ///
     264    /// \return \c true if the matching is valid.
     265    ///
     266    /// \note matchingInit() replaces init() regarding the preconditions.
     267    bool matchingInit(const BoolEdgeMap& matching) {
     268      createStructures();
     269
     270      bool valid = true;
     271      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
     272        if (matching[it]) {
     273          Node red = _bpg.redNode(it);
     274          Node blue = _bpg.blueNode(it);
     275          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
     276            valid = false;
     277          } else {
     278            _matching->set(red, it);
     279            _matching->set(blue, it);
     280          }
     281        }
     282      }
     283
     284      initElevator();
     285      return valid;
     286    }
     287
     288
     289    /// \brief Perform a push/relabel operation on a node
     290    ///
     291    /// Perform a push/relabel operation on a node.
     292    ///
     293    /// \pre \ref init() is called, and the node is active.
     294    void push_relabel(const BlueNode& chosen) {
     295      Arc new_flow = INVALID;
     296      while (_ca[chosen] != _cas[chosen].end() && new_flow == INVALID) {
     297        OutArcIt r_oi(_bpg, _bpg.source(*_ca[chosen]));
     298        while (r_oi != INVALID && new_flow == INVALID) {
     299          if ((*_elevator)[_bpg.asBlueNodeUnsafe(_bpg.target(r_oi))] ==
     300              (*_elevator)[chosen] - 1) {
     301            new_flow = r_oi;
     302          }
     303          ++r_oi;
     304        }
     305        if (new_flow == INVALID) ++_ca[chosen];
     306      }
     307
     308      // If an incident blue node is found on one level lower,
     309      // then perform a push, i. e. flip the flow from the
     310      // current arc to 'new_flow'.
     311      if (new_flow != INVALID) {
     312        _ca[chosen] = _cas[chosen].erase(_ca[chosen]);
     313        if (_cas[chosen].size() == 1) {
     314          _elevator->deactivate(chosen);
     315        }
     316        _flow[_bpg.asRedNodeUnsafe(_bpg.source(new_flow))] = new_flow;
     317        BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(new_flow));
     318        _cas[b].push_front(new_flow);
     319        if (_cas[b].size() == 2) {
     320          _elevator->activate(b);
     321        }
     322      }
     323
     324      // If none found, then relabel.
     325      else {
     326        int min = _elevator->maxLevel() + 1;
     327        for (CurrentArcIterator it1 = _cas[chosen].begin();
     328             it1 != _cas[chosen].end();
     329             ++it1) {
     330          for (OutArcIt it2(_bpg, _bpg.asRedNodeUnsafe(_bpg.source(*it1)));
     331               it2 != INVALID;
     332               ++it2) {
     333            BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(it2));
     334            if (b != chosen && min > (*_elevator)[b]) {
     335              min = (*_elevator)[b];
     336            }
     337          }
     338        }
     339
     340        int chosen_level = (*_elevator)[chosen];
     341        if (min >= _elevator->maxLevel() - 1) {
     342          _elevator->lift(chosen, _elevator->maxLevel());
     343          _elevator->deactivate(chosen);
     344        } else {
     345          _elevator->lift(chosen, min+1);
     346        }
     347        if (_elevator->emptyLevel(chosen_level)) {
     348          _elevator->liftToTop(chosen_level);
     349        }
     350
     351        // Restore the current arc iterator.
     352        _ca[chosen] = _cas[chosen].begin();
     353      }
     354    }
     355
     356    /// \brief Executes the algorithm
     357    ///
     358    /// Executes push/relabel operation on the highest active node
     359    /// until there is. After start(), writeMatching() should be called to
     360    /// remove unnecessary flow and choose a valid maximum matching.
     361    ///
     362    /// \pre \ref init() must be called before using this function.
     363    void start() {
     364      while (_elevator->highestActive() != INVALID)
     365        push_relabel(_elevator->highestActive());
     366    }
     367
     368    /// \brief Deduces and writes a matching from the preflow
     369    ///
     370    /// The underlying preflow algorithm does not use the matching map,
     371    /// so before using the matching query functions ( \ref mate(), \ref
     372    /// matchingSize(), etc.), you should explicitly call this function
     373    /// to write a matching based on the current preflow to the matching
     374    /// map. When the preflow results an ambiguous matching (a node has
     375    /// incoming flow on more than one arc), it chooses arbitrary one.
     376    void writeMatching() {
     377      for (NodeIt it(_bpg); it != INVALID; ++it) {
     378        _matching->set(it, INVALID);
     379      }
     380      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
     381        if (!_cas[b].empty()) {
     382          Arc a = _cas[b].front();
     383          _matching->set(_bpg.source(a), a);
     384          _matching->set(_bpg.target(a), a);
     385        }
     386      }
     387    }
     388
     389    /// \brief Runs the algorithm.
     390    ///
     391    /// bp1.run() is just a shorthand for:
     392    ///
     393    ///\code
     394    /// bp1.heuristicInit();
     395    /// bp1.start();
     396    /// bp1.writeMatching();
     397    ///\endcode
     398    void run() {
     399      heuristicInit();
     400      start();
     401      writeMatching();
     402    }
     403
     404    /// \brief Size of the matching
     405    ///
     406    /// Returns the size of the current matching.
     407    /// Function runs in \f$O(n)\f$ time.
     408    int matchingSize() const {
     409      int size = 0;
     410      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
     411        if ((*_matching)[it] != INVALID) ++size;
     412      }
     413      return size;
     414    }
     415
     416    /// \brief Return \c true if the given edge is in the matching.
     417    ///
     418    /// This function returns \c true if the given edge is in the current
     419    /// matching.
     420    bool matching(const Edge& edge) const {
     421      return edge == (*_matching)[_bpg.redNode(edge)];
     422    }
     423
     424    /// \brief Return the matching edge incident to the given node.
     425    ///
     426    /// This function returns the matching edge incident to the
     427    /// given node in the current matching or \c INVALID if the node is
     428    /// not covered by the matching.
     429    Edge matching(const Node& n) const {
     430      return (*_matching)[n];
     431    }
     432
     433    /// \brief Return the mate of the given node.
     434    ///
     435    /// This function returns the mate of the given node in the current
     436    /// matching or \c INVALID if the node is not covered by the matching.
     437    Node mate(const Node& n) const {
     438      return (*_matching)[n] != INVALID ?
     439        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
     440    }
     441
     442    /// \brief Query a vertex cover.
     443    ///
     444    /// This function finds a vertex cover based on the current matching,
     445    /// and in \c cover sets node in the cover \c true and \c false
     446    /// otherwise. In case of maximal matching, this will be a minimal
     447    /// vertex cover.
     448    /// Function runs in \f$O(n + e)\f$ time.
     449    ///
     450    /// \return The number of nodes in the cover.
     451    int cover(BoolNodeMap &cover) const {
     452      int count = 0;
     453      std::queue<RedNode> q;
     454
     455      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
     456        cover[b] = false;
     457
     458      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     459        if ((*_matching)[r] == INVALID) {
     460          cover[r] = false;
     461          q.push(r);
     462        } else {
     463          cover[r] = true;
     464          ++count;
     465        }
     466      }
     467
     468      while (!q.empty()) {
     469        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     470          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     471          if (cover[blue]) continue;
     472
     473          cover[blue] = true;
     474          ++count;
     475
     476          if ((*_matching)[blue] != INVALID) {
     477            RedNode nr = _bpg.redNode((*_matching)[blue]);
     478            q.push(nr);
     479            cover[nr] = false;
     480            --count;
     481          }
     482        }
     483        q.pop();
     484      }
     485
     486      return count;
     487    }
     488
     489    /// \brief Query a barrier of red nodes.
     490    ///
     491    /// This function tries to create a barrier of red nodes, which should:
     492    ///  - contain every unmatched red nodes,
     493    ///  - contain exactly that many matched red nodes as much blue nodes
     494    ///    there are in its neighbor set.
     495    ///
     496    /// Barrier exist if and only if the matching is maximal. If exist, the
     497    /// map is set true for nodes of the barrier, and false for the others.
     498    ///
     499    /// Function runs in \f$O(n + e)\f$ time.
     500    ///
     501    /// \return \c true if a barrier found.
     502    bool redBarrier(BoolRedNodeMap &barrier) const {
     503      bool exist = true;
     504      std::queue<RedNode> q;
     505
     506      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     507        if ((*_matching)[r] == INVALID) {
     508          barrier[r] = true;
     509          q.push(r);
     510        } else {
     511          barrier[r] = false;
     512        }
     513      }
     514
     515      while (!q.empty() && exist) {
     516        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     517          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     518
     519          if (exist &= ((*_matching)[blue] != INVALID)) {
     520            RedNode r = _bpg.redNode((*_matching)[blue]);
     521            if (!barrier[r]) {
     522              q.push(r);
     523              barrier[r] = true;
     524            }
     525          }
     526        }
     527        q.pop();
     528      }
     529
     530      return exist;
     531    }
     532
     533    /// \brief Query a barrier of blue nodes.
     534    ///
     535    /// This function tries to create a barrier of blue nodes, which should:
     536    ///  - contain every unmatched blue nodes,
     537    ///  - contain exactly that many matched blue nodes as much red nodes
     538    ///    there are in its neighbor set.
     539    ///
     540    /// Barrier exist if and only if the matching is maximal. If exist, the
     541    /// map is set true for nodes of the barrier, and false for the others.
     542    ///
     543    /// Function runs in \f$O(n + e)\f$ time.
     544    ///
     545    /// \return \c true if a barrier found.
     546    bool blueBarrier(BoolBlueNodeMap &barrier) const {
     547      bool exist = true;
     548      std::queue<BlueNode> q;
     549
     550      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     551        if ((*_matching)[b] == INVALID) {
     552          barrier[b] = true;
     553          q.push(b);
     554        } else {
     555          barrier[b] = false;
     556        }
     557      }
     558
     559      while (!q.empty() && exist) {
     560        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     561          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
     562
     563          if (exist &= ((*_matching)[r] != INVALID)) {
     564            BlueNode b = _bpg.blueNode((*_matching)[r]);
     565            if (!barrier[b]) {
     566              q.push(b);
     567              barrier[b] = true;
     568            }
     569          }
     570        }
     571        q.pop();
     572      }
     573
     574      return exist;
     575    }
     576  };
     577
     578}
     579#endif
  • test/CMakeLists.txt

    diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
    a b  
    4848  path_test
    4949  planarity_test
    5050  preflow_test
     51  preflow_bp_matching_test
    5152  radix_sort_test
    5253  random_test
    5354  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
    - +  
     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-2012
     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/preflow_bp_matching.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
     71const string u_lgf =
     72  "@red_nodes\n"
     73  "label\n"
     74  "1\n"
     75  "2\n"
     76  "3\n"
     77  "@blue_nodes\n"
     78  "label\n"
     79  "4\n"
     80  "5\n"
     81  "6\n"
     82  "7\n"
     83  "@edges\n"
     84  "label\n"
     85  "1 4 1\n"
     86  "1 5 2\n"
     87  "2 5 3\n"
     88  "2 6 4\n"
     89  "2 7 5\n"
     90  "3 7 6\n"
     91;
     92
     93
     94
     95void checkPreflowBpMatchingCompile() {
     96  typedef concepts::BpGraph BpGraph;
     97  BPGRAPH_TYPEDEFS(BpGraph);
     98  typedef PreflowBpMatching<BpGraph> PBM;
     99
     100  BpGraph graph;
     101  RedNode red;
     102  BlueNode blue;
     103  Node node;
     104  Edge edge;
     105  BoolEdgeMap init_matching(graph);
     106  PBM pbm_test(graph);
     107  const PBM& const_pbm_test = pbm_test;
     108  PBM::MatchingMap matching_map(graph, INVALID);
     109  PBM::Elevator_t elevator(graph);
     110  BpGraph::NodeMap<bool> cov(graph);
     111  int cover_size;
     112  BoolRedNodeMap rb(graph);
     113  BoolBlueNodeMap bb(graph);
     114
     115  pbm_test.matchingMap(matching_map);
     116  pbm_test.elevator(elevator);
     117  pbm_test.init();
     118  pbm_test.heuristicInit();
     119  pbm_test.matchingInit(init_matching);
     120
     121  pbm_test.start();
     122  pbm_test.push_relabel(blue);
     123  pbm_test.run();
     124  pbm_test.writeMatching();
     125
     126  const_pbm_test.matchingSize();
     127  const_pbm_test.matching(node);
     128  const_pbm_test.matching(edge);
     129  const_pbm_test.mate(red);
     130  const_pbm_test.mate(blue);
     131  const_pbm_test.mate(node);
     132  const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap();
     133  edge = max_matching[node];
     134  const PBM::Elevator_t& elev = const_pbm_test.elevator();
     135  blue = elev.highestActive();
     136  cover_size = const_pbm_test.cover(cov);
     137  const_pbm_test.redBarrier(rb);
     138  const_pbm_test.blueBarrier(bb);
     139}
     140
     141typedef SmartBpGraph BpGraph;
     142typedef PreflowBpMatching<BpGraph> PBM;
     143BPGRAPH_TYPEDEFS(BpGraph);
     144
     145void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) {
     146  BoolBlueNodeMap matched(bpg, false);
     147  for (RedNodeIt it(bpg); it != INVALID; ++it) {
     148    if (pbm.matching(it) != INVALID) {
     149      check(pbm.mate(pbm.mate(it)) == it, "Wrong matching");
     150      matched.set(bpg.asBlueNode(pbm.mate(it)), true);
     151    }
     152  }
     153  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
     154    // != is used as logical xor here
     155    check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching");
     156  }
     157}
     158
     159void checkCover(const BpGraph& bpg, const PBM& pbm) {
     160  int cov_size;
     161  BoolNodeMap cov(bpg);
     162  cov_size = pbm.cover(cov);
     163  BoolEdgeMap covered(bpg, false);
     164  for (NodeIt n(bpg); n!=INVALID; ++n) {
     165    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
     166      covered[e] = true;
     167    }
     168  }
     169  bool all = true;
     170  for (EdgeIt e(bpg); e!=INVALID; ++e) {
     171    all&=covered[e];
     172  }
     173
     174  check(all, "Invalid cover");
     175  check(pbm.matchingSize() == cov_size, "Suboptimal matching.");
     176}
     177
     178void checkBarriers(const BpGraph& bpg, const PBM& pbm) {
     179  BoolRedNodeMap rb(bpg);
     180  BoolBlueNodeMap bb(bpg);
     181
     182  check(pbm.redBarrier(rb), "Suboptimal matching.");
     183  check(pbm.blueBarrier(bb), "Suboptimal matching.");
     184
     185  BoolNodeMap reached(bpg, false);
     186
     187  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     188    if (!rb[r]) continue;
     189    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
     190      reached.set(bpg.target(a), true);
     191    }
     192  }
     193  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     194    if (!bb[b]) continue;
     195    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
     196      reached.set(bpg.target(a), true);
     197    }
     198  }
     199
     200  int ur = 0, ub = 0, rc = 0, bc = 0;
     201  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     202    if (rb[r]) ++rc;
     203    if (pbm.mate(r) == INVALID) ++ur;
     204    if (reached[r]) --bc;
     205  }
     206
     207  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     208    if (bb[b]) ++bc;
     209    if (pbm.mate(b) == INVALID) ++ub;
     210    if (reached[b]) --rc;
     211  }
     212
     213  check(rc == ur, "Red barrier is wrong.");
     214  check(bc == ub, "Blue barrier is wrong.");
     215}
     216
     217void checkMatching(const BpGraph& bpg, const PBM& pbm) {
     218  checkMatchingValidity(bpg, pbm);
     219  checkBarriers(bpg, pbm);
     220  checkCover(bpg, pbm);
     221}
     222
     223
     224int main() {
     225  // Checking hard wired extremal graphs
     226  for (int i=0; i<lgfn; ++i) {
     227    BpGraph bpg;
     228    istringstream lgfs(lgf[i]);
     229    bpGraphReader(bpg, lgfs).run();
     230    PBM pbm(bpg);
     231    pbm.run();
     232    checkMatching(bpg, pbm);
     233  }
     234
     235  // Checking some use cases
     236  BpGraph bpg;
     237  istringstream u_lgfs(u_lgf);
     238  bpGraphReader(bpg, u_lgfs).run();
     239
     240  {
     241    PBM pbm(bpg);
     242    pbm.init();
     243    BpGraph::BlueNode b = pbm.elevator().highestActive();
     244    if (b != INVALID) {
     245      pbm.push_relabel(b);
     246    }
     247    pbm.start();
     248    pbm.writeMatching();
     249    check(pbm.elevator().highestActive() == INVALID, "Error in elevator");
     250    checkMatching(bpg, pbm);
     251  }
     252
     253  {
     254    PBM pbm(bpg);
     255    pbm.heuristicInit();
     256    BpGraph::BlueNode b = pbm.elevator().highestActive();
     257    if (b != INVALID) {
     258      pbm.push_relabel(b);
     259    }
     260    pbm.start();
     261    checkMatching(bpg, pbm);
     262  }
     263  {
     264    PBM pbm(bpg);
     265    PBM::MatchingMap m(bpg, INVALID);
     266    pbm.matchingMap(m);
     267    pbm.run();
     268    checkMatching(bpg, pbm);
     269  }
     270  {
     271    PBM pbm(bpg);
     272    BpGraph::EdgeMap<bool> init_matching(bpg, false);
     273    BpGraph::Edge e;
     274    bpg.first(e);
     275    init_matching[e] = true;
     276    check(pbm.matchingInit(init_matching),
     277          "Wrong result from matchingInit()");
     278    pbm.start();
     279    pbm.writeMatching();
     280    checkMatching(bpg, pbm);
     281  }
     282  {
     283    PBM pbm(bpg);
     284    PBM::MatchingMap m(bpg, INVALID);
     285    pbm.matchingMap(m);
     286    BpGraph::EdgeMap<bool> init_matching(bpg, true);
     287    check(!pbm.matchingInit(init_matching),
     288          "Wrong result from matchingInit()");
     289    pbm.start();
     290    pbm.writeMatching();
     291    checkMatching(bpg, pbm);
     292  }
     293
     294  // Checking some random generated graphs
     295  const int random_test_n = 10;
     296  const int max_red_n = 1000;
     297  const int min_red_n = 10;
     298  const int max_blue_n = 1000;
     299  const int min_blue_n = 10;
     300  for (int i=0; i<random_test_n; ++i) {
     301    int red_n = rnd.integer(min_red_n, max_red_n);
     302    int blue_n = rnd.integer(min_blue_n, max_blue_n);
     303    double prob = rnd();
     304
     305    SmartBpGraph bpg;
     306    for (int j=0; j<red_n; ++j) {
     307      bpg.addRedNode();
     308    }
     309    for (int j=0; j<blue_n; ++j) {
     310      bpg.addBlueNode();
     311    }
     312    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
     313      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
     314        if (rnd() < prob/10.) {
     315          bpg.addEdge(r,b);
     316        }
     317      }
     318    }
     319
     320    PBM pbm(bpg);
     321    pbm.run();
     322    checkMatching(bpg, pbm);
     323
     324    // do not always use the highest active node
     325    PBM pbm2(bpg);
     326    PBM::Elevator_t elev(bpg);
     327    pbm2.elevator(elev);
     328    pbm2.init();
     329    while (elev.highestActive() != INVALID) {
     330      for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     331        if (elev.active(b)) pbm2.push_relabel(b);
     332      }
     333    }
     334    pbm2.writeMatching();
     335    checkMatching(bpg, pbm2);
     336  }
     337
     338
     339  return 0;
     340}
     341
     342
  • 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
    - +  
     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-2012
     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_PREFLOW_BP_MATCHING_2_H
     20#define LEMON_PREFLOW_BP_MATCHING_2_H
     21
     22#include <lemon/core.h>
     23#include <lemon/elevator.h>
     24#include <vector>
     25#include <queue>
     26
     27/// \ingroup matching
     28/// \file
     29/// \brief A preflow based bipartite matching algorithm.
     30namespace lemon {
     31
     32  /// \brief A preflow based bipartite matching algorithm
     33  ///
     34  /// Finds maximum cardinality matching in a given bipartite
     35  /// graph using specialized preflow algorithm.
     36  template <typename BPG>
     37  class PreflowBpMatching2 {
     38  public:
     39    /// Type of the bipartite graph
     40    typedef BPG BpGraph;
     41    /// Type of the matching map
     42    typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap;
     43    /// Type of the elevator.
     44    typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t;
     45
     46  private:
     47    TEMPLATE_BPGRAPH_TYPEDEFS(BPG);
     48
     49  protected:
     50    typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow;
     51    typedef typename BPG::template BlueNodeMap<int> Excess;
     52    typedef std::vector<typename BPG::BlueNode> BlueNodes;
     53    typedef std::vector<typename BPG::Arc> BlueFlow;
     54    typedef typename BPG::template BlueNodeMap<BlueFlow> IncomingFlows;
     55
     56    const BpGraph& _bpg;
     57    MatchingMap* _matching;
     58    bool _local_matching;
     59    Elevator_t* _elevator;
     60    bool _local_elevator;
     61    Preflow _flow;
     62    Excess _excess;
     63
     64    void createStructures() {
     65      if (!_matching) {
     66        _matching = new MatchingMap(_bpg, INVALID);
     67        _local_matching = true;
     68      }
     69      if (!_elevator) {
     70        _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1,
     71                                                  countBlueNodes(_bpg)));
     72        _local_elevator = true;
     73      }
     74    }
     75
     76    void destroyStructures() {
     77      if (_local_matching) {
     78        delete _matching;
     79      }
     80      if (_local_elevator) {
     81        delete _elevator;
     82      }
     83    }
     84
     85    /// \brief Initialize the elevator
     86    ///
     87    /// When this function is called, the elevator is initialized with
     88    /// the exact distance labels corresponding the current matching.
     89    void initElevator() {
     90      BlueNodes act_blue_nodes, next_blue_nodes;
     91      BoolBlueNodeMap reached(_bpg, false);
     92      IncomingFlows incoming(_bpg);
     93
     94      // Set the initial flow:
     95      // choose arbitrary outgoing arc for every unmatched red node.
     96      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     97        Arc a = (*_matching)[r] != INVALID ?
     98          _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r);
     99        if (a != INVALID) {
     100          _flow[r] = a;
     101          BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(a));
     102          incoming[b].push_back(a);
     103          if (++_excess[b] == 2) {
     104            act_blue_nodes.push_back(b);
     105            reached.set(b, true);
     106          }
     107        }
     108      }
     109      // Now we start a Bfs from the blue nodes with excess
     110      // (those which get flow from more than one red node) in the
     111      // residual graph, to supply the exact distance labels to the elevator.
     112      _elevator->initStart();
     113      while (!act_blue_nodes.empty()) {
     114        for (typename BlueNodes::iterator b_it = act_blue_nodes.begin();
     115             b_it != act_blue_nodes.end(); ++b_it) {
     116
     117          _elevator->initAddItem(*b_it);
     118          for (typename BlueFlow::iterator b_in = incoming[*b_it].begin();
     119               b_in != incoming[*b_it].end(); ++b_in) {
     120            RedNode r = _bpg.asRedNodeUnsafe(_bpg.source(*b_in));
     121            for (OutArcIt r_oi(_bpg, r); r_oi!=INVALID; ++r_oi) {
     122              BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(r_oi));
     123              if (!reached[b]) {
     124                next_blue_nodes.push_back(b);
     125                reached.set(b, true);
     126              }
     127            }
     128          }
     129        }
     130        _elevator->initNewLevel();
     131        next_blue_nodes.swap(act_blue_nodes);
     132        next_blue_nodes.clear();
     133      }
     134      _elevator->initFinish();
     135
     136      // Activate blue nodes with 0 excess, i. e. those
     137      // which does not get any flow.
     138      for (BlueNodeIt b(_bpg); b != INVALID; ++b) {
     139        if (_excess[b] == 0  && (*_elevator)[b] < _elevator->maxLevel()) {
     140          _elevator->activate(b);
     141        }
     142      }
     143    }
     144
     145    PreflowBpMatching2() {}
     146
     147  public:
     148
     149    /// \brief Constructor
     150    ///
     151    /// Constructs the class on the given bipartite graph.
     152    PreflowBpMatching2(const BpGraph& bpg) : _bpg(bpg),
     153                                             _matching(0),
     154                                             _local_matching(false),
     155                                             _elevator(0),
     156                                             _local_elevator(false),
     157                                             _flow(bpg, INVALID),
     158                                             _excess(_bpg, 0)
     159
     160    {}
     161
     162    /// \brief Destructor
     163    ///
     164    /// Destructor.
     165    ~PreflowBpMatching2() {
     166      destroyStructures();
     167    }
     168
     169    /// \brief Sets the matching map
     170    ///
     171    /// Sets the matching map. It should contain a valid matching,
     172    /// for example the empty matching is fine.
     173    /// If you don't use this function before calling \ref init(),
     174    /// an instance will be allocated automatically.
     175    /// The destructor deallocates this automatically allocated map,
     176    /// of course.
     177    ///
     178    /// \return <tt>(*this)</tt>
     179    PreflowBpMatching2& matchingMap(MatchingMap& map) {
     180      _matching = &map;
     181      _local_matching = false;
     182      return *this;
     183    }
     184
     185    /// \brief Returns a const reference to the matching map.
     186    ///
     187    /// Returns a const reference to the matching map, which contains
     188    /// the matching edge for each node (or INVALID for unmatched nodes).
     189    const MatchingMap& matchingMap() const {
     190      return *_matching;
     191    }
     192
     193    /// \brief Sets the elevator used by algorithm.
     194    ///
     195    /// Sets the elevator.
     196    /// If you don't use this function before calling \ref init(),
     197    /// an instance will be allocated automatically.
     198    /// The destructor deallocates this automatically allocated elevator,
     199    /// of course.
     200    /// \return <tt>(*this)</tt>
     201    PreflowBpMatching2& elevator(Elevator_t& elevator) {
     202      _elevator = &elevator;
     203      _local_elevator = false;
     204      return *this;
     205    }
     206
     207    /// \brief Returns a const reference to the elevator.
     208    ///
     209    /// Returns a const reference to the elevator, which keeps track
     210    /// of the distance labels of nodes.
     211    const Elevator_t& elevator() const {
     212      return *_elevator;
     213    }
     214
     215    /// \brief Initializes the algorithm
     216    ///
     217    /// Allocates the elevator and the matching map if necessary,
     218    /// and initializes the elevator.
     219    ///
     220    /// \pre \ref init() is not called.
     221    void init() {
     222      createStructures();
     223      initElevator();
     224    }
     225
     226    /// \brief Smarter initialization of the matching.
     227    ///
     228    /// Allocates the matching map if necessary, and initializes
     229    /// the algorithm with a trivially found matching: iterating
     230    /// on every edge, take it in if both endnodes are unmatched.
     231    ///
     232    /// \return Size of the so found matching.
     233    ///
     234    /// \note heuristicInit() replaces init() regarding the preconditions.
     235    int heuristicInit() {
     236      createStructures();
     237      int size = 0;
     238
     239      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     240        bool matched = false;
     241        for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) {
     242          if ((*_matching)[_bpg.target(oi)] == INVALID) {
     243            _matching->set(_bpg.source(oi), oi);
     244            _matching->set(_bpg.target(oi), oi);
     245            matched = true;
     246            ++size;
     247          }
     248        }
     249      }
     250
     251      initElevator();
     252      return size;
     253    }
     254
     255    /// \brief Initialize the matching from a map.
     256    ///
     257    /// Allocates the matching map if necessary, and
     258    /// initializes the algorithm with a matching given as a bool edgemap,
     259    /// in which an edge is true if it is in the matching.
     260    /// Also initializes the elevator.
     261    ///
     262    /// \return \c true if the matching is valid.
     263    ///
     264    /// \note matchingInit() replaces init() regarding the preconditions.
     265    bool matchingInit(const BoolEdgeMap& matching) {
     266      createStructures();
     267
     268      bool valid = true;
     269      for(EdgeIt it(_bpg); it!=INVALID; ++it) {
     270        if (matching[it]) {
     271          Node red = _bpg.redNode(it);
     272          Node blue = _bpg.blueNode(it);
     273          if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) {
     274            valid = false;
     275          } else {
     276            _matching->set(red, it);
     277            _matching->set(blue, it);
     278          }
     279        }
     280      }
     281
     282      initElevator();
     283      return valid;
     284    }
     285
     286
     287    /// \brief Perform a pull/relabel operation on a node
     288    ///
     289    /// Perform a pull/relabel operation on a node:
     290    /// search an incident blue node on one level lower,
     291    /// and flip its flow to the 'chosen', or relabel
     292    /// if none found.
     293    ///
     294    /// \pre \ref init() is called, and the node is active.
     295    void pull_relabel(const BlueNode& chosen) {
     296
     297      int min_inc_level = _elevator->maxLevel() + 1;
     298      for (OutArcIt b_oi(_bpg, chosen); b_oi!=INVALID; ++b_oi) {
     299        RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(b_oi));
     300        BlueNode lower = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r]));
     301        if ((*_elevator)[(lower)] < (*_elevator)[chosen]) {
     302          if (--_excess[lower] == 0) {
     303            _elevator->activate(lower);
     304          }
     305          _flow[r] = _bpg.direct(b_oi, true);
     306          ++_excess[chosen];
     307          _elevator->deactivate(chosen);
     308          break;
     309        } else {
     310          if ((*_elevator)[lower] < min_inc_level &&
     311              lower != chosen) {
     312            min_inc_level = (*_elevator)[lower];
     313          }
     314        }
     315      }
     316      // If none found, then relabel.
     317      if (_excess[chosen] == 0) {
     318        int chosen_level = (*_elevator)[chosen];
     319        if (min_inc_level >= _elevator->maxLevel() - 1) {
     320          _elevator->lift(chosen, _elevator->maxLevel());
     321          _elevator->deactivate(chosen);
     322        } else {
     323          _elevator->lift(chosen, min_inc_level+1);
     324        }
     325        if (_elevator->emptyLevel(chosen_level)) {
     326          _elevator->liftToTop(chosen_level);
     327        }
     328      }
     329    }
     330
     331    /// \brief Executes the algorithm
     332    ///
     333    /// Executes pull/relabel operation on the highest active node
     334    /// until there is. After start(), writeMatching() should be called to
     335    /// remove unnecessary flow and choose a valid maximum matching.
     336    ///
     337    /// \pre \ref init() must be called before using this function.
     338    void start() {
     339      while (_elevator->highestActive() != INVALID)
     340        pull_relabel(_elevator->highestActive());
     341    }
     342
     343    /// \brief Deduces and writes a matching from the preflow
     344    ///
     345    /// The underlying preflow algorithm does not use the matching map,
     346    /// so before using the matching query functions ( \ref mate(), \ref
     347    /// matchingSize(), etc.), you should explicitly call this function
     348    /// to write a matching based on the current preflow to the matching
     349    /// map. When the preflow results an ambiguous matching (a node has
     350    /// incoming flow on more than one arc), it chooses arbitrary one.
     351    void writeMatching() {
     352      for (NodeIt it(_bpg); it != INVALID; ++it) {
     353        _matching->set(it, INVALID);
     354      }
     355      for (RedNodeIt r(_bpg); r != INVALID; ++r) {
     356        if (_flow[r] != INVALID) {
     357          BlueNode mate = _bpg.blueNode(_flow[r]);
     358          if (_excess[mate] > 1) {
     359            --_excess[mate];
     360          } else {
     361            _matching->set(mate, _flow[r]);
     362            _matching->set(r, _flow[r]);
     363          }
     364        }
     365      }
     366    }
     367
     368    /// \brief Runs the algorithm.
     369    ///
     370    /// bp1.run() is just a shorthand for:
     371    ///
     372    ///\code
     373    /// bp1.heuristicInit();
     374    /// bp1.start();
     375    /// bp1.writeMatching();
     376    ///\endcode
     377    void run() {
     378      heuristicInit();
     379      start();
     380      writeMatching();
     381    }
     382
     383    /// \brief Size of the matching
     384    ///
     385    /// Returns the size of the current matching.
     386    /// Function runs in \f$O(n)\f$ time.
     387    int matchingSize() const {
     388      int size = 0;
     389      for (RedNodeIt it(_bpg); it!=INVALID; ++it) {
     390        if ((*_matching)[it] != INVALID) ++size;
     391      }
     392      return size;
     393    }
     394
     395    /// \brief Return \c true if the given edge is in the matching.
     396    ///
     397    /// This function returns \c true if the given edge is in the current
     398    /// matching.
     399    bool matching(const Edge& edge) const {
     400      return edge == (*_matching)[_bpg.redNode(edge)];
     401    }
     402
     403    /// \brief Return the matching edge incident to the given node.
     404    ///
     405    /// This function returns the matching edge incident to the
     406    /// given node in the current matching or \c INVALID if the node is
     407    /// not covered by the matching.
     408    Edge matching(const Node& n) const {
     409      return (*_matching)[n];
     410    }
     411
     412    /// \brief Return the mate of the given node.
     413    ///
     414    /// This function returns the mate of the given node in the current
     415    /// matching or \c INVALID if the node is not covered by the matching.
     416    Node mate(const Node& n) const {
     417      return (*_matching)[n] != INVALID ?
     418        _bpg.oppositeNode(n, (*_matching)[n]) : INVALID;
     419    }
     420
     421    /// \brief Query a vertex cover.
     422    ///
     423    /// This function finds a vertex cover based on the current matching,
     424    /// and in \c cover sets node in the cover \c true and \c false
     425    /// otherwise. In case of maximal matching, this will be a minimal
     426    /// vertex cover.
     427    /// Function runs in \f$O(n + e)\f$ time.
     428    ///
     429    /// \return The number of nodes in the cover.
     430    int cover(BoolNodeMap &cover) const {
     431      int count = 0;
     432      std::queue<RedNode> q;
     433
     434      for (BlueNodeIt b(_bpg); b!=INVALID; ++b)
     435        cover[b] = false;
     436
     437      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     438        if ((*_matching)[r] == INVALID) {
     439          cover[r] = false;
     440          q.push(r);
     441        } else {
     442          cover[r] = true;
     443          ++count;
     444        }
     445      }
     446
     447      while (!q.empty()) {
     448        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     449          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     450          if (cover[blue]) continue;
     451
     452          cover[blue] = true;
     453          ++count;
     454
     455          if ((*_matching)[blue] != INVALID) {
     456            RedNode nr = _bpg.redNode((*_matching)[blue]);
     457            q.push(nr);
     458            cover[nr] = false;
     459            --count;
     460          }
     461        }
     462        q.pop();
     463      }
     464
     465      return count;
     466    }
     467
     468    /// \brief Query a barrier of red nodes.
     469    ///
     470    /// This function tries to create a barrier of red nodes, which should:
     471    ///  - contain every unmatched red nodes,
     472    ///  - contain exactly that many matched red nodes as much blue nodes
     473    ///    there are in its neighbor set.
     474    ///
     475    /// Barrier exist if and only if the matching is maximal. If exist, the
     476    /// map is set true for nodes of the barrier, and false for the others.
     477    ///
     478    /// Function runs in \f$O(n + e)\f$ time.
     479    ///
     480    /// \return \c true if a barrier found.
     481    bool redBarrier(BoolRedNodeMap &barrier) const {
     482      bool exist = true;
     483      std::queue<RedNode> q;
     484
     485      for (RedNodeIt r(_bpg); r!=INVALID; ++r) {
     486        if ((*_matching)[r] == INVALID) {
     487          barrier[r] = true;
     488          q.push(r);
     489        } else {
     490          barrier[r] = false;
     491        }
     492      }
     493
     494      while (!q.empty() && exist) {
     495        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     496          BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it)));
     497
     498          if (exist &= ((*_matching)[blue] != INVALID)) {
     499            RedNode r = _bpg.redNode((*_matching)[blue]);
     500            if (!barrier[r]) {
     501              q.push(r);
     502              barrier[r] = true;
     503            }
     504          }
     505        }
     506        q.pop();
     507      }
     508
     509      return exist;
     510    }
     511
     512    /// \brief Query a barrier of blue nodes.
     513    ///
     514    /// This function tries to create a barrier of blue nodes, which should:
     515    ///  - contain every unmatched blue nodes,
     516    ///  - contain exactly that many matched blue nodes as much red nodes
     517    ///    there are in its neighbor set.
     518    ///
     519    /// Barrier exist if and only if the matching is maximal. If exist, the
     520    /// map is set true for nodes of the barrier, and false for the others.
     521    ///
     522    /// Function runs in \f$O(n + e)\f$ time.
     523    ///
     524    /// \return \c true if a barrier found.
     525    bool blueBarrier(BoolBlueNodeMap &barrier) const {
     526      bool exist = true;
     527      std::queue<BlueNode> q;
     528
     529      for (BlueNodeIt b(_bpg); b!=INVALID; ++b) {
     530        if ((*_matching)[b] == INVALID) {
     531          barrier[b] = true;
     532          q.push(b);
     533        } else {
     534          barrier[b] = false;
     535        }
     536      }
     537
     538      while (!q.empty() && exist) {
     539        for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) {
     540          RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it)));
     541
     542          if (exist &= ((*_matching)[r] != INVALID)) {
     543            BlueNode b = _bpg.blueNode((*_matching)[r]);
     544            if (!barrier[b]) {
     545              q.push(b);
     546              barrier[b] = true;
     547            }
     548          }
     549        }
     550        q.pop();
     551      }
     552
     553      return exist;
     554    }
     555  };
     556
     557}
     558#endif
  • test/CMakeLists.txt

    diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
    a b  
    4848  path_test
    4949  planarity_test
    5050  preflow_test
     51  preflow_bp_matching_2_test
    5152  preflow_bp_matching_test
    5253  radix_sort_test
    5354  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
    - +  
     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-2012
     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/preflow_bp_matching_2.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
     71const string u_lgf =
     72  "@red_nodes\n"
     73  "label\n"
     74  "1\n"
     75  "2\n"
     76  "3\n"
     77  "@blue_nodes\n"
     78  "label\n"
     79  "4\n"
     80  "5\n"
     81  "6\n"
     82  "7\n"
     83  "@edges\n"
     84  "label\n"
     85  "1 4 1\n"
     86  "1 5 2\n"
     87  "2 5 3\n"
     88  "2 6 4\n"
     89  "2 7 5\n"
     90  "3 7 6\n"
     91;
     92
     93
     94
     95void checkPreflowBpMatching2Compile() {
     96  typedef concepts::BpGraph BpGraph;
     97  BPGRAPH_TYPEDEFS(BpGraph);
     98  typedef PreflowBpMatching2<BpGraph> PBM;
     99
     100  BpGraph graph;
     101  RedNode red;
     102  BlueNode blue;
     103  Node node;
     104  Edge edge;
     105  BoolEdgeMap init_matching(graph);
     106  PBM pbm_test(graph);
     107  const PBM& const_pbm_test = pbm_test;
     108  PBM::MatchingMap matching_map(graph, INVALID);
     109  PBM::Elevator_t elevator(graph);
     110  BpGraph::NodeMap<bool> cov(graph);
     111  int cover_size;
     112  BoolRedNodeMap rb(graph);
     113  BoolBlueNodeMap bb(graph);
     114
     115  pbm_test.matchingMap(matching_map);
     116  pbm_test.elevator(elevator);
     117  pbm_test.init();
     118  pbm_test.heuristicInit();
     119  pbm_test.matchingInit(init_matching);
     120
     121  pbm_test.start();
     122  pbm_test.pull_relabel(blue);
     123  pbm_test.run();
     124  pbm_test.writeMatching();
     125
     126  const_pbm_test.matchingSize();
     127  const_pbm_test.matching(node);
     128  const_pbm_test.matching(edge);
     129  const_pbm_test.mate(red);
     130  const_pbm_test.mate(blue);
     131  const_pbm_test.mate(node);
     132  const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap();
     133  edge = max_matching[node];
     134  const PBM::Elevator_t& elev = const_pbm_test.elevator();
     135  blue = elev.highestActive();
     136  cover_size = const_pbm_test.cover(cov);
     137  const_pbm_test.redBarrier(rb);
     138  const_pbm_test.blueBarrier(bb);
     139}
     140
     141typedef SmartBpGraph BpGraph;
     142typedef PreflowBpMatching2<BpGraph> PBM;
     143BPGRAPH_TYPEDEFS(BpGraph);
     144
     145void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) {
     146  BoolBlueNodeMap matched(bpg, false);
     147  for (RedNodeIt it(bpg); it != INVALID; ++it) {
     148    if (pbm.matching(it) != INVALID) {
     149      check(pbm.mate(pbm.mate(it)) == it, "Wrong matching");
     150      matched.set(bpg.asBlueNode(pbm.mate(it)), true);
     151    }
     152  }
     153  for (BlueNodeIt it(bpg); it != INVALID; ++it) {
     154    // != is used as logical xor here
     155    check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching");
     156  }
     157}
     158
     159void checkCover(const BpGraph& bpg, const PBM& pbm) {
     160  int cov_size;
     161  BoolNodeMap cov(bpg);
     162  cov_size = pbm.cover(cov);
     163  BoolEdgeMap covered(bpg, false);
     164  for (NodeIt n(bpg); n!=INVALID; ++n) {
     165    for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) {
     166      covered[e] = true;
     167    }
     168  }
     169  bool all = true;
     170  for (EdgeIt e(bpg); e!=INVALID; ++e) {
     171    all&=covered[e];
     172  }
     173
     174  check(all, "Invalid cover");
     175  check(pbm.matchingSize() == cov_size, "Suboptimal matching.");
     176}
     177void checkBarriers(const BpGraph& bpg, const PBM& pbm) {
     178  BoolRedNodeMap rb(bpg);
     179  BoolBlueNodeMap bb(bpg);
     180
     181  check(pbm.redBarrier(rb), "Suboptimal matching.");
     182  check(pbm.blueBarrier(bb), "Suboptimal matching.");
     183
     184  BoolNodeMap reached(bpg, false);
     185
     186  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     187    if (!rb[r]) continue;
     188    for (OutArcIt a(bpg, r); a!=INVALID; ++a) {
     189      reached.set(bpg.target(a), true);
     190    }
     191  }
     192  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     193    if (!bb[b]) continue;
     194    for (OutArcIt a(bpg, b); a!=INVALID; ++a) {
     195      reached.set(bpg.target(a), true);
     196    }
     197  }
     198
     199  int ur = 0, ub = 0, rc = 0, bc = 0;
     200  for (RedNodeIt r(bpg); r!=INVALID; ++r) {
     201    if (rb[r]) ++rc;
     202    if (pbm.mate(r) == INVALID) ++ur;
     203    if (reached[r]) --bc;
     204  }
     205
     206  for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     207    if (bb[b]) ++bc;
     208    if (pbm.mate(b) == INVALID) ++ub;
     209    if (reached[b]) --rc;
     210  }
     211
     212  check(rc == ur, "Red barrier is wrong.");
     213  check(bc == ub, "Blue barrier is wrong.");
     214}
     215
     216void checkMatching(const BpGraph& bpg, const PBM& pbm) {
     217  checkMatchingValidity(bpg, pbm);
     218  checkBarriers(bpg, pbm);
     219  checkCover(bpg, pbm);
     220}
     221
     222
     223int main() {
     224  // Checking hard wired extremal graphs
     225  for (int i=0; i<lgfn; ++i) {
     226    BpGraph bpg;
     227    istringstream lgfs(lgf[i]);
     228    bpGraphReader(bpg, lgfs).run();
     229    PBM pbm(bpg);
     230    pbm.run();
     231    checkMatching(bpg, pbm);
     232  }
     233
     234  // Checking some use cases
     235  BpGraph bpg;
     236  istringstream u_lgfs(u_lgf);
     237  bpGraphReader(bpg, u_lgfs).run();
     238
     239  {
     240    PBM pbm(bpg);
     241    pbm.init();
     242    BpGraph::BlueNode b = pbm.elevator().highestActive();
     243    if (b != INVALID) {
     244      pbm.pull_relabel(b);
     245    }
     246    pbm.start();
     247    pbm.writeMatching();
     248    check(pbm.elevator().highestActive() == INVALID, "Error in elevator");
     249    checkMatching(bpg, pbm);
     250  }
     251
     252  {
     253    PBM pbm(bpg);
     254    pbm.heuristicInit();
     255    BpGraph::BlueNode b = pbm.elevator().highestActive();
     256    if (b != INVALID) {
     257      pbm.pull_relabel(b);
     258    }
     259    pbm.start();
     260    checkMatching(bpg, pbm);
     261  }
     262  {
     263    PBM pbm(bpg);
     264    PBM::MatchingMap m(bpg, INVALID);
     265    pbm.matchingMap(m);
     266    pbm.run();
     267    checkMatching(bpg, pbm);
     268
     269  }
     270  {
     271    PBM pbm(bpg);
     272    BpGraph::EdgeMap<bool> init_matching(bpg, false);
     273    BpGraph::Edge e;
     274    bpg.first(e);
     275    init_matching[e] = true;
     276    check(pbm.matchingInit(init_matching),
     277          "Wrong result from matchingInit()");
     278    pbm.start();
     279    pbm.writeMatching();
     280    checkMatching(bpg, pbm);
     281  }
     282  {
     283    PBM pbm(bpg);
     284    PBM::MatchingMap m(bpg, INVALID);
     285    pbm.matchingMap(m);
     286    BpGraph::EdgeMap<bool> init_matching(bpg, true);
     287    check(!pbm.matchingInit(init_matching),
     288          "Wrong result from matchingInit()");
     289    pbm.start();
     290    pbm.writeMatching();
     291    checkMatching(bpg, pbm);
     292  }
     293
     294  // Checking some random generated graphs
     295  const int random_test_n = 10;
     296  const int max_red_n = 1000;
     297  const int min_red_n = 10;
     298  const int max_blue_n = 1000;
     299  const int min_blue_n = 10;
     300  for (int i=0; i<random_test_n; ++i) {
     301    int red_n = rnd.integer(min_red_n, max_red_n);
     302    int blue_n = rnd.integer(min_blue_n, max_blue_n);
     303    double prob = rnd();
     304
     305    SmartBpGraph bpg;
     306    for (int j=0; j<red_n; ++j) {
     307      bpg.addRedNode();
     308    }
     309    for (int j=0; j<blue_n; ++j) {
     310      bpg.addBlueNode();
     311    }
     312    for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) {
     313      for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) {
     314        if (rnd() < prob/10.) {
     315          bpg.addEdge(r,b);
     316        }
     317      }
     318    }
     319
     320    PBM pbm(bpg);
     321    pbm.run();
     322    checkMatching(bpg, pbm);
     323
     324    // do not always use the highest active node
     325    PBM pbm2(bpg);
     326    PBM::Elevator_t elev(bpg);
     327    pbm2.elevator(elev);
     328    pbm2.init();
     329    while (elev.highestActive() != INVALID) {
     330      for (BlueNodeIt b(bpg); b!=INVALID; ++b) {
     331        if (elev.active(b)) pbm2.pull_relabel(b);
     332      }
     333    }
     334    pbm2.writeMatching();
     335    checkMatching(bpg, pbm2);
     336  }
     337
     338
     339  return 0;
     340}
     341
     342
  • 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 b  
    2323#include <lemon/hopcroft_karp.h>
    2424#include <lemon/matching.h>
    2525#include <lemon/preflow.h>
     26#include <lemon/preflow_bp_matching.h>
     27#include <lemon/preflow_bp_matching_2.h>
    2628#include <lemon/bp_matching.h>
    2729#include <lemon/cost_scaling.h>
    2830#include <lemon/capacity_scaling.h>
     
    3234using namespace std;
    3335using namespace lemon;
    3436
     37void check(bool c, const string &msg) {
     38  if (!c) {
     39    std::cerr << "Error: " << msg << "\n";
     40  }
     41}
    3542
    3643template <typename BpGraph, typename Digraph>
    3744void create_network(const BpGraph &bpg,
     
    108115///
    109116/// Depending on the input, it runs the \ref lemon::HopcroftKarp
    110117/// "Hopcroft-Karp", the \ref lemon::MaxMatching "general matching"
    111 /// and the \ref lemon::Preflow "Preflow" algorithms for finding
     118/// and three algorithms based on preflow (\ref lemon::Preflow "preflow
     119/// on extended graph", \ref lemon::PreflowBpMatching,
     120/// \ref lemon::PreflowBpMatching2") for finding
    112121/// a maximum cardinality bipartite matching, or the
    113122/// \ref lemon::MaxWeightedBpMatching "max. weighted matching
    114123/// for sparse bipartite graphs", the \ref lemon::MaxWeightedMatching
     
    118127/// simplex" algorithms to find the maximum weighted matching.
    119128int main() {
    120129  DimacsDescriptor desc = dimacsType(cin);
    121   int red_n, blue_n, edge_n = desc.edgeNum;
    122130  if (desc.type == DimacsDescriptor::UBM) {
    123131    // unweighted bipartite matching
    124     Timer hk_t(false), mm_t(false), pf_t(false);
    125     int hk_r, mm_r, pf_r;
     132    Timer hk_t(false),
     133      mm_t(false),
     134      pf_t(false),
     135      pbm_t(false),
     136      pbm2_t(false);
    126137
    127138    BpGraph *bpg = new BpGraph();
    128139    readDimacsUbm(cin, *bpg, desc);
    129     red_n = countRedNodes(*bpg);
    130     blue_n = countBlueNodes(*bpg);
    131140
    132141    {
    133142      hk_t.start();
    134143      HopcroftKarp<BpGraph> hk(*bpg);
    135144      hk.run();
    136       hk_r = hk.matchingSize();
    137145      hk_t.stop();
     146      BpGraph::NodeMap<bool> m(*bpg);
     147      check(hk.matchingSize() == hk.cover(m),
     148            "Suboptimal matching in Hopcroft-Karp.");
    138149    }
    139150    {
    140151      mm_t.start();
    141152      MaxMatching<BpGraph> mm(*bpg);
    142153      mm.run();
    143       mm_r = mm.matchingSize();
    144154      mm_t.stop();
    145155    }
     156    {
     157      pbm_t.start();
     158      PreflowBpMatching<BpGraph> pbm(*bpg);
     159      pbm.run();
     160      pbm_t.stop();
     161      BpGraph::NodeMap<bool> m(*bpg);
     162      check(pbm.matchingSize() == pbm.cover(m),
     163            "Suboptimal matching in push-relabel.");
     164    }
     165    {
     166      pbm2_t.start();
     167      PreflowBpMatching2<BpGraph> pbm2(*bpg);
     168      pbm2.run();
     169      pbm2_t.stop();
     170      BpGraph::NodeMap<bool> m(*bpg);
     171      check(pbm2.matchingSize() == pbm2.cover(m),
     172            "Suboptimal matching in pull-relabel.");
     173    }
    146174
    147175    Digraph *network = new Digraph();
    148176    Digraph::Node source, target;
    149177    create_network(*bpg, *network, source, target);
    150     delete bpg;
    151178    {
    152179      pf_t.start();
    153180      Preflow<Digraph, ConstMap<Digraph::Arc, long> >
    154181        pf(*network, constMap<Digraph::Arc, long>(1), source, target);
    155182      pf.run();
    156       pf_r = pf.flowValue();
    157183      pf_t.stop();
    158184    }
    159185    delete network;
    160186
    161187    cout << "Benchmarking in unweighted case, using a bipartite graph\n"
    162          << "with " << red_n << " red, " << blue_n << " blue nodes, and "
    163          << edge_n << " edges.\n"
     188         << "with " << countRedNodes(*bpg) << " red, "
     189         << countBlueNodes(*bpg) << " blue nodes, and "
     190         << desc.edgeNum << " edges.\n"
    164191         << "--------------------------------------------------------\n"
    165          << "Algorithm used       Matching size       Time\n"
    166          << "Hopcroft-Karp            " << hk_r
    167          << "         " << hk_t.realTime() << "\n"
    168          << "General matching         " << mm_r
    169          << "         " << mm_t.realTime() << "\n"
    170          << "Preflow                  " << pf_r
    171          << "         " << pf_t.realTime() << endl;
     192         << "Algorithm used               Time\n"
     193         << "Hopcroft-Karp            " << hk_t.realTime() << "\n"
     194         << "General matching         " << mm_t.realTime() << "\n"
     195         << "Push-relabel matching    " << pbm_t.realTime() << "\n"
     196         << "Pull-relabel matching    " << pbm2_t.realTime() << "\n"
     197         << "Preflow                  " << pf_t.realTime() << endl;
    172198
     199    delete bpg;
    173200  } else if (desc.type == DimacsDescriptor::WBM) {
    174201    // weighted bipartite matching
    175202    Timer bm_t(false),
     
    182209    BpGraph *bpg = new BpGraph();
    183210    BpGraph::EdgeMap<long> *weight = new BpGraph::EdgeMap<long>(*bpg);
    184211    readDimacsWbm(cin, *bpg, *weight, desc);
    185     red_n = countRedNodes(*bpg);
    186     blue_n = countBlueNodes(*bpg);
    187212
    188213    {
    189214      bm_t.start();
     
    241266    delete supply;
    242267
    243268    cout << "Benchmarking on weighted bipartite matching problem on a "
    244          << "graph\nwith " << red_n << " red, " << blue_n
    245          << " blue nodes and " << edge_n << " edges\n"
     269         << "graph\nwith " << countRedNodes(*bpg) << " red, "
     270         << countBlueNodes(*bpg) << " blue nodes and "
     271         << desc.edgeNum << " edges\n"
    246272         << "--------------------------------------------------------\n"
    247273         << "Algorithm used       Maximum weight       Time\n"
    248274         << "Bipartite matching       " << bm_r