COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/bipartite_matching.h @ 2462:7a096a6bf53a

Last change on this file since 2462:7a096a6bf53a was 2462:7a096a6bf53a, checked in by Balazs Dezso, 17 years ago

Common interface for bipartite matchings
Some useful query function for push-relabel based matching

The naming should be rethink for these classes
for example: pr-ap prefix for push-relabel and augmenting path
algorithms

File size: 52.7 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2007
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_BIPARTITE_MATCHING
20#define LEMON_BIPARTITE_MATCHING
21
22#include <functional>
23
24#include <lemon/bin_heap.h>
25#include <lemon/maps.h>
26
27#include <iostream>
28
29///\ingroup matching
30///\file
31///\brief Maximum matching algorithms in bipartite graphs.
32///
33///\note The pr_bipartite_matching.h file also contains algorithms to
34///solve maximum cardinality bipartite matching problems.
35
36namespace lemon {
37
38  /// \ingroup matching
39  ///
40  /// \brief Bipartite Max Cardinality Matching algorithm
41  ///
42  /// Bipartite Max Cardinality Matching algorithm. This class implements
43  /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
44  /// complexity.
45  ///
46  /// \note In several cases the push-relabel based algorithms have
47  /// better runtime performance than the augmenting path based ones.
48  ///
49  /// \see PrBipartiteMatching
50  template <typename BpUGraph>
51  class MaxBipartiteMatching {
52  protected:
53
54    typedef BpUGraph Graph;
55
56    typedef typename Graph::Node Node;
57    typedef typename Graph::ANodeIt ANodeIt;
58    typedef typename Graph::BNodeIt BNodeIt;
59    typedef typename Graph::UEdge UEdge;
60    typedef typename Graph::UEdgeIt UEdgeIt;
61    typedef typename Graph::IncEdgeIt IncEdgeIt;
62
63    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
64    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
65
66
67  public:
68
69    /// \brief Constructor.
70    ///
71    /// Constructor of the algorithm.
72    MaxBipartiteMatching(const BpUGraph& _graph)
73      : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}
74
75    /// \name Execution control
76    /// The simplest way to execute the algorithm is to use
77    /// one of the member functions called \c run().
78    /// \n
79    /// If you need more control on the execution,
80    /// first you must call \ref init() or one alternative for it.
81    /// Finally \ref start() will perform the matching computation or
82    /// with step-by-step execution you can augment the solution.
83
84    /// @{
85
86    /// \brief Initalize the data structures.
87    ///
88    /// It initalizes the data structures and creates an empty matching.
89    void init() {
90      for (ANodeIt it(*graph); it != INVALID; ++it) {
91        anode_matching[it] = INVALID;
92      }
93      for (BNodeIt it(*graph); it != INVALID; ++it) {
94        bnode_matching[it] = INVALID;
95      }
96      matching_size = 0;
97    }
98
99    /// \brief Initalize the data structures.
100    ///
101    /// It initalizes the data structures and creates a greedy
102    /// matching.  From this matching sometimes it is faster to get
103    /// the matching than from the initial empty matching.
104    void greedyInit() {
105      matching_size = 0;
106      for (BNodeIt it(*graph); it != INVALID; ++it) {
107        bnode_matching[it] = INVALID;
108      }
109      for (ANodeIt it(*graph); it != INVALID; ++it) {
110        anode_matching[it] = INVALID;
111        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
112          if (bnode_matching[graph->bNode(jt)] == INVALID) {
113            anode_matching[it] = jt;
114            bnode_matching[graph->bNode(jt)] = jt;
115            ++matching_size;
116            break;
117          }
118        }
119      }
120    }
121
122    /// \brief Initalize the data structures with an initial matching.
123    ///
124    /// It initalizes the data structures with an initial matching.
125    template <typename MatchingMap>
126    void matchingInit(const MatchingMap& mm) {
127      for (ANodeIt it(*graph); it != INVALID; ++it) {
128        anode_matching[it] = INVALID;
129      }
130      for (BNodeIt it(*graph); it != INVALID; ++it) {
131        bnode_matching[it] = INVALID;
132      }
133      matching_size = 0;
134      for (UEdgeIt it(*graph); it != INVALID; ++it) {
135        if (mm[it]) {
136          ++matching_size;
137          anode_matching[graph->aNode(it)] = it;
138          bnode_matching[graph->bNode(it)] = it;
139        }
140      }
141    }
142
143    /// \brief Initalize the data structures with an initial matching.
144    ///
145    /// It initalizes the data structures with an initial matching.
146    /// \return %True when the given map contains really a matching.
147    template <typename MatchingMap>
148    void checkedMatchingInit(const MatchingMap& mm) {
149      for (ANodeIt it(*graph); it != INVALID; ++it) {
150        anode_matching[it] = INVALID;
151      }
152      for (BNodeIt it(*graph); it != INVALID; ++it) {
153        bnode_matching[it] = INVALID;
154      }
155      matching_size = 0;
156      for (UEdgeIt it(*graph); it != INVALID; ++it) {
157        if (mm[it]) {
158          ++matching_size;
159          if (anode_matching[graph->aNode(it)] != INVALID) {
160            return false;
161          }
162          anode_matching[graph->aNode(it)] = it;
163          if (bnode_matching[graph->aNode(it)] != INVALID) {
164            return false;
165          }
166          bnode_matching[graph->bNode(it)] = it;
167        }
168      }
169      return false;
170    }
171
172    /// \brief An augmenting phase of the Hopcroft-Karp algorithm
173    ///
174    /// It runs an augmenting phase of the Hopcroft-Karp
175    /// algorithm. This phase finds maximum count of edge disjoint
176    /// augmenting paths and augments on these paths. The algorithm
177    /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
178    /// \f$ O(e) \f$ long.
179    bool augment() {
180
181      typename Graph::template ANodeMap<bool> areached(*graph, false);
182      typename Graph::template BNodeMap<bool> breached(*graph, false);
183
184      typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
185
186      std::vector<Node> queue, bqueue;
187      for (ANodeIt it(*graph); it != INVALID; ++it) {
188        if (anode_matching[it] == INVALID) {
189          queue.push_back(it);
190          areached[it] = true;
191        }
192      }
193
194      bool success = false;
195
196      while (!success && !queue.empty()) {
197        std::vector<Node> newqueue;
198        for (int i = 0; i < int(queue.size()); ++i) {
199          Node anode = queue[i];
200          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
201            Node bnode = graph->bNode(jt);
202            if (breached[bnode]) continue;
203            breached[bnode] = true;
204            bpred[bnode] = jt;
205            if (bnode_matching[bnode] == INVALID) {
206              bqueue.push_back(bnode);
207              success = true;
208            } else {           
209              Node newanode = graph->aNode(bnode_matching[bnode]);
210              if (!areached[newanode]) {
211                areached[newanode] = true;
212                newqueue.push_back(newanode);
213              }
214            }
215          }
216        }
217        queue.swap(newqueue);
218      }
219
220      if (success) {
221
222        typename Graph::template ANodeMap<bool> aused(*graph, false);
223       
224        for (int i = 0; i < int(bqueue.size()); ++i) {
225          Node bnode = bqueue[i];
226
227          bool used = false;
228
229          while (bnode != INVALID) {
230            UEdge uedge = bpred[bnode];
231            Node anode = graph->aNode(uedge);
232           
233            if (aused[anode]) {
234              used = true;
235              break;
236            }
237           
238            bnode = anode_matching[anode] != INVALID ?
239              graph->bNode(anode_matching[anode]) : INVALID;
240           
241          }
242         
243          if (used) continue;
244
245          bnode = bqueue[i];
246          while (bnode != INVALID) {
247            UEdge uedge = bpred[bnode];
248            Node anode = graph->aNode(uedge);
249           
250            bnode_matching[bnode] = uedge;
251
252            bnode = anode_matching[anode] != INVALID ?
253              graph->bNode(anode_matching[anode]) : INVALID;
254           
255            anode_matching[anode] = uedge;
256
257            aused[anode] = true;
258          }
259          ++matching_size;
260
261        }
262      }
263      return success;
264    }
265
266    /// \brief An augmenting phase of the Ford-Fulkerson algorithm
267    ///
268    /// It runs an augmenting phase of the Ford-Fulkerson
269    /// algorithm. This phase finds only one augmenting path and
270    /// augments only on this paths. The algorithm consists at most
271    /// of \f$ O(n) \f$ simple phase and one phase is at most
272    /// \f$ O(e) \f$ long.
273    bool simpleAugment() {
274
275      typename Graph::template ANodeMap<bool> areached(*graph, false);
276      typename Graph::template BNodeMap<bool> breached(*graph, false);
277
278      typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
279
280      std::vector<Node> queue;
281      for (ANodeIt it(*graph); it != INVALID; ++it) {
282        if (anode_matching[it] == INVALID) {
283          queue.push_back(it);
284          areached[it] = true;
285        }
286      }
287
288      while (!queue.empty()) {
289        std::vector<Node> newqueue;
290        for (int i = 0; i < int(queue.size()); ++i) {
291          Node anode = queue[i];
292          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
293            Node bnode = graph->bNode(jt);
294            if (breached[bnode]) continue;
295            breached[bnode] = true;
296            bpred[bnode] = jt;
297            if (bnode_matching[bnode] == INVALID) {
298              while (bnode != INVALID) {
299                UEdge uedge = bpred[bnode];
300                anode = graph->aNode(uedge);
301               
302                bnode_matching[bnode] = uedge;
303               
304                bnode = anode_matching[anode] != INVALID ?
305                  graph->bNode(anode_matching[anode]) : INVALID;
306               
307                anode_matching[anode] = uedge;
308               
309              }
310              ++matching_size;
311              return true;
312            } else {           
313              Node newanode = graph->aNode(bnode_matching[bnode]);
314              if (!areached[newanode]) {
315                areached[newanode] = true;
316                newqueue.push_back(newanode);
317              }
318            }
319          }
320        }
321        queue.swap(newqueue);
322      }
323     
324      return false;
325    }
326
327    /// \brief Starts the algorithm.
328    ///
329    /// Starts the algorithm. It runs augmenting phases until the optimal
330    /// solution reached.
331    void start() {
332      while (augment()) {}
333    }
334
335    /// \brief Runs the algorithm.
336    ///
337    /// It just initalize the algorithm and then start it.
338    void run() {
339      greedyInit();
340      start();
341    }
342
343    /// @}
344
345    /// \name Query Functions
346    /// The result of the %Matching algorithm can be obtained using these
347    /// functions.\n
348    /// Before the use of these functions,
349    /// either run() or start() must be called.
350   
351    ///@{
352
353    /// \brief Set true all matching uedge in the map.
354    ///
355    /// Set true all matching uedge in the map. It does not change the
356    /// value mapped to the other uedges.
357    /// \return The number of the matching edges.
358    template <typename MatchingMap>
359    int quickMatching(MatchingMap& mm) const {
360      for (ANodeIt it(*graph); it != INVALID; ++it) {
361        if (anode_matching[it] != INVALID) {
362          mm[anode_matching[it]] = true;
363        }
364      }
365      return matching_size;
366    }
367
368    /// \brief Set true all matching uedge in the map and the others to false.
369    ///
370    /// Set true all matching uedge in the map and the others to false.
371    /// \return The number of the matching edges.
372    template <typename MatchingMap>
373    int matching(MatchingMap& mm) const {
374      for (UEdgeIt it(*graph); it != INVALID; ++it) {
375        mm[it] = it == anode_matching[graph->aNode(it)];
376      }
377      return matching_size;
378    }
379
380
381    /// \brief Return true if the given uedge is in the matching.
382    ///
383    /// It returns true if the given uedge is in the matching.
384    bool matchingEdge(const UEdge& edge) const {
385      return anode_matching[graph->aNode(edge)] == edge;
386    }
387
388    /// \brief Returns the matching edge from the node.
389    ///
390    /// Returns the matching edge from the node. If there is not such
391    /// edge it gives back \c INVALID.
392    UEdge matchingEdge(const Node& node) const {
393      if (graph->aNode(node)) {
394        return anode_matching[node];
395      } else {
396        return bnode_matching[node];
397      }
398    }
399
400    /// \brief Gives back the number of the matching edges.
401    ///
402    /// Gives back the number of the matching edges.
403    int matchingSize() const {
404      return matching_size;
405    }
406
407    /// \brief Returns a minimum covering of the nodes.
408    ///
409    /// The minimum covering set problem is the dual solution of the
410    /// maximum bipartite matching. It provides a solution for this
411    /// problem what is proof of the optimality of the matching.
412    /// \return The size of the cover set.
413    template <typename CoverMap>
414    int coverSet(CoverMap& covering) const {
415
416      typename Graph::template ANodeMap<bool> areached(*graph, false);
417      typename Graph::template BNodeMap<bool> breached(*graph, false);
418     
419      std::vector<Node> queue;
420      for (ANodeIt it(*graph); it != INVALID; ++it) {
421        if (anode_matching[it] == INVALID) {
422          queue.push_back(it);
423        }
424      }
425
426      while (!queue.empty()) {
427        std::vector<Node> newqueue;
428        for (int i = 0; i < int(queue.size()); ++i) {
429          Node anode = queue[i];
430          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
431            Node bnode = graph->bNode(jt);
432            if (breached[bnode]) continue;
433            breached[bnode] = true;
434            if (bnode_matching[bnode] != INVALID) {
435              Node newanode = graph->aNode(bnode_matching[bnode]);
436              if (!areached[newanode]) {
437                areached[newanode] = true;
438                newqueue.push_back(newanode);
439              }
440            }
441          }
442        }
443        queue.swap(newqueue);
444      }
445
446      int size = 0;
447      for (ANodeIt it(*graph); it != INVALID; ++it) {
448        covering[it] = !areached[it] && anode_matching[it] != INVALID;
449        if (!areached[it] && anode_matching[it] != INVALID) {
450          ++size;
451        }
452      }
453      for (BNodeIt it(*graph); it != INVALID; ++it) {
454        covering[it] = breached[it];
455        if (breached[it]) {
456          ++size;
457        }
458      }
459      return size;
460    }
461
462    /// \brief Gives back a barrier on the A-nodes
463   
464    /// The barrier is s subset of the nodes on the same side of the
465    /// graph, which size minus its neighbours is exactly the
466    /// unmatched nodes on the A-side. 
467    /// \retval barrier A WriteMap on the ANodes with bool value.
468    template <typename BarrierMap>
469    void aBarrier(BarrierMap& barrier) const {
470
471      typename Graph::template ANodeMap<bool> areached(*graph, false);
472      typename Graph::template BNodeMap<bool> breached(*graph, false);
473     
474      std::vector<Node> queue;
475      for (ANodeIt it(*graph); it != INVALID; ++it) {
476        if (anode_matching[it] == INVALID) {
477          queue.push_back(it);
478        }
479      }
480
481      while (!queue.empty()) {
482        std::vector<Node> newqueue;
483        for (int i = 0; i < int(queue.size()); ++i) {
484          Node anode = queue[i];
485          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
486            Node bnode = graph->bNode(jt);
487            if (breached[bnode]) continue;
488            breached[bnode] = true;
489            if (bnode_matching[bnode] != INVALID) {
490              Node newanode = graph->aNode(bnode_matching[bnode]);
491              if (!areached[newanode]) {
492                areached[newanode] = true;
493                newqueue.push_back(newanode);
494              }
495            }
496          }
497        }
498        queue.swap(newqueue);
499      }
500
501      for (ANodeIt it(*graph); it != INVALID; ++it) {
502        barrier[it] = areached[it] || anode_matching[it] == INVALID;
503      }
504    }
505
506    /// \brief Gives back a barrier on the B-nodes
507   
508    /// The barrier is s subset of the nodes on the same side of the
509    /// graph, which size minus its neighbours is exactly the
510    /// unmatched nodes on the B-side. 
511    /// \retval barrier A WriteMap on the BNodes with bool value.
512    template <typename BarrierMap>
513    void bBarrier(BarrierMap& barrier) const {
514
515      typename Graph::template ANodeMap<bool> areached(*graph, false);
516      typename Graph::template BNodeMap<bool> breached(*graph, false);
517     
518      std::vector<Node> queue;
519      for (ANodeIt it(*graph); it != INVALID; ++it) {
520        if (anode_matching[it] == INVALID) {
521          queue.push_back(it);
522        }
523      }
524
525      while (!queue.empty()) {
526        std::vector<Node> newqueue;
527        for (int i = 0; i < int(queue.size()); ++i) {
528          Node anode = queue[i];
529          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
530            Node bnode = graph->bNode(jt);
531            if (breached[bnode]) continue;
532            breached[bnode] = true;
533            if (bnode_matching[bnode] != INVALID) {
534              Node newanode = graph->aNode(bnode_matching[bnode]);
535              if (!areached[newanode]) {
536                areached[newanode] = true;
537                newqueue.push_back(newanode);
538              }
539            }
540          }
541        }
542        queue.swap(newqueue);
543      }
544
545      for (BNodeIt it(*graph); it != INVALID; ++it) {
546        barrier[it] = !breached[it];
547      }
548    }
549
550    /// @}
551
552  private:
553
554    ANodeMatchingMap anode_matching;
555    BNodeMatchingMap bnode_matching;
556    const Graph *graph;
557
558    int matching_size;
559 
560  };
561
562  /// \ingroup matching
563  ///
564  /// \brief Maximum cardinality bipartite matching
565  ///
566  /// This function calculates the maximum cardinality matching
567  /// in a bipartite graph. It gives back the matching in an undirected
568  /// edge map.
569  ///
570  /// \param graph The bipartite graph.
571  /// \retval matching The undirected edge map which will be set to
572  /// the matching.
573  /// \return The size of the matching.
574  template <typename BpUGraph, typename MatchingMap>
575  int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) {
576    MaxBipartiteMatching<BpUGraph> bpmatching(graph);
577    bpmatching.run();
578    bpmatching.matching(matching);
579    return bpmatching.matchingSize();
580  }
581
582  /// \brief Default traits class for weighted bipartite matching algoritms.
583  ///
584  /// Default traits class for weighted bipartite matching algoritms.
585  /// \param _BpUGraph The bipartite undirected graph type.
586  /// \param _WeightMap Type of weight map.
587  template <typename _BpUGraph, typename _WeightMap>
588  struct WeightedBipartiteMatchingDefaultTraits {
589    /// \brief The type of the weight of the undirected edges.
590    typedef typename _WeightMap::Value Value;
591
592    /// The undirected bipartite graph type the algorithm runs on.
593    typedef _BpUGraph BpUGraph;
594
595    /// The map of the edges weights
596    typedef _WeightMap WeightMap;
597
598    /// \brief The cross reference type used by heap.
599    ///
600    /// The cross reference type used by heap.
601    /// Usually it is \c Graph::NodeMap<int>.
602    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
603
604    /// \brief Instantiates a HeapCrossRef.
605    ///
606    /// This function instantiates a \ref HeapCrossRef.
607    /// \param graph is the graph, to which we would like to define the
608    /// HeapCrossRef.
609    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
610      return new HeapCrossRef(graph);
611    }
612   
613    /// \brief The heap type used by weighted matching algorithms.
614    ///
615    /// The heap type used by weighted matching algorithms. It should
616    /// minimize the priorities and the heap's key type is the graph's
617    /// anode graph's node.
618    ///
619    /// \sa BinHeap
620    typedef BinHeap<Value, HeapCrossRef> Heap;
621   
622    /// \brief Instantiates a Heap.
623    ///
624    /// This function instantiates a \ref Heap.
625    /// \param crossref The cross reference of the heap.
626    static Heap *createHeap(HeapCrossRef& crossref) {
627      return new Heap(crossref);
628    }
629
630  };
631
632
633  /// \ingroup matching
634  ///
635  /// \brief Bipartite Max Weighted Matching algorithm
636  ///
637  /// This class implements the bipartite Max Weighted Matching
638  /// algorithm.  It uses the successive shortest path algorithm to
639  /// calculate the maximum weighted matching in the bipartite
640  /// graph. The algorithm can be used also to calculate the maximum
641  /// cardinality maximum weighted matching. The time complexity
642  /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
643  /// heap implementation but this can be improved to
644  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
645  ///
646  /// The algorithm also provides a potential function on the nodes
647  /// which a dual solution of the matching algorithm and it can be
648  /// used to proof the optimality of the given pimal solution.
649#ifdef DOXYGEN
650  template <typename _BpUGraph, typename _WeightMap, typename _Traits>
651#else
652  template <typename _BpUGraph,
653            typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
654            typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
655#endif
656  class MaxWeightedBipartiteMatching {
657  public:
658
659    typedef _Traits Traits;
660    typedef typename Traits::BpUGraph BpUGraph;
661    typedef typename Traits::WeightMap WeightMap;
662    typedef typename Traits::Value Value;
663
664  protected:
665
666    typedef typename Traits::HeapCrossRef HeapCrossRef;
667    typedef typename Traits::Heap Heap;
668
669   
670    typedef typename BpUGraph::Node Node;
671    typedef typename BpUGraph::ANodeIt ANodeIt;
672    typedef typename BpUGraph::BNodeIt BNodeIt;
673    typedef typename BpUGraph::UEdge UEdge;
674    typedef typename BpUGraph::UEdgeIt UEdgeIt;
675    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
676
677    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
678    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
679
680    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
681    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
682
683
684  public:
685
686    /// \brief \ref Exception for uninitialized parameters.
687    ///
688    /// This error represents problems in the initialization
689    /// of the parameters of the algorithms.
690    class UninitializedParameter : public lemon::UninitializedParameter {
691    public:
692      virtual const char* what() const throw() {
693        return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
694      }
695    };
696
697    ///\name Named template parameters
698
699    ///@{
700
701    template <class H, class CR>
702    struct DefHeapTraits : public Traits {
703      typedef CR HeapCrossRef;
704      typedef H Heap;
705      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
706        throw UninitializedParameter();
707      }
708      static Heap *createHeap(HeapCrossRef &) {
709        throw UninitializedParameter();
710      }
711    };
712
713    /// \brief \ref named-templ-param "Named parameter" for setting heap
714    /// and cross reference type
715    ///
716    /// \ref named-templ-param "Named parameter" for setting heap and cross
717    /// reference type
718    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
719    struct DefHeap
720      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
721                                            DefHeapTraits<H, CR> > {
722      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
723                                           DefHeapTraits<H, CR> > Create;
724    };
725
726    template <class H, class CR>
727    struct DefStandardHeapTraits : public Traits {
728      typedef CR HeapCrossRef;
729      typedef H Heap;
730      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
731        return new HeapCrossRef(graph);
732      }
733      static Heap *createHeap(HeapCrossRef &crossref) {
734        return new Heap(crossref);
735      }
736    };
737
738    /// \brief \ref named-templ-param "Named parameter" for setting heap and
739    /// cross reference type with automatic allocation
740    ///
741    /// \ref named-templ-param "Named parameter" for setting heap and cross
742    /// reference type. It can allocate the heap and the cross reference
743    /// object if the cross reference's constructor waits for the graph as
744    /// parameter and the heap's constructor waits for the cross reference.
745    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
746    struct DefStandardHeap
747      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
748                                            DefStandardHeapTraits<H, CR> > {
749      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
750                                           DefStandardHeapTraits<H, CR> >
751      Create;
752    };
753
754    ///@}
755
756
757    /// \brief Constructor.
758    ///
759    /// Constructor of the algorithm.
760    MaxWeightedBipartiteMatching(const BpUGraph& _graph,
761                                 const WeightMap& _weight)
762      : graph(&_graph), weight(&_weight),
763        anode_matching(_graph), bnode_matching(_graph),
764        anode_potential(_graph), bnode_potential(_graph),
765        _heap_cross_ref(0), local_heap_cross_ref(false),
766        _heap(0), local_heap(0) {}
767
768    /// \brief Destructor.
769    ///
770    /// Destructor of the algorithm.
771    ~MaxWeightedBipartiteMatching() {
772      destroyStructures();
773    }
774
775    /// \brief Sets the heap and the cross reference used by algorithm.
776    ///
777    /// Sets the heap and the cross reference used by algorithm.
778    /// If you don't use this function before calling \ref run(),
779    /// it will allocate one. The destuctor deallocates this
780    /// automatically allocated map, of course.
781    /// \return \c (*this)
782    MaxWeightedBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
783      if(local_heap_cross_ref) {
784        delete _heap_cross_ref;
785        local_heap_cross_ref = false;
786      }
787      _heap_cross_ref = &cr;
788      if(local_heap) {
789        delete _heap;
790        local_heap = false;
791      }
792      _heap = &hp;
793      return *this;
794    }
795
796    /// \name Execution control
797    /// The simplest way to execute the algorithm is to use
798    /// one of the member functions called \c run().
799    /// \n
800    /// If you need more control on the execution,
801    /// first you must call \ref init() or one alternative for it.
802    /// Finally \ref start() will perform the matching computation or
803    /// with step-by-step execution you can augment the solution.
804
805    /// @{
806
807    /// \brief Initalize the data structures.
808    ///
809    /// It initalizes the data structures and creates an empty matching.
810    void init() {
811      initStructures();
812      for (ANodeIt it(*graph); it != INVALID; ++it) {
813        anode_matching[it] = INVALID;
814        anode_potential[it] = 0;
815      }
816      for (BNodeIt it(*graph); it != INVALID; ++it) {
817        bnode_matching[it] = INVALID;
818        bnode_potential[it] = 0;
819        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
820          if ((*weight)[jt] > bnode_potential[it]) {
821            bnode_potential[it] = (*weight)[jt];
822          }
823        }
824      }
825      matching_value = 0;
826      matching_size = 0;
827    }
828
829
830    /// \brief An augmenting phase of the weighted matching algorithm
831    ///
832    /// It runs an augmenting phase of the weighted matching
833    /// algorithm. This phase finds the best augmenting path and
834    /// augments only on this paths.
835    ///
836    /// The algorithm consists at most
837    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
838    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
839    /// with binary heap.
840    /// \param decrease If the given parameter true the matching value
841    /// can be decreased in the augmenting phase. If we would like
842    /// to calculate the maximum cardinality maximum weighted matching
843    /// then we should let the algorithm to decrease the matching
844    /// value in order to increase the number of the matching edges.
845    bool augment(bool decrease = false) {
846
847      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
848      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
849
850      Node bestNode = INVALID;
851      Value bestValue = 0;
852
853      _heap->clear();
854      for (ANodeIt it(*graph); it != INVALID; ++it) {
855        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
856      }
857
858      for (ANodeIt it(*graph); it != INVALID; ++it) {
859        if (anode_matching[it] == INVALID) {
860          _heap->push(it, 0);
861        }
862      }
863
864      Value bdistMax = 0;
865      while (!_heap->empty()) {
866        Node anode = _heap->top();
867        Value avalue = _heap->prio();
868        _heap->pop();
869        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
870          if (jt == anode_matching[anode]) continue;
871          Node bnode = graph->bNode(jt);
872          Value bvalue = avalue  - (*weight)[jt] +
873            anode_potential[anode] + bnode_potential[bnode];
874          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
875            bdist[bnode] = bvalue;
876            bpred[bnode] = jt;
877          }
878          if (bvalue > bdistMax) {
879            bdistMax = bvalue;
880          }
881          if (bnode_matching[bnode] != INVALID) {
882            Node newanode = graph->aNode(bnode_matching[bnode]);
883            switch (_heap->state(newanode)) {
884            case Heap::PRE_HEAP:
885              _heap->push(newanode, bvalue);
886              break;
887            case Heap::IN_HEAP:
888              if (bvalue < (*_heap)[newanode]) {
889                _heap->decrease(newanode, bvalue);
890              }
891              break;
892            case Heap::POST_HEAP:
893              break;
894            }
895          } else {
896            if (bestNode == INVALID ||
897                bnode_potential[bnode] - bvalue > bestValue) {
898              bestValue = bnode_potential[bnode] - bvalue;
899              bestNode = bnode;
900            }
901          }
902        }
903      }
904
905      if (bestNode == INVALID || (!decrease && bestValue < 0)) {
906        return false;
907      }
908
909      matching_value += bestValue;
910      ++matching_size;
911
912      for (BNodeIt it(*graph); it != INVALID; ++it) {
913        if (bpred[it] != INVALID) {
914          bnode_potential[it] -= bdist[it];
915        } else {
916          bnode_potential[it] -= bdistMax;
917        }
918      }
919      for (ANodeIt it(*graph); it != INVALID; ++it) {
920        if (anode_matching[it] != INVALID) {
921          Node bnode = graph->bNode(anode_matching[it]);
922          if (bpred[bnode] != INVALID) {
923            anode_potential[it] += bdist[bnode];
924          } else {
925            anode_potential[it] += bdistMax;
926          }
927        }
928      }
929
930      while (bestNode != INVALID) {
931        UEdge uedge = bpred[bestNode];
932        Node anode = graph->aNode(uedge);
933       
934        bnode_matching[bestNode] = uedge;
935        if (anode_matching[anode] != INVALID) {
936          bestNode = graph->bNode(anode_matching[anode]);
937        } else {
938          bestNode = INVALID;
939        }
940        anode_matching[anode] = uedge;
941      }
942
943
944      return true;
945    }
946
947    /// \brief Starts the algorithm.
948    ///
949    /// Starts the algorithm. It runs augmenting phases until the
950    /// optimal solution reached.
951    ///
952    /// \param maxCardinality If the given value is true it will
953    /// calculate the maximum cardinality maximum matching instead of
954    /// the maximum matching.
955    void start(bool maxCardinality = false) {
956      while (augment(maxCardinality)) {}
957    }
958
959    /// \brief Runs the algorithm.
960    ///
961    /// It just initalize the algorithm and then start it.
962    ///
963    /// \param maxCardinality If the given value is true it will
964    /// calculate the maximum cardinality maximum matching instead of
965    /// the maximum matching.
966    void run(bool maxCardinality = false) {
967      init();
968      start(maxCardinality);
969    }
970
971    /// @}
972
973    /// \name Query Functions
974    /// The result of the %Matching algorithm can be obtained using these
975    /// functions.\n
976    /// Before the use of these functions,
977    /// either run() or start() must be called.
978   
979    ///@{
980
981    /// \brief Gives back the potential in the NodeMap
982    ///
983    /// Gives back the potential in the NodeMap. The matching is optimal
984    /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
985    /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
986    /// for each edges.
987    template <typename PotentialMap>
988    void potential(PotentialMap& pt) const {
989      for (ANodeIt it(*graph); it != INVALID; ++it) {
990        pt[it] = anode_potential[it];
991      }
992      for (BNodeIt it(*graph); it != INVALID; ++it) {
993        pt[it] = bnode_potential[it];
994      }
995    }
996
997    /// \brief Set true all matching uedge in the map.
998    ///
999    /// Set true all matching uedge in the map. It does not change the
1000    /// value mapped to the other uedges.
1001    /// \return The number of the matching edges.
1002    template <typename MatchingMap>
1003    int quickMatching(MatchingMap& mm) const {
1004      for (ANodeIt it(*graph); it != INVALID; ++it) {
1005        if (anode_matching[it] != INVALID) {
1006          mm[anode_matching[it]] = true;
1007        }
1008      }
1009      return matching_size;
1010    }
1011
1012    /// \brief Set true all matching uedge in the map and the others to false.
1013    ///
1014    /// Set true all matching uedge in the map and the others to false.
1015    /// \return The number of the matching edges.
1016    template <typename MatchingMap>
1017    int matching(MatchingMap& mm) const {
1018      for (UEdgeIt it(*graph); it != INVALID; ++it) {
1019        mm[it] = it == anode_matching[graph->aNode(it)];
1020      }
1021      return matching_size;
1022    }
1023
1024
1025    /// \brief Return true if the given uedge is in the matching.
1026    ///
1027    /// It returns true if the given uedge is in the matching.
1028    bool matchingEdge(const UEdge& edge) const {
1029      return anode_matching[graph->aNode(edge)] == edge;
1030    }
1031
1032    /// \brief Returns the matching edge from the node.
1033    ///
1034    /// Returns the matching edge from the node. If there is not such
1035    /// edge it gives back \c INVALID.
1036    UEdge matchingEdge(const Node& node) const {
1037      if (graph->aNode(node)) {
1038        return anode_matching[node];
1039      } else {
1040        return bnode_matching[node];
1041      }
1042    }
1043
1044    /// \brief Gives back the sum of weights of the matching edges.
1045    ///
1046    /// Gives back the sum of weights of the matching edges.
1047    Value matchingValue() const {
1048      return matching_value;
1049    }
1050
1051    /// \brief Gives back the number of the matching edges.
1052    ///
1053    /// Gives back the number of the matching edges.
1054    int matchingSize() const {
1055      return matching_size;
1056    }
1057
1058    /// @}
1059
1060  private:
1061
1062    void initStructures() {
1063      if (!_heap_cross_ref) {
1064        local_heap_cross_ref = true;
1065        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1066      }
1067      if (!_heap) {
1068        local_heap = true;
1069        _heap = Traits::createHeap(*_heap_cross_ref);
1070      }
1071    }
1072
1073    void destroyStructures() {
1074      if (local_heap_cross_ref) delete _heap_cross_ref;
1075      if (local_heap) delete _heap;
1076    }
1077
1078
1079  private:
1080   
1081    const BpUGraph *graph;
1082    const WeightMap* weight;
1083
1084    ANodeMatchingMap anode_matching;
1085    BNodeMatchingMap bnode_matching;
1086
1087    ANodePotentialMap anode_potential;
1088    BNodePotentialMap bnode_potential;
1089
1090    Value matching_value;
1091    int matching_size;
1092
1093    HeapCrossRef *_heap_cross_ref;
1094    bool local_heap_cross_ref;
1095
1096    Heap *_heap;
1097    bool local_heap;
1098 
1099  };
1100
1101  /// \ingroup matching
1102  ///
1103  /// \brief Maximum weighted bipartite matching
1104  ///
1105  /// This function calculates the maximum weighted matching
1106  /// in a bipartite graph. It gives back the matching in an undirected
1107  /// edge map.
1108  ///
1109  /// \param graph The bipartite graph.
1110  /// \param weight The undirected edge map which contains the weights.
1111  /// \retval matching The undirected edge map which will be set to
1112  /// the matching.
1113  /// \return The value of the matching.
1114  template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1115  typename WeightMap::Value
1116  maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,
1117                               MatchingMap& matching) {
1118    MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1119      bpmatching(graph, weight);
1120    bpmatching.run();
1121    bpmatching.matching(matching);
1122    return bpmatching.matchingValue();
1123  }
1124
1125  /// \ingroup matching
1126  ///
1127  /// \brief Maximum weighted maximum cardinality bipartite matching
1128  ///
1129  /// This function calculates the maximum weighted of the maximum cardinality
1130  /// matchings of a bipartite graph. It gives back the matching in an
1131  /// undirected edge map.
1132  ///
1133  /// \param graph The bipartite graph.
1134  /// \param weight The undirected edge map which contains the weights.
1135  /// \retval matching The undirected edge map which will be set to
1136  /// the matching.
1137  /// \return The value of the matching.
1138  template <typename BpUGraph, typename WeightMap, typename MatchingMap>
1139  typename WeightMap::Value
1140  maxWeightedMaxBipartiteMatching(const BpUGraph& graph,
1141                                  const WeightMap& weight,
1142                                  MatchingMap& matching) {
1143    MaxWeightedBipartiteMatching<BpUGraph, WeightMap>
1144      bpmatching(graph, weight);
1145    bpmatching.run(true);
1146    bpmatching.matching(matching);
1147    return bpmatching.matchingValue();
1148  }
1149
1150  /// \brief Default traits class for minimum cost bipartite matching
1151  /// algoritms.
1152  ///
1153  /// Default traits class for minimum cost bipartite matching
1154  /// algoritms. 
1155  ///
1156  /// \param _BpUGraph The bipartite undirected graph
1157  /// type. 
1158  ///
1159  /// \param _CostMap Type of cost map.
1160  template <typename _BpUGraph, typename _CostMap>
1161  struct MinCostMaxBipartiteMatchingDefaultTraits {
1162    /// \brief The type of the cost of the undirected edges.
1163    typedef typename _CostMap::Value Value;
1164
1165    /// The undirected bipartite graph type the algorithm runs on.
1166    typedef _BpUGraph BpUGraph;
1167
1168    /// The map of the edges costs
1169    typedef _CostMap CostMap;
1170
1171    /// \brief The cross reference type used by heap.
1172    ///
1173    /// The cross reference type used by heap.
1174    /// Usually it is \c Graph::NodeMap<int>.
1175    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
1176
1177    /// \brief Instantiates a HeapCrossRef.
1178    ///
1179    /// This function instantiates a \ref HeapCrossRef.
1180    /// \param graph is the graph, to which we would like to define the
1181    /// HeapCrossRef.
1182    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1183      return new HeapCrossRef(graph);
1184    }
1185   
1186    /// \brief The heap type used by costed matching algorithms.
1187    ///
1188    /// The heap type used by costed matching algorithms. It should
1189    /// minimize the priorities and the heap's key type is the graph's
1190    /// anode graph's node.
1191    ///
1192    /// \sa BinHeap
1193    typedef BinHeap<Value, HeapCrossRef> Heap;
1194   
1195    /// \brief Instantiates a Heap.
1196    ///
1197    /// This function instantiates a \ref Heap.
1198    /// \param crossref The cross reference of the heap.
1199    static Heap *createHeap(HeapCrossRef& crossref) {
1200      return new Heap(crossref);
1201    }
1202
1203  };
1204
1205
1206  /// \ingroup matching
1207  ///
1208  /// \brief Bipartite Min Cost Matching algorithm
1209  ///
1210  /// This class implements the bipartite Min Cost Matching algorithm.
1211  /// It uses the successive shortest path algorithm to calculate the
1212  /// minimum cost maximum matching in the bipartite graph. The time
1213  /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
1214  /// default binary heap implementation but this can be improved to
1215  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
1216  ///
1217  /// The algorithm also provides a potential function on the nodes
1218  /// which a dual solution of the matching algorithm and it can be
1219  /// used to proof the optimality of the given pimal solution.
1220#ifdef DOXYGEN
1221  template <typename _BpUGraph, typename _CostMap, typename _Traits>
1222#else
1223  template <typename _BpUGraph,
1224            typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
1225            typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
1226#endif
1227  class MinCostMaxBipartiteMatching {
1228  public:
1229
1230    typedef _Traits Traits;
1231    typedef typename Traits::BpUGraph BpUGraph;
1232    typedef typename Traits::CostMap CostMap;
1233    typedef typename Traits::Value Value;
1234
1235  protected:
1236
1237    typedef typename Traits::HeapCrossRef HeapCrossRef;
1238    typedef typename Traits::Heap Heap;
1239
1240   
1241    typedef typename BpUGraph::Node Node;
1242    typedef typename BpUGraph::ANodeIt ANodeIt;
1243    typedef typename BpUGraph::BNodeIt BNodeIt;
1244    typedef typename BpUGraph::UEdge UEdge;
1245    typedef typename BpUGraph::UEdgeIt UEdgeIt;
1246    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
1247
1248    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
1249    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
1250
1251    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
1252    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
1253
1254
1255  public:
1256
1257    /// \brief \ref Exception for uninitialized parameters.
1258    ///
1259    /// This error represents problems in the initialization
1260    /// of the parameters of the algorithms.
1261    class UninitializedParameter : public lemon::UninitializedParameter {
1262    public:
1263      virtual const char* what() const throw() {
1264        return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
1265      }
1266    };
1267
1268    ///\name Named template parameters
1269
1270    ///@{
1271
1272    template <class H, class CR>
1273    struct DefHeapTraits : public Traits {
1274      typedef CR HeapCrossRef;
1275      typedef H Heap;
1276      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
1277        throw UninitializedParameter();
1278      }
1279      static Heap *createHeap(HeapCrossRef &) {
1280        throw UninitializedParameter();
1281      }
1282    };
1283
1284    /// \brief \ref named-templ-param "Named parameter" for setting heap
1285    /// and cross reference type
1286    ///
1287    /// \ref named-templ-param "Named parameter" for setting heap and cross
1288    /// reference type
1289    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1290    struct DefHeap
1291      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1292                                            DefHeapTraits<H, CR> > {
1293      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1294                                           DefHeapTraits<H, CR> > Create;
1295    };
1296
1297    template <class H, class CR>
1298    struct DefStandardHeapTraits : public Traits {
1299      typedef CR HeapCrossRef;
1300      typedef H Heap;
1301      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1302        return new HeapCrossRef(graph);
1303      }
1304      static Heap *createHeap(HeapCrossRef &crossref) {
1305        return new Heap(crossref);
1306      }
1307    };
1308
1309    /// \brief \ref named-templ-param "Named parameter" for setting heap and
1310    /// cross reference type with automatic allocation
1311    ///
1312    /// \ref named-templ-param "Named parameter" for setting heap and cross
1313    /// reference type. It can allocate the heap and the cross reference
1314    /// object if the cross reference's constructor waits for the graph as
1315    /// parameter and the heap's constructor waits for the cross reference.
1316    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1317    struct DefStandardHeap
1318      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1319                                            DefStandardHeapTraits<H, CR> > {
1320      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1321                                           DefStandardHeapTraits<H, CR> >
1322      Create;
1323    };
1324
1325    ///@}
1326
1327
1328    /// \brief Constructor.
1329    ///
1330    /// Constructor of the algorithm.
1331    MinCostMaxBipartiteMatching(const BpUGraph& _graph,
1332                                 const CostMap& _cost)
1333      : graph(&_graph), cost(&_cost),
1334        anode_matching(_graph), bnode_matching(_graph),
1335        anode_potential(_graph), bnode_potential(_graph),
1336        _heap_cross_ref(0), local_heap_cross_ref(false),
1337        _heap(0), local_heap(0) {}
1338
1339    /// \brief Destructor.
1340    ///
1341    /// Destructor of the algorithm.
1342    ~MinCostMaxBipartiteMatching() {
1343      destroyStructures();
1344    }
1345
1346    /// \brief Sets the heap and the cross reference used by algorithm.
1347    ///
1348    /// Sets the heap and the cross reference used by algorithm.
1349    /// If you don't use this function before calling \ref run(),
1350    /// it will allocate one. The destuctor deallocates this
1351    /// automatically allocated map, of course.
1352    /// \return \c (*this)
1353    MinCostMaxBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
1354      if(local_heap_cross_ref) {
1355        delete _heap_cross_ref;
1356        local_heap_cross_ref = false;
1357      }
1358      _heap_cross_ref = &cr;
1359      if(local_heap) {
1360        delete _heap;
1361        local_heap = false;
1362      }
1363      _heap = &hp;
1364      return *this;
1365    }
1366
1367    /// \name Execution control
1368    /// The simplest way to execute the algorithm is to use
1369    /// one of the member functions called \c run().
1370    /// \n
1371    /// If you need more control on the execution,
1372    /// first you must call \ref init() or one alternative for it.
1373    /// Finally \ref start() will perform the matching computation or
1374    /// with step-by-step execution you can augment the solution.
1375
1376    /// @{
1377
1378    /// \brief Initalize the data structures.
1379    ///
1380    /// It initalizes the data structures and creates an empty matching.
1381    void init() {
1382      initStructures();
1383      for (ANodeIt it(*graph); it != INVALID; ++it) {
1384        anode_matching[it] = INVALID;
1385        anode_potential[it] = 0;
1386      }
1387      for (BNodeIt it(*graph); it != INVALID; ++it) {
1388        bnode_matching[it] = INVALID;
1389        bnode_potential[it] = 0;
1390      }
1391      matching_cost = 0;
1392      matching_size = 0;
1393    }
1394
1395
1396    /// \brief An augmenting phase of the costed matching algorithm
1397    ///
1398    /// It runs an augmenting phase of the matching algorithm. The
1399    /// phase finds the best augmenting path and augments only on this
1400    /// paths.
1401    ///
1402    /// The algorithm consists at most
1403    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
1404    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
1405    /// with binary heap.
1406    bool augment() {
1407
1408      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
1409      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1410
1411      Node bestNode = INVALID;
1412      Value bestValue = 0;
1413
1414      _heap->clear();
1415      for (ANodeIt it(*graph); it != INVALID; ++it) {
1416        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
1417      }
1418
1419      for (ANodeIt it(*graph); it != INVALID; ++it) {
1420        if (anode_matching[it] == INVALID) {
1421          _heap->push(it, 0);
1422        }
1423      }
1424      Value bdistMax = 0;
1425
1426      while (!_heap->empty()) {
1427        Node anode = _heap->top();
1428        Value avalue = _heap->prio();
1429        _heap->pop();
1430        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1431          if (jt == anode_matching[anode]) continue;
1432          Node bnode = graph->bNode(jt);
1433          Value bvalue = avalue + (*cost)[jt] +
1434            anode_potential[anode] - bnode_potential[bnode];
1435          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
1436            bdist[bnode] = bvalue;
1437            bpred[bnode] = jt;
1438          }
1439          if (bvalue > bdistMax) {
1440            bdistMax = bvalue;
1441          }
1442          if (bnode_matching[bnode] != INVALID) {
1443            Node newanode = graph->aNode(bnode_matching[bnode]);
1444            switch (_heap->state(newanode)) {
1445            case Heap::PRE_HEAP:
1446              _heap->push(newanode, bvalue);
1447              break;
1448            case Heap::IN_HEAP:
1449              if (bvalue < (*_heap)[newanode]) {
1450                _heap->decrease(newanode, bvalue);
1451              }
1452              break;
1453            case Heap::POST_HEAP:
1454              break;
1455            }
1456          } else {
1457            if (bestNode == INVALID ||
1458                bvalue + bnode_potential[bnode] < bestValue) {
1459              bestValue = bvalue + bnode_potential[bnode];
1460              bestNode = bnode;
1461            }
1462          }
1463        }
1464      }
1465
1466      if (bestNode == INVALID) {
1467        return false;
1468      }
1469
1470      matching_cost += bestValue;
1471      ++matching_size;
1472
1473      for (BNodeIt it(*graph); it != INVALID; ++it) {
1474        if (bpred[it] != INVALID) {
1475          bnode_potential[it] += bdist[it];
1476        } else {
1477          bnode_potential[it] += bdistMax;
1478        }
1479      }
1480      for (ANodeIt it(*graph); it != INVALID; ++it) {
1481        if (anode_matching[it] != INVALID) {
1482          Node bnode = graph->bNode(anode_matching[it]);
1483          if (bpred[bnode] != INVALID) {
1484            anode_potential[it] += bdist[bnode];
1485          } else {
1486            anode_potential[it] += bdistMax;
1487          }
1488        }
1489      }
1490
1491      while (bestNode != INVALID) {
1492        UEdge uedge = bpred[bestNode];
1493        Node anode = graph->aNode(uedge);
1494       
1495        bnode_matching[bestNode] = uedge;
1496        if (anode_matching[anode] != INVALID) {
1497          bestNode = graph->bNode(anode_matching[anode]);
1498        } else {
1499          bestNode = INVALID;
1500        }
1501        anode_matching[anode] = uedge;
1502      }
1503
1504
1505      return true;
1506    }
1507
1508    /// \brief Starts the algorithm.
1509    ///
1510    /// Starts the algorithm. It runs augmenting phases until the
1511    /// optimal solution reached.
1512    void start() {
1513      while (augment()) {}
1514    }
1515
1516    /// \brief Runs the algorithm.
1517    ///
1518    /// It just initalize the algorithm and then start it.
1519    void run() {
1520      init();
1521      start();
1522    }
1523
1524    /// @}
1525
1526    /// \name Query Functions
1527    /// The result of the %Matching algorithm can be obtained using these
1528    /// functions.\n
1529    /// Before the use of these functions,
1530    /// either run() or start() must be called.
1531   
1532    ///@{
1533
1534    /// \brief Gives back the potential in the NodeMap
1535    ///
1536    /// Gives back the potential in the NodeMap. The potential is optimal with
1537    /// the current number of edges if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
1538    /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
1539    /// for each edges.
1540    template <typename PotentialMap>
1541    void potential(PotentialMap& pt) const {
1542      for (ANodeIt it(*graph); it != INVALID; ++it) {
1543        pt[it] = anode_potential[it];
1544      }
1545      for (BNodeIt it(*graph); it != INVALID; ++it) {
1546        pt[it] = bnode_potential[it];
1547      }
1548    }
1549
1550    /// \brief Set true all matching uedge in the map.
1551    ///
1552    /// Set true all matching uedge in the map. It does not change the
1553    /// value mapped to the other uedges.
1554    /// \return The number of the matching edges.
1555    template <typename MatchingMap>
1556    int quickMatching(MatchingMap& mm) const {
1557      for (ANodeIt it(*graph); it != INVALID; ++it) {
1558        if (anode_matching[it] != INVALID) {
1559          mm[anode_matching[it]] = true;
1560        }
1561      }
1562      return matching_size;
1563    }
1564
1565    /// \brief Set true all matching uedge in the map and the others to false.
1566    ///
1567    /// Set true all matching uedge in the map and the others to false.
1568    /// \return The number of the matching edges.
1569    template <typename MatchingMap>
1570    int matching(MatchingMap& mm) const {
1571      for (UEdgeIt it(*graph); it != INVALID; ++it) {
1572        mm[it] = it == anode_matching[graph->aNode(it)];
1573      }
1574      return matching_size;
1575    }
1576
1577
1578    /// \brief Return true if the given uedge is in the matching.
1579    ///
1580    /// It returns true if the given uedge is in the matching.
1581    bool matchingEdge(const UEdge& edge) const {
1582      return anode_matching[graph->aNode(edge)] == edge;
1583    }
1584
1585    /// \brief Returns the matching edge from the node.
1586    ///
1587    /// Returns the matching edge from the node. If there is not such
1588    /// edge it gives back \c INVALID.
1589    UEdge matchingEdge(const Node& node) const {
1590      if (graph->aNode(node)) {
1591        return anode_matching[node];
1592      } else {
1593        return bnode_matching[node];
1594      }
1595    }
1596
1597    /// \brief Gives back the sum of costs of the matching edges.
1598    ///
1599    /// Gives back the sum of costs of the matching edges.
1600    Value matchingCost() const {
1601      return matching_cost;
1602    }
1603
1604    /// \brief Gives back the number of the matching edges.
1605    ///
1606    /// Gives back the number of the matching edges.
1607    int matchingSize() const {
1608      return matching_size;
1609    }
1610
1611    /// @}
1612
1613  private:
1614
1615    void initStructures() {
1616      if (!_heap_cross_ref) {
1617        local_heap_cross_ref = true;
1618        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1619      }
1620      if (!_heap) {
1621        local_heap = true;
1622        _heap = Traits::createHeap(*_heap_cross_ref);
1623      }
1624    }
1625
1626    void destroyStructures() {
1627      if (local_heap_cross_ref) delete _heap_cross_ref;
1628      if (local_heap) delete _heap;
1629    }
1630
1631
1632  private:
1633   
1634    const BpUGraph *graph;
1635    const CostMap* cost;
1636
1637    ANodeMatchingMap anode_matching;
1638    BNodeMatchingMap bnode_matching;
1639
1640    ANodePotentialMap anode_potential;
1641    BNodePotentialMap bnode_potential;
1642
1643    Value matching_cost;
1644    int matching_size;
1645
1646    HeapCrossRef *_heap_cross_ref;
1647    bool local_heap_cross_ref;
1648
1649    Heap *_heap;
1650    bool local_heap;
1651 
1652  };
1653
1654  /// \ingroup matching
1655  ///
1656  /// \brief Minimum cost maximum cardinality bipartite matching
1657  ///
1658  /// This function calculates the minimum cost matching of the maximum
1659  /// cardinality matchings of a bipartite graph. It gives back the matching
1660  /// in an undirected edge map.
1661  ///
1662  /// \param graph The bipartite graph.
1663  /// \param cost The undirected edge map which contains the costs.
1664  /// \retval matching The undirected edge map which will be set to
1665  /// the matching.
1666  /// \return The cost of the matching.
1667  template <typename BpUGraph, typename CostMap, typename MatchingMap>
1668  typename CostMap::Value
1669  minCostMaxBipartiteMatching(const BpUGraph& graph,
1670                              const CostMap& cost,
1671                              MatchingMap& matching) {
1672    MinCostMaxBipartiteMatching<BpUGraph, CostMap>
1673      bpmatching(graph, cost);
1674    bpmatching.run();
1675    bpmatching.matching(matching);
1676    return bpmatching.matchingCost();
1677  }
1678
1679}
1680
1681#endif
Note: See TracBrowser for help on using the repository browser.