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