COIN-OR::LEMON - Graph Library

Ticket #168: 5b4317f3ce36.patch

File 5b4317f3ce36.patch, 48.0 KB (added by Balazs Dezso, 2 years ago)
  • new file lemon/bpmatching.h

    # HG changeset patch
    # User Balazs Dezso <deba@google.com>
    # Date 1548509785 -3600
    #      Sat Jan 26 14:36:25 2019 +0100
    # Node ID 5b4317f3ce36acee1eb1ac9c20f803eb36628847
    # Parent  8c567e298d7f3fad82cc66754f2fb6c4a420366b
    Maximum weighted matching and maximum weighted perfect matching for bipartite graphs
    
    diff -r 8c567e298d7f -r 5b4317f3ce36 lemon/bpmatching.h
    - +  
     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-2019
     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_BPMATCHING_H
     20#define LEMON_BPMATCHING_H
     21
     22#include <limits>
     23
     24#include <lemon/core.h>
     25#include <lemon/unionfind.h>
     26#include <lemon/bin_heap.h>
     27#include <lemon/maps.h>
     28
     29///\ingroup matching
     30///\file
     31///\brief Maximum matching algorithms in bipartite graphs.
     32
     33namespace lemon {
     34
     35  /// \ingroup matching
     36  ///
     37  /// \brief Weighted matching in bipartite graphs
     38  ///
     39  /// This class provides an efficient implementation of multiple search tree
     40  /// augmenting path matching algorithm. The implementation is based on
     41  /// extensive use of priority queues and provides \f$O(nm\log n)\f$ time
     42  /// complexity.
     43  ///
     44  /// The maximum weighted matching problem is to find a subset of the
     45  /// edges in a bipartite graph with maximum overall weight for which
     46  /// each node has at most one incident edge.
     47  /// It can be formulated with the following linear program.
     48  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
     49  /// \f[x_e \ge 0\quad \forall e\in E\f]
     50  /// \f[\max \sum_{e\in E}x_ew_e\f]
     51  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
     52  /// \f$X\f$.
     53  ///
     54  /// The algorithm calculates an optimal matching and a proof of the
     55  /// optimality. The solution of the dual problem can be used to check
     56  /// the result of the algorithm. The dual linear problem is the
     57  /// following.
     58  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
     59  /// \f[y_u \ge 0 \quad \forall u \in V\f]
     60  /// \f[\min \sum_{u \in V}y_u \f]
     61  /// \tparam BPGR The bipartite graph type the algorithm runs on.
     62  /// \tparam WM The type edge weight map. The default type is
     63  /// \ref concepts::BpGraph::EdgeMap "BPGR::EdgeMap<int>".
     64#ifdef DOXYGEN
     65  template <typename BPGR, typename WM>
     66#else
     67  template <typename BPGR,
     68            typename WM = typename BPGR::template EdgeMap<int> >
     69#endif
     70  class MaxWeightedBpMatching {
     71  public:
     72
     73    /// The graph type of the algorithm
     74    typedef BPGR BpGraph;
     75    /// The type of the edge weight map
     76    typedef WM WeightMap;
     77    /// The value type of the edge weights
     78    typedef typename WeightMap::Value Value;
     79
     80    /// The type of the matching map
     81    typedef typename BpGraph::template NodeMap<typename BpGraph::Arc>
     82    MatchingMap;
     83
     84    /// \brief Scaling factor for dual solution
     85    ///
     86    /// Scaling factor for dual solution. It is equal to 4 or 1
     87    /// according to the value type.
     88    static const int dualScale =
     89      std::numeric_limits<Value>::is_integer ? 4 : 1;
     90
     91  private:
     92
     93    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     94
     95    typedef typename BpGraph::template NodeMap<Value> NodePotential;
     96
     97    const BpGraph& _bpgraph;
     98    const WeightMap& _weight;
     99
     100    MatchingMap* _matching;
     101
     102    NodePotential* _node_potential;
     103
     104    int _node_num;
     105
     106    enum Status {
     107      EVEN = -1, MATCHED = 0, ODD = 1
     108    };
     109
     110    typedef typename BpGraph::template NodeMap<Status> StatusMap;
     111    StatusMap* _status;
     112
     113    typedef typename BpGraph::template NodeMap<Arc> PredMap;
     114    PredMap* _pred;
     115
     116    typedef ExtendFindEnum<IntNodeMap> TreeSet;
     117    IntNodeMap *_tree_set_index;
     118    TreeSet *_tree_set;
     119
     120    IntNodeMap *_delta1_index;
     121    BinHeap<Value, IntNodeMap> *_delta1;
     122
     123    IntNodeMap *_delta2_index;
     124    BinHeap<Value, IntNodeMap> *_delta2;
     125
     126    IntEdgeMap *_delta3_index;
     127    BinHeap<Value, IntEdgeMap> *_delta3;
     128
     129    Value _delta_sum;
     130    int _unmatched;
     131
     132    void createStructures() {
     133      _node_num = countNodes(_bpgraph);
     134
     135      if (!_matching) {
     136        _matching = new MatchingMap(_bpgraph);
     137      }
     138
     139      if (!_node_potential) {
     140        _node_potential = new NodePotential(_bpgraph);
     141      }
     142
     143      if (!_status) {
     144        _status = new StatusMap(_bpgraph);
     145      }
     146
     147      if (!_pred) {
     148        _pred = new PredMap(_bpgraph);
     149      }
     150
     151      if (!_tree_set) {
     152        _tree_set_index = new IntNodeMap(_bpgraph);
     153        _tree_set = new TreeSet(*_tree_set_index);
     154      }
     155
     156      if (!_delta1) {
     157        _delta1_index = new IntNodeMap(_bpgraph);
     158        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
     159      }
     160
     161      if (!_delta2) {
     162        _delta2_index = new IntNodeMap(_bpgraph);
     163        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
     164      }
     165
     166      if (!_delta3) {
     167        _delta3_index = new IntEdgeMap(_bpgraph);
     168        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
     169      }
     170    }
     171
     172    void destroyStructures() {
     173      if (_matching) {
     174        delete _matching;
     175      }
     176      if (_node_potential) {
     177        delete _node_potential;
     178      }
     179      if (_status) {
     180        delete _status;
     181      }
     182      if (_pred) {
     183        delete _pred;
     184      }
     185      if (_tree_set) {
     186        delete _tree_set_index;
     187        delete _tree_set;
     188      }
     189      if (_delta2) {
     190        delete _delta2_index;
     191        delete _delta2;
     192      }
     193      if (_delta3) {
     194        delete _delta3_index;
     195        delete _delta3;
     196      }
     197    }
     198
     199    void matchedToEven(Node node, int tree) {
     200      _tree_set->insert(node, tree);
     201      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     202      _delta1->push(node, (*_node_potential)[node]);
     203
     204      if (_delta2->state(node) == _delta2->IN_HEAP) {
     205        _delta2->erase(node);
     206      }
     207
     208      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     209        Node v = _bpgraph.source(a);
     210        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     211          dualScale * _weight[a];
     212        if ((*_status)[v] == EVEN) {
     213          _delta3->push(a, rw / 2);
     214        } else if ((*_status)[v] == MATCHED) {
     215          if (_delta2->state(v) != _delta2->IN_HEAP) {
     216            _pred->set(v, a);
     217            _delta2->push(v, rw);
     218          } else if ((*_delta2)[v] > rw) {
     219            _pred->set(v, a);
     220            _delta2->decrease(v, rw);
     221          }
     222        }
     223      }
     224    }
     225
     226    void matchedToOdd(Node node, int tree) {
     227      _tree_set->insert(node, tree);
     228      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     229
     230      if (_delta2->state(node) == _delta2->IN_HEAP) {
     231        _delta2->erase(node);
     232      }
     233    }
     234
     235    void evenToMatched(Node node, int tree) {
     236      _delta1->erase(node);
     237      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     238      Arc min = INVALID;
     239      Value minrw = std::numeric_limits<Value>::max();
     240      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     241        Node v = _bpgraph.source(a);
     242        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     243          dualScale * _weight[a];
     244
     245        if ((*_status)[v] == EVEN) {
     246          _delta3->erase(a);
     247          if (minrw > rw) {
     248            min = _bpgraph.oppositeArc(a);
     249            minrw = rw;
     250          }
     251        } else if ((*_status)[v]  == MATCHED) {
     252          if ((*_pred)[v] == a) {
     253            Arc mina = INVALID;
     254            Value minrwa = std::numeric_limits<Value>::max();
     255            for (OutArcIt aa(_bpgraph, v); aa != INVALID; ++aa) {
     256              Node va = _bpgraph.target(aa);
     257              if ((*_status)[va] != EVEN ||
     258                  _tree_set->find(va) == tree) continue;
     259              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
     260                dualScale * _weight[aa];
     261              if (minrwa > rwa) {
     262                minrwa = rwa;
     263                mina = aa;
     264              }
     265            }
     266            if (mina != INVALID) {
     267              _pred->set(v, mina);
     268              _delta2->increase(v, minrwa);
     269            } else {
     270              _pred->set(v, INVALID);
     271              _delta2->erase(v);
     272            }
     273          }
     274        }
     275      }
     276      if (min != INVALID) {
     277        _pred->set(node, min);
     278        _delta2->push(node, minrw);
     279      } else {
     280        _pred->set(node, INVALID);
     281      }
     282    }
     283
     284    void oddToMatched(Node node) {
     285      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     286      Arc min = INVALID;
     287      Value minrw = std::numeric_limits<Value>::max();
     288      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     289        Node v = _bpgraph.source(a);
     290        if ((*_status)[v] != EVEN) continue;
     291        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     292          dualScale * _weight[a];
     293
     294        if (minrw > rw) {
     295          min = _bpgraph.oppositeArc(a);
     296          minrw = rw;
     297        }
     298      }
     299      if (min != INVALID) {
     300        _pred->set(node, min);
     301        _delta2->push(node, minrw);
     302      } else {
     303        _pred->set(node, INVALID);
     304      }
     305    }
     306
     307    void alternatePath(Node even, int tree) {
     308      Node odd;
     309
     310      _status->set(even, MATCHED);
     311      evenToMatched(even, tree);
     312
     313      Arc prev = (*_matching)[even];
     314      while (prev != INVALID) {
     315        odd = _bpgraph.target(prev);
     316        even = _bpgraph.target((*_pred)[odd]);
     317        _matching->set(odd, (*_pred)[odd]);
     318        _status->set(odd, MATCHED);
     319        oddToMatched(odd);
     320
     321        prev = (*_matching)[even];
     322        _status->set(even, MATCHED);
     323        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
     324        evenToMatched(even, tree);
     325      }
     326    }
     327
     328    void destroyTree(int tree) {
     329      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
     330        if ((*_status)[n] == EVEN) {
     331          _status->set(n, MATCHED);
     332          evenToMatched(n, tree);
     333        } else if ((*_status)[n] == ODD) {
     334          _status->set(n, MATCHED);
     335          oddToMatched(n);
     336        }
     337      }
     338      _tree_set->eraseClass(tree);
     339    }
     340
     341    void unmatchNode(const Node& node) {
     342      int tree = _tree_set->find(node);
     343
     344      alternatePath(node, tree);
     345      destroyTree(tree);
     346
     347      _matching->set(node, INVALID);
     348    }
     349
     350    void augmentOnEdge(const Edge& edge) {
     351      Node left = _bpgraph.u(edge);
     352      int left_tree = _tree_set->find(left);
     353
     354      alternatePath(left, left_tree);
     355      destroyTree(left_tree);
     356      _matching->set(left, _bpgraph.direct(edge, true));
     357
     358      Node right = _bpgraph.v(edge);
     359      int right_tree = _tree_set->find(right);
     360
     361      alternatePath(right, right_tree);
     362      destroyTree(right_tree);
     363      _matching->set(right, _bpgraph.direct(edge, false));
     364    }
     365
     366    void augmentOnArc(const Arc& arc) {
     367      Node left = _bpgraph.source(arc);
     368      _status->set(left, MATCHED);
     369      _matching->set(left, arc);
     370      _pred->set(left, arc);
     371
     372      Node right = _bpgraph.target(arc);
     373      int right_tree = _tree_set->find(right);
     374
     375      alternatePath(right, right_tree);
     376      destroyTree(right_tree);
     377      _matching->set(right, _bpgraph.oppositeArc(arc));
     378    }
     379
     380    void extendOnArc(const Arc& arc) {
     381      Node base = _bpgraph.target(arc);
     382      int tree = _tree_set->find(base);
     383
     384      Node odd = _bpgraph.source(arc);
     385      _tree_set->insert(odd, tree);
     386      _status->set(odd, ODD);
     387      matchedToOdd(odd, tree);
     388      _pred->set(odd, arc);
     389
     390      Node even = _bpgraph.target((*_matching)[odd]);
     391      _tree_set->insert(even, tree);
     392      _status->set(even, EVEN);
     393      matchedToEven(even, tree);
     394    }
     395
     396  public:
     397
     398    /// \brief Constructor
     399    ///
     400    /// Constructor.
     401    MaxWeightedBpMatching(const BpGraph& bpgraph, const WeightMap& weight)
     402      : _bpgraph(bpgraph), _weight(weight), _matching(0),
     403        _node_potential(0), _node_num(0),
     404        _status(0), _pred(0),
     405        _tree_set_index(0), _tree_set(0),
     406
     407        _delta1_index(0), _delta1(0),
     408        _delta2_index(0), _delta2(0),
     409        _delta3_index(0), _delta3(0),
     410
     411        _delta_sum(), _unmatched(0)
     412    {}
     413
     414    ~MaxWeightedBpMatching() {
     415      destroyStructures();
     416    }
     417
     418    /// \name Execution Control
     419    /// The simplest way to execute the algorithm is to use the
     420    /// \ref run() member function.
     421
     422    ///@{
     423
     424    /// \brief Initialize the algorithm
     425    ///
     426    /// This function initializes the algorithm.
     427    void init() {
     428      createStructures();
     429
     430      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     431        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     432        (*_delta2_index)[n] = _delta2->PRE_HEAP;
     433      }
     434      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     435        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     436      }
     437
     438      _delta1->clear();
     439      _delta2->clear();
     440      _delta3->clear();
     441      _tree_set->clear();
     442
     443      _unmatched = _node_num;
     444
     445      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     446        Value max = 0;
     447        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     448          if ((dualScale * _weight[a]) / 2 > max) {
     449            max = (dualScale * _weight[a]) / 2;
     450          }
     451        }
     452        _node_potential->set(n, max);
     453        _delta1->push(n, max);
     454
     455        _tree_set->insert(n);
     456
     457        _matching->set(n, INVALID);
     458        _status->set(n, EVEN);
     459      }
     460      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     461        _delta3->push(e, ((*_node_potential)[_bpgraph.u(e)] +
     462                          (*_node_potential)[_bpgraph.v(e)] -
     463                          dualScale * _weight[e]) / 2);
     464      }
     465    }
     466
     467    /// \brief Initialize the algorithm
     468    ///
     469    /// This function initializes the algorithm.
     470    void redRootInit() {
     471      createStructures();
     472
     473      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     474        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     475        (*_delta2_index)[n] = _delta2->PRE_HEAP;
     476      }
     477      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     478        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     479      }
     480
     481      _delta1->clear();
     482      _delta2->clear();
     483      _delta3->clear();
     484      _tree_set->clear();
     485
     486      _unmatched = 0;
     487      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     488        Value max = 0;
     489        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     490          if ((dualScale * _weight[a]) > max) {
     491            max = dualScale * _weight[a];
     492          }
     493        }
     494        _node_potential->set(n, max);
     495        _delta1->push(n, max);
     496
     497        _tree_set->insert(n);
     498
     499        _matching->set(n, INVALID);
     500        _status->set(n, EVEN);
     501
     502        ++_unmatched;
     503      }
     504      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     505        _matching->set(n, INVALID);
     506        _status->set(n, MATCHED);
     507
     508        Arc min = INVALID;
     509        Value minrw = std::numeric_limits<Value>::max();
     510        for (InArcIt a(_bpgraph, n); a != INVALID; ++a) {
     511          Node v = _bpgraph.source(a);
     512          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
     513                     dualScale * _weight[a];
     514
     515          if (minrw > rw) {
     516            min = _bpgraph.oppositeArc(a);
     517            minrw = rw;
     518          }
     519        }
     520        if (min != INVALID) {
     521          _pred->set(n, min);
     522          _delta2->push(n, minrw);
     523        } else {
     524          _pred->set(n, INVALID);
     525        }
     526      }
     527    }
     528
     529    /// \brief Initialize the algorithm
     530    ///
     531    /// This function initializes the algorithm.
     532    void blueRootInit() {
     533      createStructures();
     534
     535      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     536        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     537        (*_delta2_index)[n] = _delta2->PRE_HEAP;
     538      }
     539      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     540        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     541      }
     542
     543      _delta1->clear();
     544      _delta2->clear();
     545      _delta3->clear();
     546      _tree_set->clear();
     547
     548      _unmatched = 0;
     549      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     550        Value max = 0;
     551        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     552          if ((dualScale * _weight[a]) > max) {
     553            max = dualScale * _weight[a];
     554          }
     555        }
     556        _node_potential->set(n, max);
     557        _delta1->push(n, max);
     558
     559        _tree_set->insert(n);
     560
     561        _matching->set(n, INVALID);
     562        _status->set(n, EVEN);
     563
     564        ++_unmatched;
     565      }
     566      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     567        _matching->set(n, INVALID);
     568        _status->set(n, MATCHED);
     569
     570        Arc min = INVALID;
     571        Value minrw = std::numeric_limits<Value>::max();
     572        for (InArcIt a(_bpgraph, n); a != INVALID; ++a) {
     573          Node v = _bpgraph.source(a);
     574          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
     575                     dualScale * _weight[a];
     576
     577          if (minrw > rw) {
     578            min = _bpgraph.oppositeArc(a);
     579            minrw = rw;
     580          }
     581        }
     582        if (min != INVALID) {
     583          _pred->set(n, min);
     584          _delta2->push(n, minrw);
     585        } else {
     586          _pred->set(n, INVALID);
     587        }
     588      }
     589    }
     590
     591    /// \brief Start the algorithm
     592    ///
     593    /// This function starts the algorithm.
     594    ///
     595    /// \pre \ref init() must be called before using this function.
     596    void start() {
     597      enum OpType {
     598        D1, D2, D3
     599      };
     600
     601      while (_unmatched > 0) {
     602        Value d1 = !_delta1->empty() ?
     603          _delta1->prio() : std::numeric_limits<Value>::max();
     604
     605        Value d2 = !_delta2->empty() ?
     606          _delta2->prio() : std::numeric_limits<Value>::max();
     607
     608        Value d3 = !_delta3->empty() ?
     609          _delta3->prio() : std::numeric_limits<Value>::max();
     610
     611        _delta_sum = d3; OpType ot = D3;
     612        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
     613        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
     614
     615        switch (ot) {
     616        case D1:
     617          {
     618            Node n = _delta1->top();
     619            unmatchNode(n);
     620            --_unmatched;
     621          }
     622          break;
     623        case D2:
     624          {
     625            Node n = _delta2->top();
     626            Arc a = (*_pred)[n];
     627            if ((*_matching)[n] == INVALID) {
     628              augmentOnArc(a);
     629              --_unmatched;
     630            } else {
     631              extendOnArc(a);
     632            }
     633          }
     634          break;
     635        case D3:
     636          {
     637            Edge e = _delta3->top();
     638            augmentOnEdge(e);
     639            _unmatched -= 2;
     640          }
     641          break;
     642        }
     643      }
     644    }
     645
     646    /// \brief Run the algorithm.
     647    ///
     648    /// This method runs the \c %MaxWeightedBpMatching algorithm.
     649    ///
     650    /// \note mwbpm.run() is just a shortcut of the following code.
     651    /// \code
     652    ///   mwbpm.init();
     653    ///   mwbpm.start();
     654    /// \endcode
     655    void run() {
     656      init();
     657      start();
     658    }
     659
     660    /// @}
     661
     662    /// \name Primal Solution
     663    /// Functions to get the primal solution, i.e. the maximum weighted
     664    /// bipartite matching.\n
     665    /// Either \ref run() or \ref start() function should be called before
     666    /// using them.
     667
     668    /// @{
     669
     670    /// \brief Return the weight of the matching.
     671    ///
     672    /// This function returns the weight of the found matching.
     673    ///
     674    /// \pre Either run() or start() must be called before using this function.
     675    Value matchingWeight() const {
     676      Value sum = 0;
     677      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     678        if ((*_matching)[n] != INVALID) {
     679          sum += _weight[(*_matching)[n]];
     680        }
     681      }
     682      return sum;
     683    }
     684
     685    /// \brief Return the size (cardinality) of the matching.
     686    ///
     687    /// This function returns the size (cardinality) of the found matching.
     688    ///
     689    /// \pre Either run() or start() must be called before using this function.
     690    int matchingSize() const {
     691      int num = 0;
     692      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     693        if ((*_matching)[n] != INVALID) {
     694          ++num;
     695        }
     696      }
     697      return num;
     698    }
     699
     700    /// \brief Return \c true if the given edge is in the matching.
     701    ///
     702    /// This function returns \c true if the given edge is in the found
     703    /// matching.
     704    ///
     705    /// \pre Either run() or start() must be called before using this function.
     706    bool matching(const Edge& edge) const {
     707      return edge == (*_matching)[_bpgraph.u(edge)];
     708    }
     709
     710    /// \brief Return the matching arc (or edge) incident to the given node.
     711    ///
     712    /// This function returns the matching arc (or edge) incident to the
     713    /// given node in the found matching or \c INVALID if the node is
     714    /// not covered by the matching.
     715    ///
     716    /// \pre Either run() or start() must be called before using this function.
     717    Arc matching(const Node& node) const {
     718      return (*_matching)[node];
     719    }
     720
     721    /// \brief Return the mate of the given node.
     722    ///
     723    /// This function returns the mate of the given node in the found
     724    /// matching or \c INVALID if the node is not covered by the matching.
     725    ///
     726    /// \pre Either run() or start() must be called before using this function.
     727    Node mate(const Node& node) const {
     728      return (*_matching)[node] != INVALID ?
     729        _bpgraph.target((*_matching)[node]) : INVALID;
     730    }
     731
     732    /// @}
     733
     734    /// \name Dual Solution
     735    /// Functions to get the dual solution.\n
     736    /// Either \ref run() or \ref start() function should be called before
     737    /// using them.
     738
     739    /// @{
     740
     741    /// \brief Return the value of the dual solution.
     742    ///
     743    /// This function returns the value of the dual solution.
     744    /// It should be equal to the primal value scaled by \ref dualScale
     745    /// "dual scale".
     746    ///
     747    /// \pre Either run() or start() must be called before using this function.
     748    Value dualValue() const {
     749      Value sum = 0;
     750      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     751        sum += nodeValue(n);
     752      }
     753      return sum;
     754    }
     755
     756    /// \brief Return the dual value (potential) of the given node.
     757    ///
     758    /// This function returns the dual value (potential) of the given node.
     759    ///
     760    /// \pre Either run() or start() must be called before using this function.
     761    Value nodeValue(const Node& n) const {
     762      return (*_node_potential)[n];
     763    }
     764
     765    /// @}
     766  };
     767
     768  /// \ingroup matching
     769  ///
     770  /// \brief Weighted perfect matching in bipartite graphs
     771  ///
     772  /// This class provides an efficient implementation of multiple search tree
     773  /// augmenting path matching algorithm. The implementation is based on
     774  /// extensive use of priority queues and provides \f$O(nm\log n)\f$ time
     775  /// complexity.
     776  ///
     777  /// The maximum weighted matching problem is to find a subset of the
     778  /// edges in a bipartite graph with maximum overall weight for which
     779  /// each node has exactly one incident edge.
     780  /// It can be formulated with the following linear program.
     781  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
     782  /// \f[x_e \ge 0 \quad \forall e\in E\f]
     783  /// \f[\max \sum_{e\in E}x_ew_e\f]
     784  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
     785  /// \f$X\f$.
     786  ///
     787  /// The algorithm calculates an optimal matching and a proof of the
     788  /// optimality. The solution of the dual problem can be used to check
     789  /// the result of the algorithm. The dual linear problem is the
     790  /// following.
     791  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
     792  /// \f[\min \sum_{u \in V}y_u \f]
     793  /// \tparam BPGR The bipartite graph type the algorithm runs on.
     794  /// \tparam WM The type edge weight map. The default type is
     795  /// \ref concepts::BpGraph::EdgeMap "BPGR::EdgeMap<int>".
     796#ifdef DOXYGEN
     797  template <typename BPGR, typename WM>
     798#else
     799  template <typename BPGR,
     800            typename WM = typename BPGR::template EdgeMap<int> >
     801#endif
     802  class MaxWeightedPerfectBpMatching {
     803  public:
     804
     805    /// The graph type of the algorithm
     806    typedef BPGR BpGraph;
     807    /// The type of the edge weight map
     808    typedef WM WeightMap;
     809    /// The value type of the edge weights
     810    typedef typename WeightMap::Value Value;
     811
     812    /// The type of the matching map
     813    typedef typename BpGraph::template NodeMap<typename BpGraph::Arc>
     814    MatchingMap;
     815
     816    /// \brief Scaling factor for dual solution
     817    ///
     818    /// Scaling factor for dual solution. It is equal to 4 or 1
     819    /// according to the value type.
     820    static const int dualScale =
     821      std::numeric_limits<Value>::is_integer ? 4 : 1;
     822
     823  private:
     824
     825    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     826
     827    typedef typename BpGraph::template NodeMap<Value> NodePotential;
     828
     829    const BpGraph& _bpgraph;
     830    const WeightMap& _weight;
     831
     832    MatchingMap* _matching;
     833
     834    NodePotential* _node_potential;
     835
     836    int _node_num;
     837
     838    enum Status {
     839      EVEN = -1, MATCHED = 0, ODD = 1
     840    };
     841
     842    typedef typename BpGraph::template NodeMap<Status> StatusMap;
     843    StatusMap* _status;
     844
     845    typedef typename BpGraph::template NodeMap<Arc> PredMap;
     846    PredMap* _pred;
     847
     848    typedef ExtendFindEnum<IntNodeMap> TreeSet;
     849    IntNodeMap *_tree_set_index;
     850    TreeSet *_tree_set;
     851
     852    IntNodeMap *_delta2_index;
     853    BinHeap<Value, IntNodeMap> *_delta2;
     854
     855    IntEdgeMap *_delta3_index;
     856    BinHeap<Value, IntEdgeMap> *_delta3;
     857
     858    Value _delta_sum;
     859    int _unmatched;
     860
     861    void createStructures() {
     862      _node_num = countNodes(_bpgraph);
     863
     864      if (!_matching) {
     865        _matching = new MatchingMap(_bpgraph);
     866      }
     867
     868      if (!_node_potential) {
     869        _node_potential = new NodePotential(_bpgraph);
     870      }
     871
     872      if (!_status) {
     873        _status = new StatusMap(_bpgraph);
     874      }
     875
     876      if (!_pred) {
     877        _pred = new PredMap(_bpgraph);
     878      }
     879
     880      if (!_tree_set) {
     881        _tree_set_index = new IntNodeMap(_bpgraph);
     882        _tree_set = new TreeSet(*_tree_set_index);
     883      }
     884
     885      if (!_delta2) {
     886        _delta2_index = new IntNodeMap(_bpgraph);
     887        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
     888      }
     889
     890      if (!_delta3) {
     891        _delta3_index = new IntEdgeMap(_bpgraph);
     892        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
     893      }
     894    }
     895
     896    void destroyStructures() {
     897      if (_matching) {
     898        delete _matching;
     899      }
     900      if (_node_potential) {
     901        delete _node_potential;
     902      }
     903      if (_status) {
     904        delete _status;
     905      }
     906      if (_pred) {
     907        delete _pred;
     908      }
     909      if (_tree_set) {
     910        delete _tree_set_index;
     911        delete _tree_set;
     912      }
     913      if (_delta2) {
     914        delete _delta2_index;
     915        delete _delta2;
     916      }
     917      if (_delta3) {
     918        delete _delta3_index;
     919        delete _delta3;
     920      }
     921    }
     922
     923    void matchedToEven(Node node, int tree) {
     924      _tree_set->insert(node, tree);
     925      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     926
     927      if (_delta2->state(node) == _delta2->IN_HEAP) {
     928        _delta2->erase(node);
     929      }
     930
     931      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     932        Node v = _bpgraph.source(a);
     933        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     934          dualScale * _weight[a];
     935        if ((*_status)[v] == EVEN) {
     936          _delta3->push(a, rw / 2);
     937        } else if ((*_status)[v] == MATCHED) {
     938          if (_delta2->state(v) != _delta2->IN_HEAP) {
     939            _pred->set(v, a);
     940            _delta2->push(v, rw);
     941          } else if ((*_delta2)[v] > rw) {
     942            _pred->set(v, a);
     943            _delta2->decrease(v, rw);
     944          }
     945        }
     946      }
     947    }
     948
     949    void matchedToOdd(Node node, int tree) {
     950      _tree_set->insert(node, tree);
     951      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     952
     953      if (_delta2->state(node) == _delta2->IN_HEAP) {
     954        _delta2->erase(node);
     955      }
     956    }
     957
     958    void evenToMatched(Node node, int tree) {
     959      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     960      Arc min = INVALID;
     961      Value minrw = std::numeric_limits<Value>::max();
     962      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     963        Node v = _bpgraph.source(a);
     964        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     965          dualScale * _weight[a];
     966
     967        if ((*_status)[v] == EVEN) {
     968          _delta3->erase(a);
     969          if (minrw > rw) {
     970            min = _bpgraph.oppositeArc(a);
     971            minrw = rw;
     972          }
     973        } else if ((*_status)[v]  == MATCHED) {
     974          if ((*_pred)[v] == a) {
     975            Arc mina = INVALID;
     976            Value minrwa = std::numeric_limits<Value>::max();
     977            for (OutArcIt aa(_bpgraph, v); aa != INVALID; ++aa) {
     978              Node va = _bpgraph.target(aa);
     979              if ((*_status)[va] != EVEN ||
     980                  _tree_set->find(va) == tree) continue;
     981              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
     982                dualScale * _weight[aa];
     983              if (minrwa > rwa) {
     984                minrwa = rwa;
     985                mina = aa;
     986              }
     987            }
     988            if (mina != INVALID) {
     989              _pred->set(v, mina);
     990              _delta2->increase(v, minrwa);
     991            } else {
     992              _pred->set(v, INVALID);
     993              _delta2->erase(v);
     994            }
     995          }
     996        }
     997      }
     998      if (min != INVALID) {
     999        _pred->set(node, min);
     1000        _delta2->push(node, minrw);
     1001      } else {
     1002        _pred->set(node, INVALID);
     1003      }
     1004    }
     1005
     1006    void oddToMatched(Node node) {
     1007      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     1008      Arc min = INVALID;
     1009      Value minrw = std::numeric_limits<Value>::max();
     1010      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     1011        Node v = _bpgraph.source(a);
     1012        if ((*_status)[v] != EVEN) continue;
     1013        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     1014          dualScale * _weight[a];
     1015
     1016        if (minrw > rw) {
     1017          min = _bpgraph.oppositeArc(a);
     1018          minrw = rw;
     1019        }
     1020      }
     1021      if (min != INVALID) {
     1022        _pred->set(node, min);
     1023        _delta2->push(node, minrw);
     1024      } else {
     1025        _pred->set(node, INVALID);
     1026      }
     1027    }
     1028
     1029    void alternatePath(Node even, int tree) {
     1030      Node odd;
     1031
     1032      _status->set(even, MATCHED);
     1033      evenToMatched(even, tree);
     1034
     1035      Arc prev = (*_matching)[even];
     1036      while (prev != INVALID) {
     1037        odd = _bpgraph.target(prev);
     1038        even = _bpgraph.target((*_pred)[odd]);
     1039        _matching->set(odd, (*_pred)[odd]);
     1040        _status->set(odd, MATCHED);
     1041        oddToMatched(odd);
     1042
     1043        prev = (*_matching)[even];
     1044        _status->set(even, MATCHED);
     1045        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
     1046        evenToMatched(even, tree);
     1047      }
     1048    }
     1049
     1050    void destroyTree(int tree) {
     1051      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
     1052        if ((*_status)[n] == EVEN) {
     1053          _status->set(n, MATCHED);
     1054          evenToMatched(n, tree);
     1055        } else if ((*_status)[n] == ODD) {
     1056          _status->set(n, MATCHED);
     1057          oddToMatched(n);
     1058        }
     1059      }
     1060      _tree_set->eraseClass(tree);
     1061    }
     1062
     1063    void unmatchNode(const Node& node) {
     1064      int tree = _tree_set->find(node);
     1065
     1066      alternatePath(node, tree);
     1067      destroyTree(tree);
     1068
     1069      _matching->set(node, INVALID);
     1070    }
     1071
     1072    void augmentOnEdge(const Edge& edge) {
     1073      Node left = _bpgraph.u(edge);
     1074      int left_tree = _tree_set->find(left);
     1075
     1076      alternatePath(left, left_tree);
     1077      destroyTree(left_tree);
     1078      _matching->set(left, _bpgraph.direct(edge, true));
     1079
     1080      Node right = _bpgraph.v(edge);
     1081      int right_tree = _tree_set->find(right);
     1082
     1083      alternatePath(right, right_tree);
     1084      destroyTree(right_tree);
     1085      _matching->set(right, _bpgraph.direct(edge, false));
     1086    }
     1087
     1088    void augmentOnArc(const Arc& arc) {
     1089      Node left = _bpgraph.source(arc);
     1090      _status->set(left, MATCHED);
     1091      _matching->set(left, arc);
     1092      _pred->set(left, arc);
     1093
     1094      Node right = _bpgraph.target(arc);
     1095      int right_tree = _tree_set->find(right);
     1096
     1097      alternatePath(right, right_tree);
     1098      destroyTree(right_tree);
     1099      _matching->set(right, _bpgraph.oppositeArc(arc));
     1100    }
     1101
     1102    void extendOnArc(const Arc& arc) {
     1103      Node base = _bpgraph.target(arc);
     1104      int tree = _tree_set->find(base);
     1105
     1106      Node odd = _bpgraph.source(arc);
     1107      _tree_set->insert(odd, tree);
     1108      _status->set(odd, ODD);
     1109      matchedToOdd(odd, tree);
     1110      _pred->set(odd, arc);
     1111
     1112      Node even = _bpgraph.target((*_matching)[odd]);
     1113      _tree_set->insert(even, tree);
     1114      _status->set(even, EVEN);
     1115      matchedToEven(even, tree);
     1116    }
     1117
     1118  public:
     1119
     1120    /// \brief Constructor
     1121    ///
     1122    /// Constructor.
     1123    MaxWeightedPerfectBpMatching(const BpGraph& bpgraph, const WeightMap& weight)
     1124      : _bpgraph(bpgraph), _weight(weight), _matching(0),
     1125        _node_potential(0), _node_num(0),
     1126        _status(0), _pred(0),
     1127        _tree_set_index(0), _tree_set(0),
     1128
     1129        _delta2_index(0), _delta2(0),
     1130        _delta3_index(0), _delta3(0),
     1131
     1132        _delta_sum(), _unmatched(0)
     1133    {}
     1134
     1135    ~MaxWeightedPerfectBpMatching() {
     1136      destroyStructures();
     1137    }
     1138
     1139    /// \name Execution Control
     1140    /// The simplest way to execute the algorithm is to use the
     1141    /// \ref run() member function.
     1142
     1143    ///@{
     1144
     1145    /// \brief Initialize the algorithm
     1146    ///
     1147    /// This function initializes the algorithm.
     1148    ///
     1149    /// \return If it is false, then the graph does not have a perfect matching.
     1150    bool init() {
     1151      createStructures();
     1152
     1153      if (countRedNodes(_bpgraph) != countBlueNodes(_bpgraph)) {
     1154        return false;
     1155      }
     1156
     1157      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1158        (*_delta2_index)[n] = _delta2->PRE_HEAP;
     1159      }
     1160      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     1161        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     1162      }
     1163
     1164      _delta2->clear();
     1165      _delta3->clear();
     1166      _tree_set->clear();
     1167
     1168      _unmatched = _node_num;
     1169
     1170      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1171        Value max = - std::numeric_limits<Value>::max();
     1172        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     1173          if ((dualScale * _weight[a]) / 2 > max) {
     1174            max = (dualScale * _weight[a]) / 2;
     1175          }
     1176        }
     1177        if (max == - std::numeric_limits<Value>::max()) {
     1178          return false;
     1179        }
     1180        _node_potential->set(n, max);
     1181
     1182        _tree_set->insert(n);
     1183
     1184        _matching->set(n, INVALID);
     1185        _status->set(n, EVEN);
     1186      }
     1187      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     1188        _delta3->push(e, ((*_node_potential)[_bpgraph.u(e)] +
     1189                          (*_node_potential)[_bpgraph.v(e)] -
     1190                          dualScale * _weight[e]) / 2);
     1191      }
     1192      return true;
     1193    }
     1194
     1195    /// \brief Start the algorithm
     1196    ///
     1197    /// This function starts the algorithm.
     1198    ///
     1199    /// \pre \ref init() must be called before using this function.
     1200    ///
     1201    /// \return True when a perfect matching is found.
     1202    bool start() {
     1203      enum OpType {
     1204        D2, D3
     1205      };
     1206
     1207      while (_unmatched > 0) {
     1208        Value d2 = !_delta2->empty() ?
     1209          _delta2->prio() : std::numeric_limits<Value>::max();
     1210
     1211        Value d3 = !_delta3->empty() ?
     1212          _delta3->prio() : std::numeric_limits<Value>::max();
     1213
     1214        _delta_sum = d3; OpType ot = D3;
     1215        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
     1216
     1217        if (_delta_sum == std::numeric_limits<Value>::max()) {
     1218          return false;
     1219        }
     1220
     1221        switch (ot) {
     1222        case D2:
     1223          {
     1224            Node n = _delta2->top();
     1225            Arc a = (*_pred)[n];
     1226            if ((*_matching)[n] == INVALID) {
     1227              augmentOnArc(a);
     1228              --_unmatched;
     1229            } else {
     1230              extendOnArc(a);
     1231            }
     1232          }
     1233          break;
     1234        case D3:
     1235          {
     1236            Edge e = _delta3->top();
     1237            augmentOnEdge(e);
     1238            _unmatched -= 2;
     1239          }
     1240          break;
     1241        }
     1242      }
     1243      return true;
     1244    }
     1245
     1246    /// \brief Run the algorithm.
     1247    ///
     1248    /// This method runs the \c %MaxWeightedPerfectBpMatching algorithm.
     1249    ///
     1250    /// \note mwpbpm.run() is just a shortcut of the following code.
     1251    /// \code
     1252    ///   return mwpbpm.init() && mwpbpm.start();
     1253    /// \endcode
     1254    ///
     1255    /// \return True when a perfect matching is found.
     1256    bool run() {
     1257      return init() && start();
     1258    }
     1259
     1260    /// @}
     1261
     1262    /// \name Primal Solution
     1263    /// Functions to get the primal solution, i.e. the maximum weighted
     1264    /// bipartite matching.\n
     1265    /// Either \ref run() or \ref start() function should be called before
     1266    /// using them and a matching should be found.
     1267
     1268    /// @{
     1269
     1270    /// \brief Return the weight of the matching.
     1271    ///
     1272    /// This function returns the weight of the found matching.
     1273    ///
     1274    /// \pre A perfect matching has been found.
     1275    Value matchingWeight() const {
     1276      Value sum = 0;
     1277      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     1278        sum += _weight[(*_matching)[n]];
     1279      }
     1280      return sum;
     1281    }
     1282
     1283    /// \brief Return the size (cardinality) of the matching.
     1284    ///
     1285    /// This function returns the size (cardinality) of the found matching.
     1286    ///
     1287    /// \pre A perfect matching has been found.
     1288    int matchingSize() const {
     1289      int num = 0;
     1290      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     1291          ++num;
     1292      }
     1293      return num;
     1294    }
     1295
     1296    /// \brief Return \c true if the given edge is in the matching.
     1297    ///
     1298    /// This function returns \c true if the given edge is in the found
     1299    /// matching.
     1300    ///
     1301    /// \pre A perfect matching has been found.
     1302    bool matching(const Edge& edge) const {
     1303      return edge == (*_matching)[_bpgraph.u(edge)];
     1304    }
     1305
     1306    /// \brief Return the matching arc (or edge) incident to the given node.
     1307    ///
     1308    /// This function returns the matching arc (or edge) incident to the
     1309    /// given node in the found matching or \c INVALID if the node is
     1310    /// not covered by the matching.
     1311    ///
     1312    /// \pre A perfect matching has been found.
     1313    Arc matching(const Node& node) const {
     1314      return (*_matching)[node];
     1315    }
     1316
     1317    /// \brief Return the mate of the given node.
     1318    ///
     1319    /// This function returns the mate of the given node in the found
     1320    /// matching or \c INVALID if the node is not covered by the matching.
     1321    ///
     1322    /// \pre Either run() or start() must be called before using this function.
     1323    Node mate(const Node& node) const {
     1324      return (*_matching)[node] != INVALID ?
     1325        _bpgraph.target((*_matching)[node]) : INVALID;
     1326    }
     1327
     1328    /// @}
     1329
     1330    /// \name Dual Solution
     1331    /// Functions to get the dual solution.\n
     1332    /// Either \ref run() or \ref start() function should be called before
     1333    /// using them and a matching should be found.
     1334
     1335    /// @{
     1336
     1337    /// \brief Return the value of the dual solution.
     1338    ///
     1339    /// This function returns the value of the dual solution.
     1340    /// It should be equal to the primal value scaled by \ref dualScale
     1341    /// "dual scale".
     1342    ///
     1343    /// \pre A perfect matching has been found.
     1344    Value dualValue() const {
     1345      Value sum = 0;
     1346      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1347        sum += nodeValue(n);
     1348      }
     1349      return sum;
     1350    }
     1351
     1352    /// \brief Return the dual value (potential) of the given node.
     1353    ///
     1354    /// This function returns the dual value (potential) of the given node.
     1355    ///
     1356    /// \pre A perfect matching has been found.
     1357    Value nodeValue(const Node& n) const {
     1358      return (*_node_potential)[n];
     1359    }
     1360
     1361    /// @}
     1362  };
     1363
     1364} //END OF NAMESPACE LEMON
     1365
     1366#endif //LEMON_BPMATCHING_H
  • test/CMakeLists.txt

    diff -r 8c567e298d7f -r 5b4317f3ce36 test/CMakeLists.txt
    a b  
    1717  bellman_ford_test
    1818  bfs_test
    1919  bpgraph_test
     20  bpmatching_test
    2021  circulation_test
    2122  connectivity_test
    2223  counter_test
  • new file test/bpmatching_test.cc

    diff -r 8c567e298d7f -r 5b4317f3ce36 test/bpmatching_test.cc
    - +  
     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/bpmatching.h>
     25#include <lemon/matching.h>
     26#include <lemon/smart_graph.h>
     27#include <lemon/concepts/bpgraph.h>
     28#include <lemon/concepts/maps.h>
     29#include <lemon/lgf_reader.h>
     30
     31#include "test_tools.h"
     32
     33using namespace std;
     34using namespace lemon;
     35
     36const int lgfn = 7;
     37const string lgf[lgfn] = {
     38  // Empty graph
     39  "@red_nodes\n"
     40  "label\n"
     41  "\n"
     42  "@blue_nodes\n"
     43  "label\n"
     44  "\n"
     45  "@edges\n"
     46  "  label weight\n"
     47  "\n",
     48
     49  // No red nodes
     50  "@red_nodes\n"
     51  "label\n"
     52  "\n"
     53  "@blue_nodes\n"
     54  "label\n"
     55  "1\n"
     56  "@edges\n"
     57  "  label weight\n"
     58  "\n",
     59
     60  // No blue nodes
     61  "@red_nodes\n"
     62  "label\n"
     63  "1\n"
     64  "\n"
     65  "@blue_nodes\n"
     66  "label\n"
     67  "\n"
     68  "@edges\n"
     69  "  label weight\n"
     70  "\n",
     71
     72  "@red_nodes\n"
     73  "label\n"
     74  "1\n"
     75  "@blue_nodes\n"
     76  "label\n"
     77  "2\n"
     78  "\n"
     79  "@edges\n"
     80  "  label weight\n"
     81  "1 2 1 1\n"
     82  "\n",
     83
     84  "@red_nodes\n"
     85  "label\n"
     86  "1\n"
     87  "2\n"
     88  "\n"
     89  "@blue_nodes\n"
     90  "label\n"
     91  "3\n"
     92  "4\n"
     93  "\n"
     94  "@edges\n"
     95  "  label weight\n"
     96  "1 3 1 5\n"
     97  "1 4 2 3\n"
     98  "2 3 3 1\n"
     99  "2 4 4 2\n"
     100  "\n",
     101
     102  "@red_nodes\n"
     103  "label\n"
     104  "1\n"
     105  "2\n"
     106  "\n"
     107  "@blue_nodes\n"
     108  "label\n"
     109  "3\n"
     110  "4\n"
     111  "\n"
     112  "@edges\n"
     113  "  label weight\n"
     114  "1 3 1 3\n"
     115  "1 4 2 2\n"
     116  "2 3 3 2\n"
     117  "\n",
     118
     119  "@red_nodes\n"
     120  "label\n"
     121  "1\n"
     122  "2\n"
     123  "\n"
     124  "@blue_nodes\n"
     125  "label\n"
     126  "3\n"
     127  "4\n"
     128  "\n"
     129  "@edges\n"
     130  "  label weight\n"
     131  "1 3 1 5\n"
     132  "1 4 2 2\n"
     133  "2 3 3 2\n"
     134  "\n",
     135};
     136
     137void checkMaxWeightedBpMatchingCompiles() {
     138  typedef concepts::BpGraph BpGraph;
     139  typedef BpGraph::Node Node;
     140  typedef BpGraph::Edge Edge;
     141  typedef BpGraph::EdgeMap<int> IntEdgeMap;
     142  typedef MaxWeightedBpMatching<BpGraph, IntEdgeMap> MWBPM;
     143
     144  int v;
     145  Node n;
     146  Edge e;
     147  bool b;
     148
     149  BpGraph graph;
     150  IntEdgeMap weight(graph);
     151  MWBPM mwbpm(graph, weight);
     152  const MWBPM& cmwbpm = mwbpm;
     153
     154  mwbpm.init();
     155  mwbpm.redRootInit();
     156  mwbpm.blueRootInit();
     157  mwbpm.start();
     158  mwbpm.run();
     159
     160  v = cmwbpm.matchingWeight();
     161  v = cmwbpm.matchingSize();
     162  e = cmwbpm.matching(n);
     163  b = cmwbpm.matching(e);
     164  n = cmwbpm.mate(n);
     165  v = cmwbpm.dualValue();
     166  v = cmwbpm.nodeValue(n);
     167
     168  ::lemon::ignore_unused_variable_warning(v, b);
     169}
     170
     171void checkMaxWeightedPerfectBpMatchingCompiles() {
     172  typedef concepts::BpGraph BpGraph;
     173  typedef BpGraph::Node Node;
     174  typedef BpGraph::Edge Edge;
     175  typedef BpGraph::EdgeMap<int> IntEdgeMap;
     176  typedef MaxWeightedPerfectBpMatching<BpGraph, IntEdgeMap> MWPBPM;
     177
     178  int v;
     179  Node n;
     180  Edge e;
     181  bool b;
     182
     183  BpGraph graph;
     184  IntEdgeMap weight(graph);
     185  MWPBPM mwpbpm(graph, weight);
     186  const MWPBPM& cmwpbpm = mwpbpm;
     187
     188  mwpbpm.init();
     189  mwpbpm.start();
     190  mwpbpm.run();
     191
     192  v = cmwpbpm.matchingWeight();
     193  v = cmwpbpm.matchingSize();
     194  e = cmwpbpm.matching(n);
     195  b = cmwpbpm.matching(e);
     196  n = cmwpbpm.mate(n);
     197  v = cmwpbpm.dualValue();
     198  v = cmwpbpm.nodeValue(n);
     199
     200  ::lemon::ignore_unused_variable_warning(v, b);
     201}
     202
     203
     204void checkMaxWeightedBpMatching(
     205    const SmartBpGraph& bpgraph,
     206    const SmartBpGraph::EdgeMap<int>& weight,
     207    const MaxWeightedBpMatching<SmartBpGraph>& mwbpm) {
     208  int pv = 0;
     209  for (SmartBpGraph::EdgeIt e(bpgraph); e != INVALID; ++e) {
     210    int rw = mwbpm.nodeValue(bpgraph.redNode(e)) +
     211             mwbpm.nodeValue(bpgraph.blueNode(e));
     212    rw -= weight[e] * mwbpm.dualScale;
     213    check(rw >= 0, "Negative reduced weight");
     214    check(rw == 0 || !mwbpm.matching(e),
     215          "Non-zero reduced weight on matching edge");
     216    if (mwbpm.matching(e)) {
     217      pv += weight[e];
     218    }
     219  }
     220  int dv = 0;
     221  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
     222    if (mwbpm.matching(n) != INVALID) {
     223      check(mwbpm.nodeValue(n) >= 0, "Invalid node value");
     224      SmartBpGraph::Node o = bpgraph.target(mwbpm.matching(n));
     225      check(mwbpm.mate(n) == o, "Invalid matching");
     226      check(mwbpm.matching(n) == bpgraph.oppositeArc(mwbpm.matching(o)),
     227            "Invalid matching");
     228    } else {
     229      check(mwbpm.mate(n) == INVALID, "Invalid matching");
     230      check(mwbpm.nodeValue(n) == 0, "Invalid matching");
     231    }
     232  }
     233
     234  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
     235    dv += mwbpm.nodeValue(n);
     236  }
     237
     238  check(pv * mwbpm.dualScale == dv, "Wrong duality");
     239}
     240
     241void checkPerfectBpMatchingExists(
     242    const SmartBpGraph& bpgraph,
     243    bool exists) {
     244  if (countRedNodes(bpgraph) == countBlueNodes(bpgraph)) {
     245    // TODO: Use a bipartite maximum cardinality matching algorithm here.
     246    MaxMatching<SmartBpGraph> mbpm(bpgraph);
     247    mbpm.run();
     248    check(exists == (countRedNodes(bpgraph) == mbpm.matchingSize()),
     249          "Unexpected maximum cardinality matching size");
     250  } else {
     251    check(!exists, "Different size for red and blue node set");
     252  }
     253}
     254
     255void checkMaxWeightedPerfectBpMatching(
     256    const SmartBpGraph& bpgraph,
     257    const SmartBpGraph::EdgeMap<int>& weight,
     258    const MaxWeightedPerfectBpMatching<SmartBpGraph>& mwpbpm) {
     259  int pv = 0;
     260  for (SmartBpGraph::EdgeIt e(bpgraph); e != INVALID; ++e) {
     261    int rw = mwpbpm.nodeValue(bpgraph.redNode(e)) +
     262             mwpbpm.nodeValue(bpgraph.blueNode(e));
     263    rw -= weight[e] * mwpbpm.dualScale;
     264    check(rw >= 0, "Negative reduced weight");
     265    check(rw == 0 || !mwpbpm.matching(e),
     266          "Non-zero reduced weight on matching edge");
     267    if (mwpbpm.matching(e)) {
     268      pv += weight[e];
     269    }
     270  }
     271  int dv = 0;
     272  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
     273    check(mwpbpm.matching(n) != INVALID, "Matching is not perfect");
     274    SmartBpGraph::Node o = bpgraph.target(mwpbpm.matching(n));
     275    check(mwpbpm.mate(n) == o, "Invalid matching");
     276    check(mwpbpm.matching(n) == bpgraph.oppositeArc(mwpbpm.matching(o)),
     277          "Invalid matching");
     278  }
     279
     280  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
     281    dv += mwpbpm.nodeValue(n);
     282  }
     283
     284  check(pv * mwpbpm.dualScale == dv, "Wrong duality");
     285}
     286
     287int main() {
     288  // Checking hard wired extremal graphs
     289  for (int i=0; i<lgfn; ++i) {
     290    SmartBpGraph bpgraph;
     291    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     292    istringstream lgfs(lgf[i]);
     293    bpGraphReader(bpgraph, lgfs)
     294        .edgeMap("weight", weight).run();
     295
     296    {
     297      MaxWeightedBpMatching<SmartBpGraph> mwbpm(bpgraph, weight);
     298      mwbpm.run();
     299      checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     300    }
     301    {
     302      MaxWeightedBpMatching<SmartBpGraph> mwbpm(bpgraph, weight);
     303      mwbpm.redRootInit();
     304      mwbpm.start();
     305      checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     306    }
     307    {
     308      MaxWeightedBpMatching<SmartBpGraph> mwbpm(bpgraph, weight);
     309      mwbpm.blueRootInit();
     310      mwbpm.start();
     311      checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     312    }
     313
     314    {
     315      MaxWeightedPerfectBpMatching<SmartBpGraph> mwpbpm(bpgraph, weight);
     316      bool exists = mwpbpm.run();
     317      checkPerfectBpMatchingExists(bpgraph, exists);
     318      if (exists) {
     319        checkMaxWeightedPerfectBpMatching(bpgraph, weight, mwpbpm);
     320      }
     321    }
     322  }
     323
     324  return 0;
     325}