COIN-OR::LEMON - Graph Library

Ticket #168: bp_matching_b08198bdd613.patch

File bp_matching_b08198bdd613.patch, 31.6 KB (added by Mohammed El-Kebir, 11 years ago)

Maximum Weight Bipartite Matching (b08198bdd613)

  • lemon/Makefile.am

    exporting patch:
    # HG changeset patch
    # User Mohammed El-Kebir <m.elkebir@gmail.com>
    # Date 1291391575 -3600
    # Node ID b08198bdd6138ae9e8905579f81756961247a839
    # Parent  20eec7f2b0d6f2dcfcb2b6bab9fb652d777d5a9b
    Implemented two Maximum Weight Bipartite Matching algorithms.
    
    diff --git a/lemon/Makefile.am b/lemon/Makefile.am
    a b  
    6161        lemon/bfs.h \
    6262        lemon/bin_heap.h \
    6363        lemon/binomial_heap.h \
     64        lemon/bp_matching.h \
    6465        lemon/bucket_heap.h \
    6566        lemon/capacity_scaling.h \
    6667        lemon/cbc.h \
  • new file lemon/bp_matching.h

    diff --git a/lemon/bp_matching.h b/lemon/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-2010
     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 BP_MATCHING_H
     20#define BP_MATCHING_H
     21
     22#include <limits>
     23#include <list>
     24#include <algorithm>
     25#include <assert.h>
     26#include <queue>
     27
     28#include <lemon/core.h>
     29#include <lemon/bin_heap.h>
     30
     31///\ingroup matching
     32///\file
     33///\brief Maximum weight matching algorithms in bipartite graphs.
     34
     35namespace lemon {
     36
     37  /// \ingroup matching
     38  ///
     39  /// \brief Maximum weight matching in (sparse) bipartite graphs
     40  ///
     41  /// This class implements a successive shortest path algorithm for finding
     42  /// a maximum weight matching in an undirected bipartite graph.
     43  /// Let \f$G = (X \cup Y, E)\f$ be an undirected bipartite graph. The
     44  /// following linear program corresponds to a maximum weight matching
     45  /// in the graph \f$G\f$.
     46  ///
     47  /** \f$\begin{array}{rrcll} \
     48      \max & \displaystyle\sum_{(i,j) \in E} c_{ij} x_{ij}\\ \
     49      \mbox{s.t.} & \displaystyle\sum_{i \in X} x_{ij} & \leq & 1, \
     50      & \forall j \in \{ j^\prime \in Y \mid (i,j^\prime) \in E \}\\ \
     51      & \displaystyle\sum_{j \in Y} x_{ij} & \leq & 1, \
     52      & \forall i \in \{ i^\prime \in X \mid (i^\prime,j) \in E \}\\ \
     53      & x_{ij}        & \geq & 0, & \forall (i,j) \in E\\\end{array}\f$
     54  */
     55  ///
     56  /// where \f$c_{ij}\f$ is the weight of edge \f$(i,j)\f$. The dual problem
     57  /// is:
     58  ///
     59  /** \f$\begin{array}{rrcll}\min & \displaystyle\sum_{v \in X \cup Y} p_v\\ \
     60      \mbox{s.t.} & p_i + p_j & \geq & c_{ij}, & \forall (i,j) \in E\\ \
     61      & p_v & \geq & 0, & \forall v \in X \cup Y \end{array}\f$
     62  */
     63  ///
     64  /// A maximum weight matching is constructed by iteratively considering the
     65  /// vertices in \f$X = \{x_1, \ldots, x_n\}\f$. In every iteration \f$k\f$
     66  /// we establish primal and dual complementary slackness for the subgraph
     67  /// \f$G[X_k \cup Y]\f$ where \f$X_k = \{x_1, \ldots, x_k\}\f$.
     68  /// So after the final iteration the primal and dual solution will be equal,
     69  /// and we will thus have a maximum weight matching. The time complexity of
     70  /// this method is \f$O(n(n + m)\log n)\f$.
     71  ///
     72  /// In case the bipartite graph is dense, it is better to use
     73  /// \ref MaxWeightedDenseBipartiteMatching, which has a time complexity of
     74  /// \f$O(n^3)\f$.
     75  ///
     76  /// \tparam BGR The bipartite graph type the algorithm runs on.
     77  /// \tparam WM The type edge weight map. The default type is
     78  /// \ref concepts::Graph::EdgeMap "BGR::EdgeMap<int>".
     79#ifdef DOXYGEN
     80  template <typename BGR, typename WM>
     81#else
     82  template <typename BGR,
     83            typename WM = typename BGR::template EdgeMap<int> >
     84#endif
     85  class MaxWeightedBipartiteMatching
     86  {
     87  public:
     88    /// The graph type of the algorithm
     89    typedef BGR BpGraph;
     90    /// The type of the edge weight map
     91    typedef WM WeightMap;
     92    /// The value type of the edge weights
     93    typedef typename WeightMap::Value Value;
     94    /// The type of the matching map
     95    typedef typename BpGraph::
     96      template NodeMap<typename BpGraph::Edge> MatchingMap;
     97
     98  private:
     99    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     100    typedef typename BpGraph::template NodeMap<Value> PotMap;
     101    typedef std::list<RedNode> RedNodeList;
     102    typedef std::list<BlueNode> BlueNodeList;
     103    typedef typename BpGraph::template NodeMap<Value> DistMap;
     104    typedef typename BpGraph::template BlueMap<int> HeapCrossRef;
     105    typedef BinHeap<Value, HeapCrossRef> Heap;
     106    typedef typename BpGraph::template NodeMap<Arc> PredMap;
     107
     108    const BpGraph& _graph;
     109    const WeightMap& _weight;
     110    PotMap* _pPot;
     111    MatchingMap* _pMatchingMap;
     112    Value _matchingWeight;
     113    int _matchingSize;
     114
     115    void createStructures()
     116    {
     117      _pPot = new PotMap(_graph, 0);
     118      _pMatchingMap = new MatchingMap(_graph, INVALID);
     119    }
     120
     121    void destroyStructures()
     122    {
     123      delete _pPot;
     124      delete _pMatchingMap;
     125    }
     126
     127    bool isFree(const Node& v)
     128    {
     129      return (*_pMatchingMap)[v] == INVALID;
     130    }
     131
     132    void augmentPath(Arc a, bool matching, const PredMap& pred)
     133    {
     134      // M' = M ^ EP
     135      while (a != INVALID)
     136      {
     137        if (!matching)
     138        {
     139          _pMatchingMap->set(_graph.source(a), a);
     140          _pMatchingMap->set(_graph.target(a), a);
     141        }
     142
     143        matching = !matching;
     144        a = pred[_graph.source(a)];
     145      }
     146    }
     147
     148    void augment(const Node& x, DistMap& dist, PredMap& pred)
     149    {
     150      assert(isFree(x));
     151
     152      /**
     153       * In case maxCardinality == false, we also need to consider
     154       * augmenting paths starting from x and ending in a matched
     155       * node x' in X. Augmenting such a path does *not* increase
     156       * the cardinality of the matching. It may, however, increase
     157       * the weight of the matching.
     158       *
     159       * Along with a shortest path starting from x and ending in
     160       * a free vertex y in Y, we also determine x' such that
     161       *   y' = pred[x'],
     162       *   (pot[x] + pot[y'] - dist[x, y']) - w(y', x') is maximal
     163       *
     164       * Since (y', x') is part of the matching,
     165       * by primal complementary slackness we have that
     166       *   pot[y'] + pot[x'] = w(y', x').
     167       *
     168       * Hence
     169       * x' = arg max_{x' \in X} { pot[x] + pot[y'] - dist[x, y']) -w(y', x') }
     170       *    = arg max_{x' \in X} { pot[x] - dist[x, y'] - pot[x'] }
     171       *    = arg max_{x' \in X} { -dist[x, y'] - pot[x'] }
     172       *    = arg min_{x' \in X} { dist[x, y'] + [x'] }
     173       *
     174       * We only augment x ->* x' if dist(x,y) > dist[x, y'] + pot[x']
     175       * Otherwise we augment x ->* y.
     176       */
     177
     178      Value UB = (*_pPot)[x];
     179      dist[x] = 0;
     180
     181      RedNodeList visitedX;
     182      BlueNodeList visitedY;
     183
     184      // heap only contains nodes in Y
     185      HeapCrossRef heapCrossRef(_graph, Heap::PRE_HEAP);
     186      Heap heap(heapCrossRef);
     187
     188      // add nodes adjacent to x to heap, and update UB
     189      visitedX.push_back(x);
     190      for (IncEdgeIt e(_graph, x); e != INVALID; ++e)
     191      {
     192        const BlueNode y = _graph.blueNode(e);
     193        Value dist_y = (*_pPot)[x] + (*_pPot)[y] - _weight[e];
     194
     195        if (dist_y >= UB)
     196          continue;
     197
     198        if (isFree(y))
     199          UB = dist_y;
     200
     201        dist[y] = dist_y;
     202        pred[y] = _graph.direct(e, x);
     203
     204        assert(heap.state(y) == Heap::PRE_HEAP);
     205        heap.push(y, dist_y);
     206      }
     207
     208      Node x_min = x;
     209      Value min_dist = 0, x_min_dist = (*_pPot)[x];
     210
     211      while (true)
     212      {
     213        assert(heap.empty() || heap.prio() == dist[heap.top()]);
     214
     215        if (heap.empty() || heap.prio() >= x_min_dist)
     216        {
     217          min_dist = x_min_dist;
     218
     219          if (x_min != x)
     220          {
     221            // we have an augmenting path between x and x_min
     222            // that doesn't increase the matching size
     223            _matchingWeight += (*_pPot)[x] - x_min_dist;
     224
     225            // x_min becomes free, and will always remain free
     226            (*_pMatchingMap)[x_min] = INVALID;
     227            augmentPath(pred[x_min], true, pred);
     228          }
     229          break;
     230        }
     231
     232        const BlueNode y = heap.top();
     233        const Value dist_y = heap.prio();
     234        heap.pop();
     235
     236        visitedY.push_back(y);
     237        if (isFree(y))
     238        {
     239          // we have an augmenting path between x and y
     240          augmentPath(pred[y], false, pred);
     241          _matchingSize++;
     242
     243          assert((*_pPot)[y] == 0);
     244          _matchingWeight += (*_pPot)[x] - dist_y;
     245
     246          min_dist = dist_y;
     247          break;
     248        }
     249        else
     250        {
     251          // y is not free, so there *must* be only one arc pointing toward X
     252          const Edge e = (*_pMatchingMap)[y];
     253          assert(_graph.blueNode(e) == y);
     254
     255          const RedNode x2 = _graph.redNode(e);
     256          pred[x2] = _graph.direct(e, y);
     257          visitedX.push_back(x2);
     258          dist[x2] = dist_y; // matched edges have a reduced weight of 0
     259
     260          if (dist_y + (*_pPot)[x2] < x_min_dist)
     261          {
     262            x_min = x2;
     263            x_min_dist = dist_y + (*_pPot)[x2];
     264
     265            // we have a better criterion now
     266            if (UB > x_min_dist)
     267              UB = x_min_dist;
     268          }
     269
     270          for (IncEdgeIt e2(_graph, x2); e2 != INVALID; ++e2)
     271          {
     272            if (static_cast<const Edge>(e2) == e) continue;
     273
     274            const BlueNode y2 = _graph.blueNode(e2);
     275
     276            Value dist_y2 = dist_y + (*_pPot)[x2] + (*_pPot)[y2] - _weight[e2];
     277
     278            if (dist_y2 >= UB)
     279              continue;
     280
     281            if (isFree(y2))
     282              UB = dist_y2;
     283
     284            if (heap.state(y2) == Heap::PRE_HEAP)
     285            {
     286              dist[y2] = dist_y2;
     287              pred[y2] = _graph.direct(e2, x2);
     288              heap.push(y2, dist_y2);
     289            }
     290            else if (dist_y2 < dist[y2])
     291            {
     292              dist[y2] = dist_y2;
     293              pred[y2] = _graph.direct(e2, x2);
     294              heap.decrease(y2, dist_y2);
     295            }
     296          }
     297        }
     298      }
     299
     300      // restore primal and dual complementary slackness
     301      for (typename RedNodeList::const_iterator itX = visitedX.begin();
     302        itX != visitedX.end(); itX++)
     303      {
     304        const RedNode& x = *itX;
     305        assert(min_dist - dist[x] >= 0);
     306        (*_pPot)[x] -= min_dist - dist[x];
     307        assert((*_pPot)[x] >= 0);
     308      }
     309
     310      for (typename BlueNodeList::const_iterator itY = visitedY.begin();
     311        itY != visitedY.end(); itY++)
     312      {
     313        const BlueNode& y = *itY;
     314        assert(min_dist - dist[y] >= 0);
     315        (*_pPot)[y] += min_dist - dist[y];
     316        assert((*_pPot)[y] >= 0);
     317      }
     318    }
     319
     320  public:
     321    /// \brief Constructor
     322    ///
     323    /// Constructor.
     324    ///
     325    /// \param graph is the input graph
     326    /// \param weight are the edge weights
     327    MaxWeightedBipartiteMatching(const BpGraph& graph, const WeightMap& weight)
     328      : _graph(graph)
     329      , _weight(weight)
     330      , _pPot(NULL)
     331      , _pMatchingMap(NULL)
     332      , _matchingWeight(0)
     333      , _matchingSize(0)
     334    {
     335    }
     336
     337    ~MaxWeightedBipartiteMatching()
     338    {
     339      destroyStructures();
     340    }
     341
     342    /// \brief Initialize the algorithm
     343    ///
     344    /// This function initializes the algorithm.
     345    ///
     346    /// \param greedy indicates whether a nonempty initial matching
     347    /// should be used; this might be faster in some cases.
     348    void init(bool greedy = true)
     349    {
     350      destroyStructures();
     351      createStructures();
     352      _matchingWeight = 0;
     353      _matchingSize = 0;
     354
     355      // pot[x] is set to maximum incident edge weight
     356      for (RedIt x(_graph); x != INVALID; ++x)
     357      {
     358        Value max_weight = 0;
     359        Edge e_max = INVALID;
     360        for (IncEdgeIt e(_graph, x); e != INVALID; ++e)
     361        {
     362          // pot[y] = 0 for all y \in Y
     363          assert((*_pPot)[_graph.blueNode(e)] ==  0);
     364
     365          if (_weight[e] > max_weight)
     366          {
     367            max_weight = _weight[e];
     368            e_max = e;
     369          }
     370        }
     371
     372        if (e_max != INVALID)
     373        {
     374          _pPot->set(x, max_weight);
     375
     376          const Node y = _graph.blueNode(e_max);
     377          if (greedy && isFree(y))
     378          {
     379            _matchingWeight += max_weight;
     380            _matchingSize++;
     381            _pMatchingMap->set(x, e_max);
     382            _pMatchingMap->set(y, e_max);
     383          }
     384        }
     385      }
     386    }
     387
     388    /// \brief Start the algorithm
     389    ///
     390    /// This function starts the algorithm.
     391    ///
     392    /// \pre \ref init() must have been called before using this function.
     393    void start()
     394    {
     395      DistMap dist(_graph, 0);
     396      PredMap pred(_graph, INVALID);
     397
     398      for (RedIt x(_graph); x != INVALID; ++x)
     399      {
     400        if (isFree(x))
     401          augment(x, dist, pred);
     402      }
     403    }
     404
     405    /// \brief Run the algorithm.
     406    ///
     407    /// This method runs the \c %MaxWeightedBipartiteMatching algorithm.
     408    ///
     409    /// \param greedy indicates whether a nonempty initial matching
     410    /// should be used; this might be faster in some cases.
     411    ///
     412    /// \note mwbm.run() is just a shortcut of the following code.
     413    /// \code
     414    ///   mwbm.init();
     415    ///   mwbm.start();
     416    /// \endcode
     417    void run(bool greedy = true)
     418    {
     419      init(greedy);
     420      start();
     421    }
     422
     423    /// \brief Check whether the solution is optimal
     424    ///
     425    /// Check using the dual solution whether the primal solution is optimal.
     426    ///
     427    /// \return \c true if the solution is optimal.
     428    bool checkOptimality() const
     429    {
     430      assert(_pMatchingMap && _pPot);
     431
     432      /*
     433       * Primal:
     434       *   max  \sum_{i,j} c_{ij} x_{ij}
     435       *   s.t. \sum_i x_{ij} <= 1
     436       *        \sum_j x_{ij} <= 1
     437       *               x_{ij} >= 0
     438       *
     439       * Dual:
     440       *   min  \sum_j p_j + \sum_i r_i
     441       *   s.t. p_j + r_i >= c_{ij}
     442       *              p_j >= 0
     443       *              r_i >= 0
     444       *
     445       * Solution is optimal iff:
     446       * - Primal complementary slackness:
     447       *   - x_{ij} = 1  =>  p_j + r_i = c_{ij}
     448       * - Dual complementary slackness:
     449       *   - p_j != 0    =>  \sum_i x_{ij} = 1
     450       *   - r_i != 0    =>  \sum_j x_{ij} = 1
     451       */
     452
     453      // check primal solution
     454      for (NodeIt n(_graph); n != INVALID; ++n)
     455      {
     456        const Edge e = (*_pMatchingMap)[n];
     457
     458        if (e != INVALID)
     459        {
     460          const Node u = _graph.u(e);
     461          const Node v = _graph.v(e);
     462
     463          if (n != u && n != v)
     464            return false; // e must be incident to n
     465          if ((*_pMatchingMap)[u] != (*_pMatchingMap)[v])
     466            return false; // primal feasibility
     467          if ((*_pPot)[u] + (*_pPot)[v] != _weight[e])
     468            return false; // primal complementary slackness
     469        }
     470      }
     471
     472      // check dual solution
     473      for (NodeIt n(_graph); n != INVALID; ++n)
     474      {
     475        const Value pot_n = (*_pPot)[n];
     476        if (pot_n < 0)
     477          return false; // dual feasibility
     478        if ((*_pMatchingMap)[n] == INVALID && pot_n != 0)
     479          return false; // dual complementary slackness
     480      }
     481      for (EdgeIt e(_graph); e != INVALID; ++e)
     482      {
     483        if ((*_pPot)[_graph.u(e)] + (*_pPot)[_graph.v(e)] < _weight[e])
     484          return false; // dual feasibility
     485      }
     486
     487      return true;
     488    }
     489
     490    /// \brief Return the dual value of the given node
     491    ///
     492    /// This function returns the potential of the given node
     493    ///
     494    /// \pre init() must have been called before using this function
     495    const Value pot(const Node& n) const
     496    {
     497      assert(_pPot);
     498      return (*_pPot)[n];
     499    }
     500
     501    /// \brief Return a const reference to the matching map.
     502    ///
     503    /// This function returns a const reference to a node map that stores
     504    /// the matching edge incident to each node.
     505    ///
     506    /// \pre init() must have been called before using this function.
     507    const MatchingMap& matchingMap() const
     508    {
     509      assert(_pMatchingMap);
     510      return *_pMatchingMap;
     511    }
     512
     513    /// \brief Return the weight of the matching.
     514    ///
     515    /// This function returns the weight of the found matching.
     516    ///
     517    /// \pre init() must have been called before using this function.
     518    Value matchingWeight() const
     519    {
     520      return _matchingWeight;
     521    }
     522
     523    /// \brief Return the number of edges in the matching.
     524    ///
     525    /// This function returns the number of edges in the matching.
     526    int matchingSize() const
     527    {
     528      return _matchingSize;
     529    }
     530
     531    /// \brief Return \c true if the given edge is in the matching.
     532    ///
     533    /// This function returns \c true if the given edge is in the found
     534    /// matching.
     535    ///
     536    /// \pre init() must have been been called before using this function.
     537    bool matching(const Edge& e) const
     538    {
     539      assert(_pMatchingMap);
     540      return e != INVALID && (*_pMatchingMap)[_graph.u(e)] != INVALID;
     541    }
     542
     543    /// \brief Return the matching edge incident to the given node.
     544    ///
     545    /// This function returns the matching edge incident to the
     546    /// given node in the found matching or \c INVALID if the node is
     547    /// not covered by the matching.
     548    ///
     549    /// \pre init() must have been been called before using this function.
     550    Edge matching(const Node& n) const
     551    {
     552      assert(_pMatchingMap);
     553      return (*_pMatchingMap)[n];
     554    }
     555
     556    /// \brief Return the mate of the given node.
     557    ///
     558    /// This function returns the mate of the given node in the found
     559    /// matching or \c INVALID if the node is not covered by the matching.
     560    ///
     561    /// \pre init() must have been been called before using this function.
     562    Node mate(const Node& n) const
     563    {
     564      assert(_pMatchingMap);
     565
     566      const Edge e = (*_pMatchingMap)[n];
     567
     568      if (e == INVALID)
     569        return INVALID;
     570      else
     571        return _graph.oppositeNode(n, e);
     572    }
     573  };
     574
     575  /// \ingroup matching
     576  ///
     577  /// \brief Maximum weight matching in (dense) bipartite graphs
     578  ///
     579  /// This class provides an implementation of the classical Hungarian
     580  /// algorithm for finding a maximum weight matching in an undirected
     581  /// bipartite graph. This algorithm follows the primal-dual schema.
     582  /// The time complexity is \f$O(n^3)\f$. In case the bipartite graph is
     583  /// sparse, it is better to use \ref MaxWeightedBipartiteMatching, which
     584  /// has a time complexity of \f$O(n^2 \log n)\f$ for sparse graphs.
     585  ///
     586  /// \tparam BGR The bipartite graph type the algorithm runs on.
     587  /// \tparam WM The type edge weight map. The default type is
     588  /// \ref concepts::Graph::EdgeMap "BGR::EdgeMap<int>".
     589#ifdef DOXYGEN
     590  template <typename BGR, typename WM>
     591#else
     592  template <typename BGR,
     593            typename WM = typename BGR::template EdgeMap<int> >
     594#endif
     595  class MaxWeightedDenseBipartiteMatching
     596  {
     597  public:
     598    /// The graph type of the algorithm
     599    typedef BGR BpGraph;
     600    /// The type of the edge weight map
     601    typedef WM WeightMap;
     602    /// The value type of the edge weights
     603    typedef typename WeightMap::Value Value;
     604    /// The type of the matching map
     605    typedef typename BpGraph::
     606      template NodeMap<typename BpGraph::Edge> MatchingMap;
     607
     608  private:
     609    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     610
     611    typedef typename BpGraph::template NodeMap<int> IdMap;
     612    typedef std::vector<int> MateVector;
     613    typedef std::vector<Value> WeightVector;
     614    typedef std::vector<bool> BoolVector;
     615
     616    class BpEdgeT
     617    {
     618    private:
     619      Value _weight;
     620      Edge _edge;
     621
     622    public:
     623      BpEdgeT()
     624        : _weight(0)
     625        , _edge(INVALID)
     626      {
     627      }
     628
     629      void setWeight(Value weight)
     630      {
     631        _weight = weight;
     632      }
     633
     634      Value getWeight() const
     635      {
     636        return _weight;
     637      }
     638
     639      void setEdge(const Edge& edge)
     640      {
     641        _edge = edge;
     642      }
     643
     644      const Edge& getEdge() const
     645      {
     646        return _edge;
     647      }
     648    };
     649
     650    typedef std::vector<std::vector<BpEdgeT> > AdjacencyMatrixType;
     651
     652    const BpGraph& _graph;
     653    const WeightMap& _weight;
     654    IdMap _idMap;
     655    MatchingMap _matchingMap;
     656
     657    AdjacencyMatrixType _adjacencyMatrix;
     658    WeightVector _labelMapX;
     659    WeightVector _labelMapY;
     660    MateVector _mateMapX;
     661    MateVector _mateMapY;
     662    int _nX;
     663    int _nY;
     664    int _matchingSize;
     665    Value _matchingWeight;
     666
     667    static const Value _minValue;
     668    static const Value _maxValue;
     669
     670    void buildMatchingMap()
     671    {
     672      _matchingWeight = 0;
     673      _matchingSize = 0;
     674
     675      for (int x = 0; x < _nX; x++)
     676      {
     677        assert(_mateMapX[x] != -1);
     678        int y = _mateMapX[x];
     679
     680        const Edge& e = _adjacencyMatrix[x][y].getEdge();
     681        if (e != INVALID)
     682        {
     683          // only edges that where present
     684          // in the original graph count in the matching
     685          _matchingMap[_graph.u(e)] = _matchingMap[_graph.v(e)] = e;
     686          _matchingSize++;
     687          _matchingWeight += _weight[e];
     688        }
     689      }
     690    }
     691
     692    void updateSlacks(WeightVector& slack, int x)
     693    {
     694      Value lx = _labelMapX[x];
     695      for (int y = 0; y < _nY; y++)
     696      {
     697        // slack[y] = min_{x \in S} [l(x) + l(y) - w(x, y)]
     698        Value val = lx + _labelMapY[y] - _adjacencyMatrix[x][y].getWeight();
     699        if (slack[y] > val)
     700          slack[y] = val;
     701      }
     702    }
     703
     704    void updateLabels(const BoolVector& setS,
     705                      const BoolVector& setT, WeightVector& slack)
     706    {
     707      // recall that slack[y] = min_{x \in S} [l(x) + l(y) - w(x,y)]
     708
     709      // delta = min_{y \not \in T} (slack[y])
     710      Value delta = _maxValue;
     711      for (int y = 0; y < _nY; y++)
     712      {
     713        if (!setT[y] && slack[y] < delta)
     714          delta = slack[y];
     715      }
     716
     717      // update labels in X
     718      for (int x = 0; x < _nX; x++)
     719      {
     720        if (setS[x])
     721          _labelMapX[x] -= delta;
     722      }
     723
     724      // update labels in Y
     725      for (int y = 0; y < _nY; y++)
     726      {
     727        if (setT[y])
     728          _labelMapY[y] += delta;
     729        else
     730        {
     731          // update slacks
     732          // remember that l(x) + l(y) hasn't changed for x \in S and y \in T
     733          // the only thing that has changed is"
     734          //   l(x) + l(y) for x \in S and y \not \in T
     735          slack[y] -= delta;
     736        }
     737      }
     738    }
     739
     740  public:
     741    /// \brief Constructor
     742    ///
     743    /// Constructor.
     744    ///
     745    /// \param graph is the input graph
     746    /// \param weight are the edge weights
     747    MaxWeightedDenseBipartiteMatching(const BpGraph& graph,
     748                                      const WeightMap& weight)
     749      : _graph(graph)
     750      , _weight(weight)
     751      , _idMap(graph, -1)
     752      , _matchingMap(graph, INVALID)
     753      , _adjacencyMatrix()
     754      , _labelMapX()
     755      , _labelMapY()
     756      , _mateMapX()
     757      , _mateMapY()
     758      , _nX(0)
     759      , _nY(0)
     760      , _matchingSize(0)
     761      , _matchingWeight(0)
     762    {
     763
     764    }
     765
     766    /// \brief Initialize the algorithm
     767    ///
     768    /// This function initializes the algorithm.
     769    void init()
     770    {
     771      // construct _idMap
     772      int id_x = 0;
     773      for (RedIt x(_graph); x != INVALID; ++x, ++id_x)
     774      {
     775        _idMap[x] = id_x;
     776      }
     777      _nX = id_x;
     778
     779      int id_y = 0;
     780      for (BlueIt y(_graph); y != INVALID; ++y, ++id_y)
     781      {
     782        _idMap[y] = id_y;
     783      }
     784      _nY = id_y;
     785
     786      assert(_nX <= _nY);
     787
     788      // init matching is empty
     789      _mateMapX = MateVector(_nX, -1);
     790      _mateMapY = MateVector(_nY, -1);
     791
     792      // labels of nodes in X are initialized to 0,
     793      // these will be updated during initAdjacencyMatrix()
     794      _labelMapX = WeightVector(_nX, 0);
     795
     796      // labels of nodes in Y are initialized to 0,
     797      // these won't be updated during initAdjacencyMatrix()
     798      _labelMapY = WeightVector(_nY, 0);
     799
     800      // adjacency matrix has dimensions |X| * |Y|,
     801      // every entry in this matrix is initialized to (0, INVALID)
     802      _adjacencyMatrix = AdjacencyMatrixType(_nX,
     803        std::vector<BpEdgeT>(_nY, BpEdgeT()));
     804
     805      _matchingWeight = 0;
     806      _matchingSize = 0;
     807
     808      for (RedIt x(_graph); x != INVALID; ++x)
     809      {
     810        id_x = _idMap[x];
     811        for (IncEdgeIt e(_graph, x); e != INVALID; ++e)
     812        {
     813          Node y = _graph.blueNode(e);
     814          id_y = _idMap[y];
     815
     816          Value w = _weight[e];
     817
     818          BpEdgeT& item = _adjacencyMatrix[id_x][id_y];
     819          item.setEdge(e);
     820          item.setWeight(w);
     821
     822          // label of a node x in X is initialized to maximum weight
     823          // of edges incident to x
     824          if (w > _labelMapX[id_x])
     825            _labelMapX[id_x] = w;
     826        }
     827      }
     828    }
     829
     830    /// \brief Run the algorithm.
     831    ///
     832    /// This method runs the \c %MaxWeightedDenseBipartiteMatching algorithm.
     833    ///
     834    /// \note mwdbm.run() is just a shortcut of the following code.
     835    /// \code
     836    ///   mwdbm.init()
     837    ///   mwdbm.start();
     838    /// \endcode
     839    void run()
     840    {
     841      init();
     842      start();
     843    }
     844
     845    /// \brief Start the algorithm
     846    ///
     847    /// This function starts the algorithm.
     848    ///
     849    /// \pre \ref init() must have been called before using this function.
     850    void start()
     851    {
     852      // maps y in Y to x in X by which it was discovered
     853      MateVector discoveredY(_nY, -1);
     854
     855      // pick a root
     856      for (int r = 0; r < _nX; )
     857      {
     858        assert(_mateMapX[r] == -1);
     859
     860        // clear slack map, i.e. set all slacks to +INF
     861        WeightVector slack(_nY, _maxValue);
     862
     863        // initially T = {}
     864        BoolVector setT(_nY, false);
     865
     866        // initially S = {r}
     867        BoolVector setS(_nX, false);
     868        setS[r] = true;
     869
     870        std::queue<int> queue;
     871        queue.push(r);
     872
     873        updateSlacks(slack, r);
     874
     875        bool augmented = false;
     876        while (!queue.empty() && !augmented)
     877        {
     878          int x = queue.front();
     879          queue.pop();
     880
     881          for (int y = 0; y < _nY; y++)
     882          {
     883            if (!setT[y] &&
     884              _labelMapX[x] + _labelMapY[y] ==
     885              _adjacencyMatrix[x][y].getWeight())
     886            {
     887              // y was (first) discovered by x
     888              discoveredY[y] = x;
     889
     890              if (_mateMapY[y] != -1) // y is matched, extend alternating tree
     891              {
     892                int z = _mateMapY[y];
     893
     894                // add z to queue if not in S
     895                if (!setS[z])
     896                {
     897                  setS[z] = true;
     898                  queue.push(z);
     899                  updateSlacks(slack, z);
     900                }
     901
     902                setT[y] = true;
     903              }
     904              else // y is free, we have an augmenting path between r and y
     905              {
     906                int cx, ty, cy = y;
     907                do {
     908                  cx = discoveredY[cy];
     909                  ty = _mateMapX[cx];
     910
     911                  _mateMapX[cx] = cy;
     912                  _mateMapY[cy] = cx;
     913
     914                  cy = ty;
     915                } while (cx != r);
     916
     917                // we found an augmenting path,
     918                // start a new iteration of the first for loop
     919                augmented = true;
     920                break; // break for y
     921              }
     922            }
     923          } // y \not in T such that (r,y) in E_l
     924        } // queue
     925
     926        if (!augmented)
     927          updateLabels(setS, setT, slack);
     928        else
     929          r++;
     930      }
     931
     932      buildMatchingMap();
     933    }
     934
     935    /// \brief Return the dual value of the given node
     936    ///
     937    /// This function returns the potential of the given node
     938    ///
     939    /// \pre init() must have been called before using this function
     940    const Value pot(const Node& n) const
     941    {
     942      if (_graph.red(n))
     943        return _labelMapX[_idMap[n]];
     944      else
     945        return _labelMapY[_idMap[n]];
     946    }
     947
     948    /// \brief Return the weight of the matching.
     949    ///
     950    /// This function returns the weight of the found matching.
     951    ///
     952    /// \pre init() must have been called before using this function.
     953    Value matchingWeight() const
     954    {
     955      return _matchingWeight;
     956    }
     957
     958    /// \brief Return the number of edges in the matching.
     959    ///
     960    /// This function returns the number of edges in the matching.
     961    int matchingSize() const
     962    {
     963      return _matchingSize;
     964    }
     965
     966    /// \brief Return \c true if the given edge is in the matching.
     967    ///
     968    /// This function returns \c true if the given edge is in the found
     969    /// matching.
     970    ///
     971    /// \pre init() must have been been called before using this function.
     972    bool matching(const Edge& e) const
     973    {
     974      return _matchingMap[_graph.u(e)] != INVALID;
     975    }
     976
     977    /// \brief Return the matching edge incident to the given node.
     978    ///
     979    /// This function returns the matching edge incident to the
     980    /// given node in the found matching or \c INVALID if the node is
     981    /// not covered by the matching.
     982    ///
     983    /// \pre init() must have been been called before using this function.
     984    Edge matching(const Node& n) const
     985    {
     986      return _matchingMap[n];
     987    }
     988
     989    /// \brief Return the mate of the given node.
     990    ///
     991    /// This function returns the mate of the given node in the found
     992    /// matching or \c INVALID if the node is not covered by the matching.
     993    ///
     994    /// \pre init() must have been been called before using this function.
     995    Node mate(const Node& n) const
     996    {
     997      return _graph.oppositeNode(n, _matchingMap[n]);
     998    }
     999
     1000    /// \brief Return a const reference to the matching map.
     1001    ///
     1002    /// This function returns a const reference to a node map that stores
     1003    /// the matching edge incident to each node.
     1004    ///
     1005    /// \pre init() must have been called before using this function.
     1006    const MatchingMap& matchingMap() const
     1007    {
     1008      return _matchingMap;
     1009    }
     1010
     1011    /// \brief Checks whether the solution is optimal
     1012    ///
     1013    /// Checks using the dual solution whether the primal solution is optimal.
     1014    ///
     1015    /// \return \c true if the solution is optimal.
     1016    bool checkOptimality() const
     1017    {
     1018      for (RedIt x(_graph); x != INVALID; ++x)
     1019      {
     1020        if (_labelMapX[_idMap[x]] < 0)
     1021          return false; // feasibility of dual solution
     1022
     1023        const Edge e = _matchingMap[x];
     1024
     1025        if (_mateMapX[_idMap[x]] == -1)
     1026          return false; // all nodes in X must be matched in the complete graph
     1027        else if (e != INVALID)
     1028        {
     1029          const Node y = _graph.blueNode(e);
     1030
     1031          if (_matchingMap[y] != e)
     1032            return false; // if x is matched via e then so must y
     1033          if (_labelMapX[_idMap[x]] + _labelMapY[_idMap[y]] !=
     1034              _adjacencyMatrix[_idMap[x]][_idMap[y]].getWeight())
     1035            return false; // primal complementary slackness
     1036        }
     1037      }
     1038
     1039      for (BlueIt y(_graph); y != INVALID; ++y)
     1040      {
     1041        if (_labelMapY[_idMap[y]] < 0)
     1042          return false; // feasibility of dual solution
     1043
     1044        if (_matchingMap[y] == INVALID && _labelMapY[_idMap[y]] != 0)
     1045          return false; // dual complementary slackness
     1046      }
     1047
     1048      for (EdgeIt e(_graph); e != INVALID; ++e)
     1049      {
     1050        // feasibility of dual solution
     1051        if (_labelMapX[_idMap[_graph.redNode(e)]] +
     1052            _labelMapY[_idMap[_graph.blueNode(e)]] < _weight[e])
     1053            return false;
     1054      }
     1055
     1056      return true;
     1057    }
     1058  };
     1059
     1060  template<typename BGR, typename WM>
     1061  const typename WM::Value MaxWeightedDenseBipartiteMatching<BGR, WM>::
     1062    _maxValue = std::numeric_limits<typename WM::Value>::max();
     1063}
     1064
     1065#endif //BP_MATCHING_H
  • test/CMakeLists.txt

    diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
    a b  
    1111  adaptors_test
    1212  bellman_ford_test
    1313  bfs_test
     14  bp_matching_test
    1415  bpgraph_test
    1516  circulation_test
    1617  connectivity_test