COIN-OR::LEMON - Graph Library

Ticket #168: bp_matching.h

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