deba@2040: /* -*- C++ -*-
deba@2040:  *
deba@2040:  * This file is a part of LEMON, a generic C++ optimization library
deba@2040:  *
alpar@2553:  * Copyright (C) 2003-2008
deba@2040:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@2040:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@2040:  *
deba@2040:  * Permission to use, modify and distribute this software is granted
deba@2040:  * provided that this copyright notice appears in all copies. For
deba@2040:  * precise terms see the accompanying LICENSE file.
deba@2040:  *
deba@2040:  * This software is provided "AS IS" with no warranty of any kind,
deba@2040:  * express or implied, and with no claim as to its suitability for any
deba@2040:  * purpose.
deba@2040:  *
deba@2040:  */
deba@2040: 
deba@2040: #ifndef LEMON_BIPARTITE_MATCHING
deba@2040: #define LEMON_BIPARTITE_MATCHING
deba@2040: 
deba@2051: #include <functional>
deba@2051: 
deba@2051: #include <lemon/bin_heap.h>
deba@2051: #include <lemon/maps.h>
deba@2040: 
deba@2040: #include <iostream>
deba@2040: 
deba@2040: ///\ingroup matching
deba@2040: ///\file
deba@2040: ///\brief Maximum matching algorithms in bipartite graphs.
deba@2462: ///
deba@2462: ///\note The pr_bipartite_matching.h file also contains algorithms to
deba@2462: ///solve maximum cardinality bipartite matching problems.
deba@2040: 
deba@2040: namespace lemon {
deba@2040: 
deba@2040:   /// \ingroup matching
deba@2040:   ///
deba@2040:   /// \brief Bipartite Max Cardinality Matching algorithm
deba@2040:   ///
deba@2040:   /// Bipartite Max Cardinality Matching algorithm. This class implements
deba@2051:   /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
deba@2040:   /// complexity.
deba@2462:   ///
deba@2462:   /// \note In several cases the push-relabel based algorithms have
deba@2462:   /// better runtime performance than the augmenting path based ones. 
deba@2462:   ///
deba@2462:   /// \see PrBipartiteMatching
deba@2040:   template <typename BpUGraph>
deba@2040:   class MaxBipartiteMatching {
deba@2040:   protected:
deba@2040: 
deba@2040:     typedef BpUGraph Graph;
deba@2040: 
deba@2040:     typedef typename Graph::Node Node;
deba@2040:     typedef typename Graph::ANodeIt ANodeIt;
deba@2040:     typedef typename Graph::BNodeIt BNodeIt;
deba@2040:     typedef typename Graph::UEdge UEdge;
deba@2040:     typedef typename Graph::UEdgeIt UEdgeIt;
deba@2040:     typedef typename Graph::IncEdgeIt IncEdgeIt;
deba@2040: 
deba@2040:     typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
deba@2040:     typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
deba@2040: 
deba@2040: 
deba@2040:   public:
deba@2040: 
deba@2040:     /// \brief Constructor.
deba@2040:     ///
deba@2040:     /// Constructor of the algorithm. 
deba@2466:     MaxBipartiteMatching(const BpUGraph& graph) 
deba@2466:       : _matching(graph), _rmatching(graph), _reached(graph), _graph(&graph) {}
deba@2040: 
deba@2040:     /// \name Execution control
deba@2040:     /// The simplest way to execute the algorithm is to use
deba@2040:     /// one of the member functions called \c run().
deba@2040:     /// \n
deba@2040:     /// If you need more control on the execution,
deba@2040:     /// first you must call \ref init() or one alternative for it.
deba@2040:     /// Finally \ref start() will perform the matching computation or
deba@2040:     /// with step-by-step execution you can augment the solution.
deba@2040: 
deba@2040:     /// @{
deba@2040: 
deba@2040:     /// \brief Initalize the data structures.
deba@2040:     ///
deba@2040:     /// It initalizes the data structures and creates an empty matching.
deba@2040:     void init() {
deba@2466:       for (ANodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         _matching.set(it, INVALID);
deba@2040:       }
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         _rmatching.set(it, INVALID);
deba@2466: 	_reached.set(it, -1);
deba@2040:       }
deba@2466:       _size = 0;
deba@2466:       _phase = -1;
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief Initalize the data structures.
deba@2040:     ///
deba@2040:     /// It initalizes the data structures and creates a greedy
deba@2040:     /// matching.  From this matching sometimes it is faster to get
deba@2040:     /// the matching than from the initial empty matching.
deba@2040:     void greedyInit() {
deba@2466:       _size = 0;
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         _rmatching.set(it, INVALID);
deba@2466: 	_reached.set(it, 0);
deba@2040:       }
deba@2466:       for (ANodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         _matching[it] = INVALID;
deba@2466:         for (IncEdgeIt jt(*_graph, it); jt != INVALID; ++jt) {
deba@2466:           if (_rmatching[_graph->bNode(jt)] == INVALID) {
deba@2466:             _matching.set(it, jt);
deba@2466: 	    _rmatching.set(_graph->bNode(jt), jt);
deba@2488: 	    _reached.set(_graph->bNode(jt), -1);
deba@2466:             ++_size;
deba@2040:             break;
deba@2040:           }
deba@2040:         }
deba@2040:       }
deba@2466:       _phase = 0;
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief Initalize the data structures with an initial matching.
deba@2040:     ///
deba@2040:     /// It initalizes the data structures with an initial matching.
deba@2040:     template <typename MatchingMap>
deba@2386:     void matchingInit(const MatchingMap& mm) {
deba@2466:       for (ANodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         _matching.set(it, INVALID);
deba@2040:       }
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         _rmatching.set(it, INVALID);
deba@2466: 	_reached.set(it, 0);
deba@2040:       }
deba@2466:       _size = 0;
deba@2466:       for (UEdgeIt it(*_graph); it != INVALID; ++it) {
deba@2386:         if (mm[it]) {
deba@2466:           ++_size;
deba@2466:           _matching.set(_graph->aNode(it), it);
deba@2466:           _rmatching.set(_graph->bNode(it), it);
deba@2466: 	  _reached.set(it, 0);
deba@2040:         }
deba@2040:       }
deba@2466:       _phase = 0;
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief Initalize the data structures with an initial matching.
deba@2040:     ///
deba@2040:     /// It initalizes the data structures with an initial matching.
deba@2040:     /// \return %True when the given map contains really a matching.
deba@2040:     template <typename MatchingMap>
deba@2466:     bool checkedMatchingInit(const MatchingMap& mm) {
deba@2466:       for (ANodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         _matching.set(it, INVALID);
deba@2040:       }
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         _rmatching.set(it, INVALID);
deba@2466: 	_reached.set(it, 0);
deba@2040:       }
deba@2466:       _size = 0;
deba@2466:       for (UEdgeIt it(*_graph); it != INVALID; ++it) {
deba@2386:         if (mm[it]) {
deba@2466:           ++_size;
deba@2466:           if (_matching[_graph->aNode(it)] != INVALID) {
deba@2040:             return false;
deba@2040:           }
deba@2466:           _matching.set(_graph->aNode(it), it);
deba@2466:           if (_matching[_graph->bNode(it)] != INVALID) {
deba@2040:             return false;
deba@2040:           }
deba@2466:           _matching.set(_graph->bNode(it), it);
deba@2466: 	  _reached.set(_graph->bNode(it), -1);
deba@2040:         }
deba@2040:       }
deba@2466:       _phase = 0;
deba@2466:       return true;
deba@2466:     }
deba@2466: 
deba@2466:   private:
deba@2466:     
deba@2466:     bool _find_path(Node anode, int maxlevel,
deba@2466: 		    typename Graph::template BNodeMap<int>& level) {
deba@2466:       for (IncEdgeIt it(*_graph, anode); it != INVALID; ++it) {
deba@2466: 	Node bnode = _graph->bNode(it); 
deba@2466: 	if (level[bnode] == maxlevel) {
deba@2466: 	  level.set(bnode, -1);
deba@2466: 	  if (maxlevel == 0) {
deba@2466: 	    _matching.set(anode, it);
deba@2466: 	    _rmatching.set(bnode, it);
deba@2466: 	    return true;
deba@2466: 	  } else {
deba@2466: 	    Node nnode = _graph->aNode(_rmatching[bnode]);
deba@2466: 	    if (_find_path(nnode, maxlevel - 1, level)) {
deba@2466: 	      _matching.set(anode, it);
deba@2466: 	      _rmatching.set(bnode, it);
deba@2466: 	      return true;
deba@2466: 	    }
deba@2466: 	  }
deba@2466: 	}
deba@2466:       }
deba@2040:       return false;
deba@2040:     }
deba@2040: 
deba@2466:   public:
deba@2466: 
deba@2040:     /// \brief An augmenting phase of the Hopcroft-Karp algorithm
deba@2040:     ///
deba@2040:     /// It runs an augmenting phase of the Hopcroft-Karp
deba@2466:     /// algorithm. This phase finds maximal edge disjoint augmenting
deba@2466:     /// paths and augments on these paths. The algorithm consists at
deba@2466:     /// most of \f$ O(\sqrt{n}) \f$ phase and one phase is \f$ O(e)
deba@2466:     /// \f$ long.
deba@2040:     bool augment() {
deba@2040: 
deba@2466:       ++_phase;
deba@2466:       
deba@2466:       typename Graph::template BNodeMap<int> _level(*_graph, -1);
deba@2466:       typename Graph::template ANodeMap<bool> _found(*_graph, false);
deba@2040: 
deba@2466:       std::vector<Node> queue, aqueue;
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         if (_rmatching[it] == INVALID) {
deba@2040:           queue.push_back(it);
deba@2466:           _reached.set(it, _phase);
deba@2466: 	  _level.set(it, 0);
deba@2040:         }
deba@2040:       }
deba@2040: 
deba@2040:       bool success = false;
deba@2040: 
deba@2466:       int level = 0;
deba@2040:       while (!success && !queue.empty()) {
deba@2466:         std::vector<Node> nqueue;
deba@2386:         for (int i = 0; i < int(queue.size()); ++i) {
deba@2466:           Node bnode = queue[i];
deba@2466:           for (IncEdgeIt jt(*_graph, bnode); jt != INVALID; ++jt) {
deba@2466:             Node anode = _graph->aNode(jt);
deba@2466:             if (_matching[anode] == INVALID) {
deba@2466: 
deba@2466: 	      if (!_found[anode]) {
deba@2466: 		if (_find_path(anode, level, _level)) {
deba@2466: 		  ++_size;
deba@2466: 		}
deba@2466: 		_found.set(anode, true);
deba@2466: 	      }
deba@2040:               success = true;
deba@2040:             } else {           
deba@2466:               Node nnode = _graph->bNode(_matching[anode]);
deba@2466:               if (_reached[nnode] != _phase) {
deba@2466:                 _reached.set(nnode, _phase);
deba@2466:                 nqueue.push_back(nnode);
deba@2466: 		_level.set(nnode, level + 1);
deba@2040:               }
deba@2040:             }
deba@2040:           }
deba@2040:         }
deba@2466: 	++level;
deba@2466:         queue.swap(nqueue);
deba@2040:       }
deba@2466:       
deba@2040:       return success;
deba@2040:     }
deba@2466:   private:
deba@2466:     
deba@2466:     void _find_path_bfs(Node anode,
deba@2466: 			typename Graph::template ANodeMap<UEdge>& pred) {
deba@2466:       while (true) {
deba@2466: 	UEdge uedge = pred[anode];
deba@2466: 	Node bnode = _graph->bNode(uedge);
deba@2040: 
deba@2466: 	UEdge nedge = _rmatching[bnode];
deba@2466: 
deba@2466: 	_matching.set(anode, uedge);
deba@2466: 	_rmatching.set(bnode, uedge);
deba@2466: 
deba@2466: 	if (nedge == INVALID) break;
deba@2466: 	anode = _graph->aNode(nedge);
deba@2466:       }
deba@2466:     }
deba@2466: 
deba@2466:   public:
deba@2466: 
deba@2466:     /// \brief An augmenting phase with single path augementing
deba@2040:     ///
deba@2466:     /// This phase finds only one augmenting paths and augments on
deba@2466:     /// these paths. The algorithm consists at most of \f$ O(n) \f$
deba@2466:     /// phase and one phase is \f$ O(e) \f$ long.
deba@2466:     bool simpleAugment() { 
deba@2466:       ++_phase;
deba@2466:       
deba@2466:       typename Graph::template ANodeMap<UEdge> _pred(*_graph);
deba@2040: 
deba@2466:       std::vector<Node> queue, aqueue;
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         if (_rmatching[it] == INVALID) {
deba@2040:           queue.push_back(it);
deba@2466:           _reached.set(it, _phase);
deba@2040:         }
deba@2040:       }
deba@2040: 
deba@2466:       bool success = false;
deba@2466: 
deba@2466:       int level = 0;
deba@2466:       while (!success && !queue.empty()) {
deba@2466:         std::vector<Node> nqueue;
deba@2386:         for (int i = 0; i < int(queue.size()); ++i) {
deba@2466:           Node bnode = queue[i];
deba@2466:           for (IncEdgeIt jt(*_graph, bnode); jt != INVALID; ++jt) {
deba@2466:             Node anode = _graph->aNode(jt);
deba@2466:             if (_matching[anode] == INVALID) {
deba@2466: 	      _pred.set(anode, jt);
deba@2466: 	      _find_path_bfs(anode, _pred);
deba@2466: 	      ++_size;
deba@2466: 	      return true;
deba@2040:             } else {           
deba@2466:               Node nnode = _graph->bNode(_matching[anode]);
deba@2466:               if (_reached[nnode] != _phase) {
deba@2466: 		_pred.set(anode, jt);
deba@2466: 		_reached.set(nnode, _phase);
deba@2466:                 nqueue.push_back(nnode);
deba@2040:               }
deba@2040:             }
deba@2040:           }
deba@2040:         }
deba@2466: 	++level;
deba@2466:         queue.swap(nqueue);
deba@2040:       }
deba@2040:       
deba@2466:       return success;
deba@2040:     }
deba@2040: 
deba@2466: 
deba@2466: 
deba@2040:     /// \brief Starts the algorithm.
deba@2040:     ///
deba@2040:     /// Starts the algorithm. It runs augmenting phases until the optimal
deba@2040:     /// solution reached.
deba@2040:     void start() {
deba@2040:       while (augment()) {}
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief Runs the algorithm.
deba@2040:     ///
deba@2040:     /// It just initalize the algorithm and then start it.
deba@2040:     void run() {
deba@2058:       greedyInit();
deba@2040:       start();
deba@2040:     }
deba@2040: 
deba@2040:     /// @}
deba@2040: 
deba@2040:     /// \name Query Functions
deba@2040:     /// The result of the %Matching algorithm can be obtained using these
deba@2040:     /// functions.\n
deba@2040:     /// Before the use of these functions,
deba@2040:     /// either run() or start() must be called.
deba@2040:     
deba@2040:     ///@{
deba@2040: 
deba@2466:     /// \brief Return true if the given uedge is in the matching.
deba@2466:     /// 
deba@2466:     /// It returns true if the given uedge is in the matching.
deba@2466:     bool matchingEdge(const UEdge& edge) const {
deba@2466:       return _matching[_graph->aNode(edge)] == edge;
deba@2466:     }
deba@2466: 
deba@2466:     /// \brief Returns the matching edge from the node.
deba@2466:     /// 
deba@2466:     /// Returns the matching edge from the node. If there is not such
deba@2466:     /// edge it gives back \c INVALID.
deba@2466:     /// \note If the parameter node is a B-node then the running time is
deba@2466:     /// propotional to the degree of the node.
deba@2466:     UEdge matchingEdge(const Node& node) const {
deba@2466:       if (_graph->aNode(node)) {
deba@2466:         return _matching[node];
deba@2466:       } else {
deba@2466:         return _rmatching[node];
deba@2466:       }
deba@2466:     }
deba@2466: 
deba@2462:     /// \brief Set true all matching uedge in the map.
deba@2462:     /// 
deba@2462:     /// Set true all matching uedge in the map. It does not change the
deba@2462:     /// value mapped to the other uedges.
deba@2462:     /// \return The number of the matching edges.
deba@2462:     template <typename MatchingMap>
deba@2462:     int quickMatching(MatchingMap& mm) const {
deba@2466:       for (ANodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         if (_matching[it] != INVALID) {
deba@2466:           mm.set(_matching[it], true);
deba@2462:         }
deba@2462:       }
deba@2466:       return _size;
deba@2462:     }
deba@2462: 
deba@2462:     /// \brief Set true all matching uedge in the map and the others to false.
deba@2462:     /// 
deba@2462:     /// Set true all matching uedge in the map and the others to false.
deba@2462:     /// \return The number of the matching edges.
deba@2462:     template <typename MatchingMap>
deba@2462:     int matching(MatchingMap& mm) const {
deba@2466:       for (UEdgeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         mm.set(it, it == _matching[_graph->aNode(it)]);
deba@2462:       }
deba@2466:       return _size;
deba@2462:     }
deba@2462: 
deba@2463:     ///Gives back the matching in an ANodeMap.
deba@2463: 
deba@2463:     ///Gives back the matching in an ANodeMap. The parameter should
deba@2463:     ///be a write ANodeMap of UEdge values.
deba@2463:     ///\return The number of the matching edges.
deba@2463:     template<class MatchingMap>
deba@2463:     int aMatching(MatchingMap& mm) const {
deba@2466:       for (ANodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         mm.set(it, _matching[it]);
deba@2463:       }
deba@2466:       return _size;
deba@2463:     }
deba@2463: 
deba@2463:     ///Gives back the matching in a BNodeMap.
deba@2463: 
deba@2463:     ///Gives back the matching in a BNodeMap. The parameter should
deba@2463:     ///be a write BNodeMap of UEdge values.
deba@2463:     ///\return The number of the matching edges.
deba@2463:     template<class MatchingMap>
deba@2463:     int bMatching(MatchingMap& mm) const {
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         mm.set(it, _rmatching[it]);
deba@2463:       }
deba@2466:       return _size;
deba@2462:     }
deba@2462: 
deba@2462:     /// \brief Returns a minimum covering of the nodes.
deba@2040:     ///
deba@2040:     /// The minimum covering set problem is the dual solution of the
deba@2462:     /// maximum bipartite matching. It provides a solution for this
deba@2040:     /// problem what is proof of the optimality of the matching.
deba@2040:     /// \return The size of the cover set.
deba@2040:     template <typename CoverMap>
deba@2058:     int coverSet(CoverMap& covering) const {
deba@2040: 
deba@2466:       int size = 0;
deba@2466:       for (ANodeIt it(*_graph); it != INVALID; ++it) {
deba@2466: 	bool cn = _matching[it] != INVALID && 
deba@2466: 	  _reached[_graph->bNode(_matching[it])] == _phase;
deba@2466:         covering.set(it, cn);
deba@2466:         if (cn) ++size;
deba@2040:       }
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466: 	bool cn = _reached[it] != _phase;
deba@2466:         covering.set(it, cn);
deba@2466:         if (cn) ++size;
deba@2040:       }
deba@2040:       return size;
deba@2040:     }
deba@2040: 
deba@2462:     /// \brief Gives back a barrier on the A-nodes
deba@2466:     ///    
deba@2462:     /// The barrier is s subset of the nodes on the same side of the
deba@2462:     /// graph, which size minus its neighbours is exactly the
deba@2462:     /// unmatched nodes on the A-side.  
deba@2462:     /// \retval barrier A WriteMap on the ANodes with bool value.
deba@2462:     template <typename BarrierMap>
deba@2462:     void aBarrier(BarrierMap& barrier) const {
deba@2462: 
deba@2466:       for (ANodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         barrier.set(it, _matching[it] == INVALID || 
deba@2466: 		    _reached[_graph->bNode(_matching[it])] != _phase);
deba@2040:       }
deba@2040:     }
deba@2040: 
deba@2462:     /// \brief Gives back a barrier on the B-nodes
deba@2466:     ///    
deba@2462:     /// The barrier is s subset of the nodes on the same side of the
deba@2462:     /// graph, which size minus its neighbours is exactly the
deba@2462:     /// unmatched nodes on the B-side.  
deba@2462:     /// \retval barrier A WriteMap on the BNodes with bool value.
deba@2462:     template <typename BarrierMap>
deba@2462:     void bBarrier(BarrierMap& barrier) const {
deba@2462: 
deba@2466:       for (BNodeIt it(*_graph); it != INVALID; ++it) {
deba@2466:         barrier.set(it, _reached[it] == _phase);
deba@2462:       }
deba@2466:     }
deba@2462: 
deba@2466:     /// \brief Gives back the number of the matching edges.
deba@2466:     ///
deba@2466:     /// Gives back the number of the matching edges.
deba@2466:     int matchingSize() const {
deba@2466:       return _size;
deba@2040:     }
deba@2040: 
deba@2040:     /// @}
deba@2040: 
deba@2040:   private:
deba@2040: 
deba@2466:     typename BpUGraph::template ANodeMap<UEdge> _matching;
deba@2466:     typename BpUGraph::template BNodeMap<UEdge> _rmatching;
deba@2040: 
deba@2466:     typename BpUGraph::template BNodeMap<int> _reached;
deba@2466: 
deba@2466:     int _phase;
deba@2466:     const Graph *_graph;
deba@2466: 
deba@2466:     int _size;
deba@2051:   
deba@2051:   };
deba@2051: 
deba@2058:   /// \ingroup matching
deba@2058:   ///
deba@2058:   /// \brief Maximum cardinality bipartite matching
deba@2058:   ///
deba@2058:   /// This function calculates the maximum cardinality matching
deba@2058:   /// in a bipartite graph. It gives back the matching in an undirected
deba@2058:   /// edge map.
deba@2058:   ///
deba@2058:   /// \param graph The bipartite graph.
deba@2463:   /// \return The size of the matching.
deba@2463:   template <typename BpUGraph>
deba@2463:   int maxBipartiteMatching(const BpUGraph& graph) {
deba@2463:     MaxBipartiteMatching<BpUGraph> bpmatching(graph);
deba@2463:     bpmatching.run();
deba@2463:     return bpmatching.matchingSize();
deba@2463:   }
deba@2463: 
deba@2463:   /// \ingroup matching
deba@2463:   ///
deba@2463:   /// \brief Maximum cardinality bipartite matching
deba@2463:   ///
deba@2463:   /// This function calculates the maximum cardinality matching
deba@2463:   /// in a bipartite graph. It gives back the matching in an undirected
deba@2463:   /// edge map.
deba@2463:   ///
deba@2463:   /// \param graph The bipartite graph.
deba@2463:   /// \retval matching The ANodeMap of UEdges which will be set to covered
deba@2463:   /// matching undirected edge.
deba@2058:   /// \return The size of the matching.
deba@2058:   template <typename BpUGraph, typename MatchingMap>
deba@2058:   int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) {
deba@2058:     MaxBipartiteMatching<BpUGraph> bpmatching(graph);
deba@2058:     bpmatching.run();
deba@2463:     bpmatching.aMatching(matching);
deba@2463:     return bpmatching.matchingSize();
deba@2463:   }
deba@2463: 
deba@2463:   /// \ingroup matching
deba@2463:   ///
deba@2463:   /// \brief Maximum cardinality bipartite matching
deba@2463:   ///
deba@2463:   /// This function calculates the maximum cardinality matching
deba@2463:   /// in a bipartite graph. It gives back the matching in an undirected
deba@2463:   /// edge map.
deba@2463:   ///
deba@2463:   /// \param graph The bipartite graph.
deba@2463:   /// \retval matching The ANodeMap of UEdges which will be set to covered
deba@2463:   /// matching undirected edge.
deba@2463:   /// \retval barrier The BNodeMap of bools which will be set to a barrier
deba@2463:   /// of the BNode-set.
deba@2463:   /// \return The size of the matching.
deba@2463:   template <typename BpUGraph, typename MatchingMap, typename BarrierMap>
deba@2463:   int maxBipartiteMatching(const BpUGraph& graph, 
deba@2463: 			   MatchingMap& matching, BarrierMap& barrier) {
deba@2463:     MaxBipartiteMatching<BpUGraph> bpmatching(graph);
deba@2463:     bpmatching.run();
deba@2463:     bpmatching.aMatching(matching);
deba@2463:     bpmatching.bBarrier(barrier);
deba@2058:     return bpmatching.matchingSize();
deba@2058:   }
deba@2058: 
deba@2051:   /// \brief Default traits class for weighted bipartite matching algoritms.
deba@2051:   ///
deba@2051:   /// Default traits class for weighted bipartite matching algoritms.
deba@2051:   /// \param _BpUGraph The bipartite undirected graph type.
deba@2051:   /// \param _WeightMap Type of weight map.
deba@2051:   template <typename _BpUGraph, typename _WeightMap>
deba@2550:   struct MaxWeightedBipartiteMatchingDefaultTraits {
deba@2051:     /// \brief The type of the weight of the undirected edges.
deba@2051:     typedef typename _WeightMap::Value Value;
deba@2051: 
deba@2051:     /// The undirected bipartite graph type the algorithm runs on. 
deba@2051:     typedef _BpUGraph BpUGraph;
deba@2051: 
deba@2051:     /// The map of the edges weights
deba@2051:     typedef _WeightMap WeightMap;
deba@2051: 
deba@2051:     /// \brief The cross reference type used by heap.
deba@2051:     ///
deba@2051:     /// The cross reference type used by heap.
deba@2550:     /// Usually it is \c Graph::ANodeMap<int>.
deba@2550:     typedef typename BpUGraph::template ANodeMap<int> HeapCrossRef;
deba@2051: 
deba@2051:     /// \brief Instantiates a HeapCrossRef.
deba@2051:     ///
deba@2051:     /// This function instantiates a \ref HeapCrossRef. 
deba@2051:     /// \param graph is the graph, to which we would like to define the 
deba@2051:     /// HeapCrossRef.
deba@2051:     static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
deba@2051:       return new HeapCrossRef(graph);
deba@2051:     }
deba@2051:     
deba@2051:     /// \brief The heap type used by weighted matching algorithms.
deba@2051:     ///
deba@2051:     /// The heap type used by weighted matching algorithms. It should
deba@2051:     /// minimize the priorities and the heap's key type is the graph's
deba@2051:     /// anode graph's node.
deba@2051:     ///
deba@2051:     /// \sa BinHeap
mqrelly@2263:     typedef BinHeap<Value, HeapCrossRef> Heap;
deba@2051:     
deba@2051:     /// \brief Instantiates a Heap.
deba@2051:     ///
deba@2051:     /// This function instantiates a \ref Heap. 
deba@2051:     /// \param crossref The cross reference of the heap.
deba@2051:     static Heap *createHeap(HeapCrossRef& crossref) {
deba@2051:       return new Heap(crossref);
deba@2051:     }
deba@2051: 
deba@2051:   };
deba@2051: 
deba@2051: 
deba@2051:   /// \ingroup matching
deba@2051:   ///
deba@2051:   /// \brief Bipartite Max Weighted Matching algorithm
deba@2051:   ///
deba@2051:   /// This class implements the bipartite Max Weighted Matching
deba@2051:   /// algorithm.  It uses the successive shortest path algorithm to
deba@2051:   /// calculate the maximum weighted matching in the bipartite
deba@2051:   /// graph. The algorithm can be used also to calculate the maximum
deba@2051:   /// cardinality maximum weighted matching. The time complexity
deba@2051:   /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
deba@2051:   /// heap implementation but this can be improved to 
deba@2051:   /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
deba@2051:   ///
deba@2051:   /// The algorithm also provides a potential function on the nodes
deba@2051:   /// which a dual solution of the matching algorithm and it can be
deba@2051:   /// used to proof the optimality of the given pimal solution.
deba@2051: #ifdef DOXYGEN
deba@2051:   template <typename _BpUGraph, typename _WeightMap, typename _Traits>
deba@2051: #else
deba@2051:   template <typename _BpUGraph, 
deba@2051:             typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
deba@2550:             typename _Traits = 
deba@2550: 	    MaxWeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
deba@2051: #endif
deba@2051:   class MaxWeightedBipartiteMatching {
deba@2051:   public:
deba@2051: 
deba@2051:     typedef _Traits Traits;
deba@2051:     typedef typename Traits::BpUGraph BpUGraph;
deba@2051:     typedef typename Traits::WeightMap WeightMap;
deba@2051:     typedef typename Traits::Value Value;
deba@2051: 
deba@2051:   protected:
deba@2051: 
deba@2051:     typedef typename Traits::HeapCrossRef HeapCrossRef;
deba@2051:     typedef typename Traits::Heap Heap; 
deba@2051: 
deba@2051:     
deba@2051:     typedef typename BpUGraph::Node Node;
deba@2051:     typedef typename BpUGraph::ANodeIt ANodeIt;
deba@2051:     typedef typename BpUGraph::BNodeIt BNodeIt;
deba@2051:     typedef typename BpUGraph::UEdge UEdge;
deba@2051:     typedef typename BpUGraph::UEdgeIt UEdgeIt;
deba@2051:     typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
deba@2051: 
deba@2051:     typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
deba@2051:     typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
deba@2051: 
deba@2051:     typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
deba@2051:     typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
deba@2051: 
deba@2051: 
deba@2051:   public:
deba@2051: 
deba@2051:     /// \brief \ref Exception for uninitialized parameters.
deba@2051:     ///
deba@2051:     /// This error represents problems in the initialization
deba@2051:     /// of the parameters of the algorithms.
deba@2051:     class UninitializedParameter : public lemon::UninitializedParameter {
deba@2051:     public:
alpar@2151:       virtual const char* what() const throw() {
deba@2051: 	return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
deba@2051:       }
deba@2051:     };
deba@2051: 
deba@2051:     ///\name Named template parameters
deba@2051: 
deba@2051:     ///@{
deba@2051: 
deba@2051:     template <class H, class CR>
deba@2051:     struct DefHeapTraits : public Traits {
deba@2051:       typedef CR HeapCrossRef;
deba@2051:       typedef H Heap;
deba@2051:       static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
deba@2051: 	throw UninitializedParameter();
deba@2051:       }
deba@2051:       static Heap *createHeap(HeapCrossRef &) {
deba@2051: 	throw UninitializedParameter();
deba@2051:       }
deba@2051:     };
deba@2051: 
deba@2051:     /// \brief \ref named-templ-param "Named parameter" for setting heap 
deba@2051:     /// and cross reference type
deba@2051:     ///
deba@2051:     /// \ref named-templ-param "Named parameter" for setting heap and cross 
deba@2051:     /// reference type
deba@2051:     template <class H, class CR = typename BpUGraph::template NodeMap<int> >
deba@2051:     struct DefHeap
deba@2051:       : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
deba@2051:                                             DefHeapTraits<H, CR> > { 
deba@2051:       typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
deba@2051:                                            DefHeapTraits<H, CR> > Create;
deba@2051:     };
deba@2051: 
deba@2051:     template <class H, class CR>
deba@2051:     struct DefStandardHeapTraits : public Traits {
deba@2051:       typedef CR HeapCrossRef;
deba@2051:       typedef H Heap;
deba@2051:       static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
deba@2051: 	return new HeapCrossRef(graph);
deba@2051:       }
deba@2051:       static Heap *createHeap(HeapCrossRef &crossref) {
deba@2051: 	return new Heap(crossref);
deba@2051:       }
deba@2051:     };
deba@2051: 
deba@2051:     /// \brief \ref named-templ-param "Named parameter" for setting heap and 
deba@2051:     /// cross reference type with automatic allocation
deba@2051:     ///
deba@2051:     /// \ref named-templ-param "Named parameter" for setting heap and cross 
deba@2051:     /// reference type. It can allocate the heap and the cross reference 
deba@2051:     /// object if the cross reference's constructor waits for the graph as 
deba@2051:     /// parameter and the heap's constructor waits for the cross reference.
deba@2051:     template <class H, class CR = typename BpUGraph::template NodeMap<int> >
deba@2051:     struct DefStandardHeap
deba@2051:       : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
deba@2051:                                             DefStandardHeapTraits<H, CR> > { 
deba@2051:       typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
deba@2051:                                            DefStandardHeapTraits<H, CR> > 
deba@2051:       Create;
deba@2051:     };
deba@2051: 
deba@2051:     ///@}
deba@2051: 
deba@2051: 
deba@2051:     /// \brief Constructor.
deba@2051:     ///
deba@2051:     /// Constructor of the algorithm. 
deba@2051:     MaxWeightedBipartiteMatching(const BpUGraph& _graph, 
deba@2051:                                  const WeightMap& _weight) 
deba@2051:       : graph(&_graph), weight(&_weight),
deba@2051:         anode_matching(_graph), bnode_matching(_graph),
deba@2051:         anode_potential(_graph), bnode_potential(_graph),
deba@2051:         _heap_cross_ref(0), local_heap_cross_ref(false),
deba@2051:         _heap(0), local_heap(0) {}
deba@2051: 
deba@2051:     /// \brief Destructor.
deba@2051:     ///
deba@2051:     /// Destructor of the algorithm.
deba@2051:     ~MaxWeightedBipartiteMatching() {
deba@2051:       destroyStructures();
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Sets the heap and the cross reference used by algorithm.
deba@2051:     ///
deba@2051:     /// Sets the heap and the cross reference used by algorithm.
deba@2051:     /// If you don't use this function before calling \ref run(),
deba@2051:     /// it will allocate one. The destuctor deallocates this
deba@2051:     /// automatically allocated map, of course.
deba@2051:     /// \return \c (*this)
deba@2386:     MaxWeightedBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
deba@2051:       if(local_heap_cross_ref) {
deba@2051: 	delete _heap_cross_ref;
deba@2051: 	local_heap_cross_ref = false;
deba@2051:       }
deba@2386:       _heap_cross_ref = &cr;
deba@2051:       if(local_heap) {
deba@2051: 	delete _heap;
deba@2051: 	local_heap = false;
deba@2051:       }
deba@2386:       _heap = &hp;
deba@2051:       return *this;
deba@2051:     }
deba@2051: 
deba@2051:     /// \name Execution control
deba@2051:     /// The simplest way to execute the algorithm is to use
deba@2051:     /// one of the member functions called \c run().
deba@2051:     /// \n
deba@2051:     /// If you need more control on the execution,
deba@2051:     /// first you must call \ref init() or one alternative for it.
deba@2051:     /// Finally \ref start() will perform the matching computation or
deba@2051:     /// with step-by-step execution you can augment the solution.
deba@2051: 
deba@2051:     /// @{
deba@2051: 
deba@2051:     /// \brief Initalize the data structures.
deba@2051:     ///
deba@2051:     /// It initalizes the data structures and creates an empty matching.
deba@2051:     void init() {
deba@2051:       initStructures();
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         anode_matching[it] = INVALID;
deba@2051:         anode_potential[it] = 0;
deba@2051:       }
deba@2051:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         bnode_matching[it] = INVALID;
deba@2051:         bnode_potential[it] = 0;
deba@2051:         for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
deba@2058:           if ((*weight)[jt] > bnode_potential[it]) {
deba@2058:             bnode_potential[it] = (*weight)[jt];
deba@2051:           }
deba@2051:         }
deba@2051:       }
deba@2051:       matching_value = 0;
deba@2051:       matching_size = 0;
deba@2051:     }
deba@2051: 
deba@2051: 
deba@2051:     /// \brief An augmenting phase of the weighted matching algorithm
deba@2051:     ///
deba@2051:     /// It runs an augmenting phase of the weighted matching 
alpar@2352:     /// algorithm. This phase finds the best augmenting path and 
deba@2051:     /// augments only on this paths. 
deba@2051:     ///
deba@2051:     /// The algorithm consists at most 
deba@2051:     /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 
deba@2051:     /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 
deba@2051:     /// with binary heap.
deba@2051:     /// \param decrease If the given parameter true the matching value
deba@2051:     /// can be decreased in the augmenting phase. If we would like
deba@2051:     /// to calculate the maximum cardinality maximum weighted matching
deba@2051:     /// then we should let the algorithm to decrease the matching
deba@2051:     /// value in order to increase the number of the matching edges.
deba@2051:     bool augment(bool decrease = false) {
deba@2051: 
deba@2051:       typename BpUGraph::template BNodeMap<Value> bdist(*graph);
deba@2051:       typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
deba@2051: 
deba@2051:       Node bestNode = INVALID;
deba@2051:       Value bestValue = 0;
deba@2051: 
deba@2051:       _heap->clear();
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
deba@2051:       }
deba@2051: 
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (anode_matching[it] == INVALID) {
deba@2051:           _heap->push(it, 0);
deba@2051:         }
deba@2051:       }
deba@2051: 
deba@2051:       Value bdistMax = 0;
deba@2051:       while (!_heap->empty()) {
deba@2051:         Node anode = _heap->top();
deba@2051:         Value avalue = _heap->prio();
deba@2051:         _heap->pop();
deba@2051:         for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
deba@2051:           if (jt == anode_matching[anode]) continue;
deba@2051:           Node bnode = graph->bNode(jt);
deba@2058:           Value bvalue = avalue  - (*weight)[jt] +
deba@2058:             anode_potential[anode] + bnode_potential[bnode];
deba@2550:           if (bvalue > bdistMax) {
deba@2550:             bdistMax = bvalue;
deba@2550:           }
deba@2051:           if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
deba@2051:             bdist[bnode] = bvalue;
deba@2051:             bpred[bnode] = jt;
deba@2550:           } else continue;
deba@2051:           if (bnode_matching[bnode] != INVALID) {
deba@2051:             Node newanode = graph->aNode(bnode_matching[bnode]);
deba@2051:             switch (_heap->state(newanode)) {
deba@2051:             case Heap::PRE_HEAP:
deba@2051:               _heap->push(newanode, bvalue);
deba@2051:               break;
deba@2051:             case Heap::IN_HEAP:
deba@2051:               if (bvalue < (*_heap)[newanode]) {
deba@2051:                 _heap->decrease(newanode, bvalue);
deba@2051:               }
deba@2051:               break;
deba@2051:             case Heap::POST_HEAP:
deba@2051:               break;
deba@2051:             }
deba@2051:           } else {
deba@2051:             if (bestNode == INVALID || 
deba@2058:                 bnode_potential[bnode] - bvalue > bestValue) {
deba@2058:               bestValue = bnode_potential[bnode] - bvalue;
deba@2051:               bestNode = bnode;
deba@2051:             }
deba@2051:           }
deba@2051:         }
deba@2051:       }
deba@2051: 
deba@2051:       if (bestNode == INVALID || (!decrease && bestValue < 0)) {
deba@2051:         return false;
deba@2051:       }
deba@2051: 
deba@2051:       matching_value += bestValue;
deba@2051:       ++matching_size;
deba@2051: 
deba@2051:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (bpred[it] != INVALID) {
deba@2058:           bnode_potential[it] -= bdist[it];
deba@2051:         } else {
deba@2058:           bnode_potential[it] -= bdistMax;
deba@2051:         }
deba@2051:       }
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (anode_matching[it] != INVALID) {
deba@2051:           Node bnode = graph->bNode(anode_matching[it]);
deba@2051:           if (bpred[bnode] != INVALID) {
deba@2051:             anode_potential[it] += bdist[bnode];
deba@2051:           } else {
deba@2051:             anode_potential[it] += bdistMax;
deba@2051:           }
deba@2051:         }
deba@2051:       }
deba@2051: 
deba@2051:       while (bestNode != INVALID) {
deba@2051:         UEdge uedge = bpred[bestNode];
deba@2051:         Node anode = graph->aNode(uedge);
deba@2051:         
deba@2051:         bnode_matching[bestNode] = uedge;
deba@2051:         if (anode_matching[anode] != INVALID) {
deba@2051:           bestNode = graph->bNode(anode_matching[anode]);
deba@2051:         } else {
deba@2051:           bestNode = INVALID;
deba@2051:         }
deba@2051:         anode_matching[anode] = uedge;
deba@2051:       }
deba@2051: 
deba@2051: 
deba@2051:       return true;
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Starts the algorithm.
deba@2051:     ///
deba@2051:     /// Starts the algorithm. It runs augmenting phases until the
deba@2051:     /// optimal solution reached.
deba@2051:     ///
deba@2051:     /// \param maxCardinality If the given value is true it will
deba@2051:     /// calculate the maximum cardinality maximum matching instead of
deba@2051:     /// the maximum matching.
deba@2051:     void start(bool maxCardinality = false) {
deba@2051:       while (augment(maxCardinality)) {}
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Runs the algorithm.
deba@2051:     ///
deba@2051:     /// It just initalize the algorithm and then start it.
deba@2051:     ///
deba@2051:     /// \param maxCardinality If the given value is true it will
deba@2051:     /// calculate the maximum cardinality maximum matching instead of
deba@2051:     /// the maximum matching.
deba@2051:     void run(bool maxCardinality = false) {
deba@2051:       init();
deba@2051:       start(maxCardinality);
deba@2051:     }
deba@2051: 
deba@2051:     /// @}
deba@2051: 
deba@2051:     /// \name Query Functions
deba@2051:     /// The result of the %Matching algorithm can be obtained using these
deba@2051:     /// functions.\n
deba@2051:     /// Before the use of these functions,
deba@2051:     /// either run() or start() must be called.
deba@2051:     
deba@2051:     ///@{
deba@2051: 
deba@2051:     /// \brief Gives back the potential in the NodeMap
deba@2051:     ///
deba@2058:     /// Gives back the potential in the NodeMap. The matching is optimal
deba@2058:     /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
deba@2058:     /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
deba@2058:     /// for each edges. 
deba@2051:     template <typename PotentialMap>
deba@2386:     void potential(PotentialMap& pt) const {
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2463:         pt.set(it, anode_potential[it]);
deba@2051:       }
deba@2051:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2463:         pt.set(it, bnode_potential[it]);
deba@2051:       }
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Set true all matching uedge in the map.
deba@2051:     /// 
deba@2051:     /// Set true all matching uedge in the map. It does not change the
deba@2051:     /// value mapped to the other uedges.
deba@2051:     /// \return The number of the matching edges.
deba@2051:     template <typename MatchingMap>
deba@2386:     int quickMatching(MatchingMap& mm) const {
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (anode_matching[it] != INVALID) {
deba@2463:           mm.set(anode_matching[it], true);
deba@2051:         }
deba@2051:       }
deba@2051:       return matching_size;
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Set true all matching uedge in the map and the others to false.
deba@2051:     /// 
deba@2051:     /// Set true all matching uedge in the map and the others to false.
deba@2051:     /// \return The number of the matching edges.
deba@2051:     template <typename MatchingMap>
deba@2386:     int matching(MatchingMap& mm) const {
deba@2051:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2463:         mm.set(it, it == anode_matching[graph->aNode(it)]);
deba@2463:       }
deba@2463:       return matching_size;
deba@2463:     }
deba@2463: 
deba@2463:     ///Gives back the matching in an ANodeMap.
deba@2463: 
deba@2463:     ///Gives back the matching in an ANodeMap. The parameter should
deba@2463:     ///be a write ANodeMap of UEdge values.
deba@2463:     ///\return The number of the matching edges.
deba@2463:     template<class MatchingMap>
deba@2463:     int aMatching(MatchingMap& mm) const {
deba@2463:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2463:         mm.set(it, anode_matching[it]);
deba@2463:       }
deba@2463:       return matching_size;
deba@2463:     }
deba@2463: 
deba@2463:     ///Gives back the matching in a BNodeMap.
deba@2463: 
deba@2463:     ///Gives back the matching in a BNodeMap. The parameter should
deba@2463:     ///be a write BNodeMap of UEdge values.
deba@2463:     ///\return The number of the matching edges.
deba@2463:     template<class MatchingMap>
deba@2463:     int bMatching(MatchingMap& mm) const {
deba@2463:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2463:         mm.set(it, bnode_matching[it]);
deba@2051:       }
deba@2051:       return matching_size;
deba@2051:     }
deba@2051: 
deba@2051: 
deba@2051:     /// \brief Return true if the given uedge is in the matching.
deba@2051:     /// 
deba@2051:     /// It returns true if the given uedge is in the matching.
deba@2058:     bool matchingEdge(const UEdge& edge) const {
deba@2051:       return anode_matching[graph->aNode(edge)] == edge;
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Returns the matching edge from the node.
deba@2051:     /// 
deba@2051:     /// Returns the matching edge from the node. If there is not such
deba@2051:     /// edge it gives back \c INVALID.
deba@2058:     UEdge matchingEdge(const Node& node) const {
deba@2051:       if (graph->aNode(node)) {
deba@2051:         return anode_matching[node];
deba@2051:       } else {
deba@2051:         return bnode_matching[node];
deba@2051:       }
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Gives back the sum of weights of the matching edges.
deba@2051:     ///
deba@2051:     /// Gives back the sum of weights of the matching edges.
deba@2051:     Value matchingValue() const {
deba@2051:       return matching_value;
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Gives back the number of the matching edges.
deba@2051:     ///
deba@2051:     /// Gives back the number of the matching edges.
deba@2051:     int matchingSize() const {
deba@2051:       return matching_size;
deba@2051:     }
deba@2051: 
deba@2051:     /// @}
deba@2051: 
deba@2051:   private:
deba@2051: 
deba@2051:     void initStructures() {
deba@2051:       if (!_heap_cross_ref) {
deba@2051: 	local_heap_cross_ref = true;
deba@2051: 	_heap_cross_ref = Traits::createHeapCrossRef(*graph);
deba@2051:       }
deba@2051:       if (!_heap) {
deba@2051: 	local_heap = true;
deba@2051: 	_heap = Traits::createHeap(*_heap_cross_ref);
deba@2051:       }
deba@2051:     }
deba@2051: 
deba@2051:     void destroyStructures() {
deba@2051:       if (local_heap_cross_ref) delete _heap_cross_ref;
deba@2051:       if (local_heap) delete _heap;
deba@2051:     }
deba@2051: 
deba@2051: 
deba@2051:   private:
deba@2051:     
deba@2051:     const BpUGraph *graph;
deba@2051:     const WeightMap* weight;
deba@2051: 
deba@2051:     ANodeMatchingMap anode_matching;
deba@2051:     BNodeMatchingMap bnode_matching;
deba@2051: 
deba@2051:     ANodePotentialMap anode_potential;
deba@2051:     BNodePotentialMap bnode_potential;
deba@2051: 
deba@2051:     Value matching_value;
deba@2051:     int matching_size;
deba@2051: 
deba@2051:     HeapCrossRef *_heap_cross_ref;
deba@2051:     bool local_heap_cross_ref;
deba@2051: 
deba@2051:     Heap *_heap;
deba@2051:     bool local_heap;
deba@2051:   
deba@2051:   };
deba@2051: 
deba@2058:   /// \ingroup matching
deba@2058:   ///
deba@2058:   /// \brief Maximum weighted bipartite matching
deba@2058:   ///
deba@2058:   /// This function calculates the maximum weighted matching
deba@2058:   /// in a bipartite graph. It gives back the matching in an undirected
deba@2058:   /// edge map.
deba@2058:   ///
deba@2058:   /// \param graph The bipartite graph.
deba@2058:   /// \param weight The undirected edge map which contains the weights.
deba@2058:   /// \retval matching The undirected edge map which will be set to 
deba@2058:   /// the matching.
deba@2058:   /// \return The value of the matching.
deba@2058:   template <typename BpUGraph, typename WeightMap, typename MatchingMap>
deba@2058:   typename WeightMap::Value 
deba@2058:   maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,
deba@2058:                                MatchingMap& matching) {
deba@2058:     MaxWeightedBipartiteMatching<BpUGraph, WeightMap> 
deba@2058:       bpmatching(graph, weight);
deba@2058:     bpmatching.run();
deba@2058:     bpmatching.matching(matching);
deba@2058:     return bpmatching.matchingValue();
deba@2058:   }
deba@2058: 
deba@2058:   /// \ingroup matching
deba@2058:   ///
deba@2058:   /// \brief Maximum weighted maximum cardinality bipartite matching
deba@2058:   ///
deba@2058:   /// This function calculates the maximum weighted of the maximum cardinality
deba@2058:   /// matchings of a bipartite graph. It gives back the matching in an 
deba@2058:   /// undirected edge map.
deba@2058:   ///
deba@2058:   /// \param graph The bipartite graph.
deba@2058:   /// \param weight The undirected edge map which contains the weights.
deba@2058:   /// \retval matching The undirected edge map which will be set to 
deba@2058:   /// the matching.
deba@2058:   /// \return The value of the matching.
deba@2058:   template <typename BpUGraph, typename WeightMap, typename MatchingMap>
deba@2058:   typename WeightMap::Value 
deba@2058:   maxWeightedMaxBipartiteMatching(const BpUGraph& graph, 
deba@2058:                                   const WeightMap& weight,
deba@2058:                                   MatchingMap& matching) {
deba@2058:     MaxWeightedBipartiteMatching<BpUGraph, WeightMap> 
deba@2058:       bpmatching(graph, weight);
deba@2058:     bpmatching.run(true);
deba@2058:     bpmatching.matching(matching);
deba@2058:     return bpmatching.matchingValue();
deba@2058:   }
deba@2058: 
deba@2051:   /// \brief Default traits class for minimum cost bipartite matching
deba@2051:   /// algoritms.
deba@2051:   ///
deba@2051:   /// Default traits class for minimum cost bipartite matching
deba@2051:   /// algoritms.  
deba@2051:   ///
deba@2051:   /// \param _BpUGraph The bipartite undirected graph
deba@2051:   /// type.  
deba@2051:   ///
deba@2051:   /// \param _CostMap Type of cost map.
deba@2051:   template <typename _BpUGraph, typename _CostMap>
deba@2051:   struct MinCostMaxBipartiteMatchingDefaultTraits {
deba@2051:     /// \brief The type of the cost of the undirected edges.
deba@2051:     typedef typename _CostMap::Value Value;
deba@2051: 
deba@2051:     /// The undirected bipartite graph type the algorithm runs on. 
deba@2051:     typedef _BpUGraph BpUGraph;
deba@2051: 
deba@2051:     /// The map of the edges costs
deba@2051:     typedef _CostMap CostMap;
deba@2051: 
deba@2051:     /// \brief The cross reference type used by heap.
deba@2051:     ///
deba@2051:     /// The cross reference type used by heap.
deba@2051:     /// Usually it is \c Graph::NodeMap<int>.
deba@2051:     typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
deba@2051: 
deba@2051:     /// \brief Instantiates a HeapCrossRef.
deba@2051:     ///
deba@2051:     /// This function instantiates a \ref HeapCrossRef. 
deba@2051:     /// \param graph is the graph, to which we would like to define the 
deba@2051:     /// HeapCrossRef.
deba@2051:     static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
deba@2051:       return new HeapCrossRef(graph);
deba@2051:     }
deba@2051:     
deba@2051:     /// \brief The heap type used by costed matching algorithms.
deba@2051:     ///
deba@2051:     /// The heap type used by costed matching algorithms. It should
deba@2051:     /// minimize the priorities and the heap's key type is the graph's
deba@2051:     /// anode graph's node.
deba@2051:     ///
deba@2051:     /// \sa BinHeap
deba@2269:     typedef BinHeap<Value, HeapCrossRef> Heap;
deba@2051:     
deba@2051:     /// \brief Instantiates a Heap.
deba@2051:     ///
deba@2051:     /// This function instantiates a \ref Heap. 
deba@2051:     /// \param crossref The cross reference of the heap.
deba@2051:     static Heap *createHeap(HeapCrossRef& crossref) {
deba@2051:       return new Heap(crossref);
deba@2051:     }
deba@2051: 
deba@2051:   };
deba@2051: 
deba@2051: 
deba@2051:   /// \ingroup matching
deba@2051:   ///
deba@2051:   /// \brief Bipartite Min Cost Matching algorithm
deba@2051:   ///
deba@2051:   /// This class implements the bipartite Min Cost Matching algorithm.
deba@2051:   /// It uses the successive shortest path algorithm to calculate the
deba@2051:   /// minimum cost maximum matching in the bipartite graph. The time
deba@2051:   /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
deba@2051:   /// default binary heap implementation but this can be improved to
deba@2051:   /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
deba@2051:   ///
deba@2051:   /// The algorithm also provides a potential function on the nodes
deba@2051:   /// which a dual solution of the matching algorithm and it can be
deba@2051:   /// used to proof the optimality of the given pimal solution.
deba@2051: #ifdef DOXYGEN
deba@2051:   template <typename _BpUGraph, typename _CostMap, typename _Traits>
deba@2051: #else
deba@2051:   template <typename _BpUGraph, 
deba@2051:             typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
deba@2550:             typename _Traits = 
deba@2550: 	    MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
deba@2051: #endif
deba@2051:   class MinCostMaxBipartiteMatching {
deba@2051:   public:
deba@2051: 
deba@2051:     typedef _Traits Traits;
deba@2051:     typedef typename Traits::BpUGraph BpUGraph;
deba@2051:     typedef typename Traits::CostMap CostMap;
deba@2051:     typedef typename Traits::Value Value;
deba@2051: 
deba@2051:   protected:
deba@2051: 
deba@2051:     typedef typename Traits::HeapCrossRef HeapCrossRef;
deba@2051:     typedef typename Traits::Heap Heap; 
deba@2051: 
deba@2051:     
deba@2051:     typedef typename BpUGraph::Node Node;
deba@2051:     typedef typename BpUGraph::ANodeIt ANodeIt;
deba@2051:     typedef typename BpUGraph::BNodeIt BNodeIt;
deba@2051:     typedef typename BpUGraph::UEdge UEdge;
deba@2051:     typedef typename BpUGraph::UEdgeIt UEdgeIt;
deba@2051:     typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
deba@2051: 
deba@2051:     typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
deba@2051:     typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
deba@2051: 
deba@2051:     typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
deba@2051:     typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
deba@2051: 
deba@2051: 
deba@2051:   public:
deba@2051: 
deba@2051:     /// \brief \ref Exception for uninitialized parameters.
deba@2051:     ///
deba@2051:     /// This error represents problems in the initialization
deba@2051:     /// of the parameters of the algorithms.
deba@2051:     class UninitializedParameter : public lemon::UninitializedParameter {
deba@2051:     public:
alpar@2151:       virtual const char* what() const throw() {
deba@2051: 	return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
deba@2051:       }
deba@2051:     };
deba@2051: 
deba@2051:     ///\name Named template parameters
deba@2051: 
deba@2051:     ///@{
deba@2051: 
deba@2051:     template <class H, class CR>
deba@2051:     struct DefHeapTraits : public Traits {
deba@2051:       typedef CR HeapCrossRef;
deba@2051:       typedef H Heap;
deba@2051:       static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
deba@2051: 	throw UninitializedParameter();
deba@2051:       }
deba@2051:       static Heap *createHeap(HeapCrossRef &) {
deba@2051: 	throw UninitializedParameter();
deba@2051:       }
deba@2051:     };
deba@2051: 
deba@2051:     /// \brief \ref named-templ-param "Named parameter" for setting heap 
deba@2051:     /// and cross reference type
deba@2051:     ///
deba@2051:     /// \ref named-templ-param "Named parameter" for setting heap and cross 
deba@2051:     /// reference type
deba@2051:     template <class H, class CR = typename BpUGraph::template NodeMap<int> >
deba@2051:     struct DefHeap
deba@2051:       : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
deba@2051:                                             DefHeapTraits<H, CR> > { 
deba@2051:       typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
deba@2051:                                            DefHeapTraits<H, CR> > Create;
deba@2051:     };
deba@2051: 
deba@2051:     template <class H, class CR>
deba@2051:     struct DefStandardHeapTraits : public Traits {
deba@2051:       typedef CR HeapCrossRef;
deba@2051:       typedef H Heap;
deba@2051:       static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
deba@2051: 	return new HeapCrossRef(graph);
deba@2051:       }
deba@2051:       static Heap *createHeap(HeapCrossRef &crossref) {
deba@2051: 	return new Heap(crossref);
deba@2051:       }
deba@2051:     };
deba@2051: 
deba@2051:     /// \brief \ref named-templ-param "Named parameter" for setting heap and 
deba@2051:     /// cross reference type with automatic allocation
deba@2051:     ///
deba@2051:     /// \ref named-templ-param "Named parameter" for setting heap and cross 
deba@2051:     /// reference type. It can allocate the heap and the cross reference 
deba@2051:     /// object if the cross reference's constructor waits for the graph as 
deba@2051:     /// parameter and the heap's constructor waits for the cross reference.
deba@2051:     template <class H, class CR = typename BpUGraph::template NodeMap<int> >
deba@2051:     struct DefStandardHeap
deba@2051:       : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
deba@2051:                                             DefStandardHeapTraits<H, CR> > { 
deba@2051:       typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
deba@2051:                                            DefStandardHeapTraits<H, CR> > 
deba@2051:       Create;
deba@2051:     };
deba@2051: 
deba@2051:     ///@}
deba@2051: 
deba@2051: 
deba@2051:     /// \brief Constructor.
deba@2051:     ///
deba@2051:     /// Constructor of the algorithm. 
deba@2051:     MinCostMaxBipartiteMatching(const BpUGraph& _graph, 
deba@2051:                                  const CostMap& _cost) 
deba@2051:       : graph(&_graph), cost(&_cost),
deba@2051:         anode_matching(_graph), bnode_matching(_graph),
deba@2051:         anode_potential(_graph), bnode_potential(_graph),
deba@2051:         _heap_cross_ref(0), local_heap_cross_ref(false),
deba@2051:         _heap(0), local_heap(0) {}
deba@2051: 
deba@2051:     /// \brief Destructor.
deba@2051:     ///
deba@2051:     /// Destructor of the algorithm.
deba@2051:     ~MinCostMaxBipartiteMatching() {
deba@2051:       destroyStructures();
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Sets the heap and the cross reference used by algorithm.
deba@2051:     ///
deba@2051:     /// Sets the heap and the cross reference used by algorithm.
deba@2051:     /// If you don't use this function before calling \ref run(),
deba@2051:     /// it will allocate one. The destuctor deallocates this
deba@2051:     /// automatically allocated map, of course.
deba@2051:     /// \return \c (*this)
deba@2386:     MinCostMaxBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
deba@2051:       if(local_heap_cross_ref) {
deba@2051: 	delete _heap_cross_ref;
deba@2051: 	local_heap_cross_ref = false;
deba@2051:       }
deba@2386:       _heap_cross_ref = &cr;
deba@2051:       if(local_heap) {
deba@2051: 	delete _heap;
deba@2051: 	local_heap = false;
deba@2051:       }
deba@2386:       _heap = &hp;
deba@2051:       return *this;
deba@2051:     }
deba@2051: 
deba@2051:     /// \name Execution control
deba@2051:     /// The simplest way to execute the algorithm is to use
deba@2051:     /// one of the member functions called \c run().
deba@2051:     /// \n
deba@2051:     /// If you need more control on the execution,
deba@2051:     /// first you must call \ref init() or one alternative for it.
deba@2051:     /// Finally \ref start() will perform the matching computation or
deba@2051:     /// with step-by-step execution you can augment the solution.
deba@2051: 
deba@2051:     /// @{
deba@2051: 
deba@2051:     /// \brief Initalize the data structures.
deba@2051:     ///
deba@2051:     /// It initalizes the data structures and creates an empty matching.
deba@2051:     void init() {
deba@2051:       initStructures();
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         anode_matching[it] = INVALID;
deba@2051:         anode_potential[it] = 0;
deba@2051:       }
deba@2051:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         bnode_matching[it] = INVALID;
deba@2051:         bnode_potential[it] = 0;
deba@2051:       }
deba@2051:       matching_cost = 0;
deba@2051:       matching_size = 0;
deba@2051:     }
deba@2051: 
deba@2051: 
deba@2051:     /// \brief An augmenting phase of the costed matching algorithm
deba@2051:     ///
deba@2051:     /// It runs an augmenting phase of the matching algorithm. The
deba@2051:     /// phase finds the best augmenting path and augments only on this
deba@2051:     /// paths.
deba@2051:     ///
deba@2051:     /// The algorithm consists at most 
deba@2051:     /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 
deba@2051:     /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 
deba@2051:     /// with binary heap.
deba@2051:     bool augment() {
deba@2051: 
deba@2051:       typename BpUGraph::template BNodeMap<Value> bdist(*graph);
deba@2051:       typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
deba@2051: 
deba@2051:       Node bestNode = INVALID;
deba@2051:       Value bestValue = 0;
deba@2051: 
deba@2051:       _heap->clear();
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
deba@2051:       }
deba@2051: 
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (anode_matching[it] == INVALID) {
deba@2051:           _heap->push(it, 0);
deba@2051:         }
deba@2051:       }
deba@2136:       Value bdistMax = 0;
deba@2051: 
deba@2051:       while (!_heap->empty()) {
deba@2051:         Node anode = _heap->top();
deba@2051:         Value avalue = _heap->prio();
deba@2051:         _heap->pop();
deba@2051:         for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
deba@2051:           if (jt == anode_matching[anode]) continue;
deba@2051:           Node bnode = graph->bNode(jt);
deba@2051:           Value bvalue = avalue + (*cost)[jt] + 
deba@2051:             anode_potential[anode] - bnode_potential[bnode];
deba@2550:           if (bvalue > bdistMax) {
deba@2550:             bdistMax = bvalue;
deba@2550:           }
deba@2051:           if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
deba@2051:             bdist[bnode] = bvalue;
deba@2051:             bpred[bnode] = jt;
deba@2550:           } else continue;
deba@2051:           if (bnode_matching[bnode] != INVALID) {
deba@2051:             Node newanode = graph->aNode(bnode_matching[bnode]);
deba@2051:             switch (_heap->state(newanode)) {
deba@2051:             case Heap::PRE_HEAP:
deba@2051:               _heap->push(newanode, bvalue);
deba@2051:               break;
deba@2051:             case Heap::IN_HEAP:
deba@2051:               if (bvalue < (*_heap)[newanode]) {
deba@2051:                 _heap->decrease(newanode, bvalue);
deba@2051:               }
deba@2051:               break;
deba@2051:             case Heap::POST_HEAP:
deba@2051:               break;
deba@2051:             }
deba@2051:           } else {
deba@2051:             if (bestNode == INVALID || 
deba@2051:                 bvalue + bnode_potential[bnode] < bestValue) {
deba@2051:               bestValue = bvalue + bnode_potential[bnode];
deba@2051:               bestNode = bnode;
deba@2051:             }
deba@2051:           }
deba@2051:         }
deba@2051:       }
deba@2051: 
deba@2051:       if (bestNode == INVALID) {
deba@2051:         return false;
deba@2051:       }
deba@2051: 
deba@2051:       matching_cost += bestValue;
deba@2051:       ++matching_size;
deba@2051: 
deba@2051:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (bpred[it] != INVALID) {
deba@2051:           bnode_potential[it] += bdist[it];
deba@2136:         } else {
deba@2136:           bnode_potential[it] += bdistMax;
deba@2051:         }
deba@2051:       }
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (anode_matching[it] != INVALID) {
deba@2051:           Node bnode = graph->bNode(anode_matching[it]);
deba@2051:           if (bpred[bnode] != INVALID) {
deba@2051:             anode_potential[it] += bdist[bnode];
deba@2136:           } else {
deba@2136:             anode_potential[it] += bdistMax;
deba@2051:           }
deba@2051:         }
deba@2051:       }
deba@2051: 
deba@2051:       while (bestNode != INVALID) {
deba@2051:         UEdge uedge = bpred[bestNode];
deba@2051:         Node anode = graph->aNode(uedge);
deba@2051:         
deba@2051:         bnode_matching[bestNode] = uedge;
deba@2051:         if (anode_matching[anode] != INVALID) {
deba@2051:           bestNode = graph->bNode(anode_matching[anode]);
deba@2051:         } else {
deba@2051:           bestNode = INVALID;
deba@2051:         }
deba@2051:         anode_matching[anode] = uedge;
deba@2051:       }
deba@2051: 
deba@2051: 
deba@2051:       return true;
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Starts the algorithm.
deba@2051:     ///
deba@2051:     /// Starts the algorithm. It runs augmenting phases until the
deba@2051:     /// optimal solution reached.
deba@2051:     void start() {
deba@2051:       while (augment()) {}
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Runs the algorithm.
deba@2051:     ///
deba@2051:     /// It just initalize the algorithm and then start it.
deba@2051:     void run() {
deba@2051:       init();
deba@2051:       start();
deba@2051:     }
deba@2051: 
deba@2051:     /// @}
deba@2051: 
deba@2051:     /// \name Query Functions
deba@2051:     /// The result of the %Matching algorithm can be obtained using these
deba@2051:     /// functions.\n
deba@2051:     /// Before the use of these functions,
deba@2051:     /// either run() or start() must be called.
deba@2051:     
deba@2051:     ///@{
deba@2051: 
deba@2051:     /// \brief Gives back the potential in the NodeMap
deba@2051:     ///
deba@2463:     /// Gives back the potential in the NodeMap. The matching is optimal
deba@2463:     /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
deba@2463:     /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
deba@2463:     /// for each edges. 
deba@2051:     template <typename PotentialMap>
deba@2386:     void potential(PotentialMap& pt) const {
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2463:         pt.set(it, anode_potential[it]);
deba@2051:       }
deba@2051:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2463:         pt.set(it, bnode_potential[it]);
deba@2051:       }
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Set true all matching uedge in the map.
deba@2051:     /// 
deba@2051:     /// Set true all matching uedge in the map. It does not change the
deba@2051:     /// value mapped to the other uedges.
deba@2051:     /// \return The number of the matching edges.
deba@2051:     template <typename MatchingMap>
deba@2386:     int quickMatching(MatchingMap& mm) const {
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (anode_matching[it] != INVALID) {
deba@2463:           mm.set(anode_matching[it], true);
deba@2051:         }
deba@2051:       }
deba@2051:       return matching_size;
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Set true all matching uedge in the map and the others to false.
deba@2051:     /// 
deba@2051:     /// Set true all matching uedge in the map and the others to false.
deba@2051:     /// \return The number of the matching edges.
deba@2051:     template <typename MatchingMap>
deba@2386:     int matching(MatchingMap& mm) const {
deba@2051:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2463:         mm.set(it, it == anode_matching[graph->aNode(it)]);
deba@2051:       }
deba@2051:       return matching_size;
deba@2051:     }
deba@2051: 
deba@2463:     /// \brief Gives back the matching in an ANodeMap.
deba@2463:     ///
deba@2463:     /// Gives back the matching in an ANodeMap. The parameter should
deba@2463:     /// be a write ANodeMap of UEdge values.
deba@2463:     /// \return The number of the matching edges.
deba@2463:     template<class MatchingMap>
deba@2463:     int aMatching(MatchingMap& mm) const {
deba@2463:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2463:         mm.set(it, anode_matching[it]);
deba@2463:       }
deba@2463:       return matching_size;
deba@2463:     }
deba@2463: 
deba@2463:     /// \brief Gives back the matching in a BNodeMap.
deba@2463:     ///
deba@2463:     /// Gives back the matching in a BNodeMap. The parameter should
deba@2463:     /// be a write BNodeMap of UEdge values.
deba@2463:     /// \return The number of the matching edges.
deba@2463:     template<class MatchingMap>
deba@2463:     int bMatching(MatchingMap& mm) const {
deba@2463:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2463:         mm.set(it, bnode_matching[it]);
deba@2463:       }
deba@2463:       return matching_size;
deba@2463:     }
deba@2051: 
deba@2051:     /// \brief Return true if the given uedge is in the matching.
deba@2051:     /// 
deba@2051:     /// It returns true if the given uedge is in the matching.
deba@2058:     bool matchingEdge(const UEdge& edge) const {
deba@2051:       return anode_matching[graph->aNode(edge)] == edge;
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Returns the matching edge from the node.
deba@2051:     /// 
deba@2051:     /// Returns the matching edge from the node. If there is not such
deba@2051:     /// edge it gives back \c INVALID.
deba@2058:     UEdge matchingEdge(const Node& node) const {
deba@2051:       if (graph->aNode(node)) {
deba@2051:         return anode_matching[node];
deba@2051:       } else {
deba@2051:         return bnode_matching[node];
deba@2051:       }
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Gives back the sum of costs of the matching edges.
deba@2051:     ///
deba@2051:     /// Gives back the sum of costs of the matching edges.
deba@2051:     Value matchingCost() const {
deba@2051:       return matching_cost;
deba@2051:     }
deba@2051: 
deba@2051:     /// \brief Gives back the number of the matching edges.
deba@2051:     ///
deba@2051:     /// Gives back the number of the matching edges.
deba@2051:     int matchingSize() const {
deba@2051:       return matching_size;
deba@2051:     }
deba@2051: 
deba@2051:     /// @}
deba@2051: 
deba@2051:   private:
deba@2051: 
deba@2051:     void initStructures() {
deba@2051:       if (!_heap_cross_ref) {
deba@2051: 	local_heap_cross_ref = true;
deba@2051: 	_heap_cross_ref = Traits::createHeapCrossRef(*graph);
deba@2051:       }
deba@2051:       if (!_heap) {
deba@2051: 	local_heap = true;
deba@2051: 	_heap = Traits::createHeap(*_heap_cross_ref);
deba@2051:       }
deba@2051:     }
deba@2051: 
deba@2051:     void destroyStructures() {
deba@2051:       if (local_heap_cross_ref) delete _heap_cross_ref;
deba@2051:       if (local_heap) delete _heap;
deba@2051:     }
deba@2051: 
deba@2051: 
deba@2051:   private:
deba@2051:     
deba@2051:     const BpUGraph *graph;
deba@2051:     const CostMap* cost;
deba@2051: 
deba@2051:     ANodeMatchingMap anode_matching;
deba@2051:     BNodeMatchingMap bnode_matching;
deba@2051: 
deba@2051:     ANodePotentialMap anode_potential;
deba@2051:     BNodePotentialMap bnode_potential;
deba@2051: 
deba@2051:     Value matching_cost;
deba@2051:     int matching_size;
deba@2051: 
deba@2051:     HeapCrossRef *_heap_cross_ref;
deba@2051:     bool local_heap_cross_ref;
deba@2051: 
deba@2051:     Heap *_heap;
deba@2051:     bool local_heap;
deba@2040:   
deba@2040:   };
deba@2040: 
deba@2058:   /// \ingroup matching
deba@2058:   ///
deba@2058:   /// \brief Minimum cost maximum cardinality bipartite matching
deba@2058:   ///
deba@2463:   /// This function calculates the maximum cardinality matching with
deba@2463:   /// minimum cost of a bipartite graph. It gives back the matching in
deba@2463:   /// an undirected edge map.
deba@2058:   ///
deba@2058:   /// \param graph The bipartite graph.
deba@2058:   /// \param cost The undirected edge map which contains the costs.
deba@2058:   /// \retval matching The undirected edge map which will be set to 
deba@2058:   /// the matching.
deba@2058:   /// \return The cost of the matching.
deba@2058:   template <typename BpUGraph, typename CostMap, typename MatchingMap>
deba@2058:   typename CostMap::Value 
deba@2058:   minCostMaxBipartiteMatching(const BpUGraph& graph, 
deba@2058:                               const CostMap& cost,
deba@2058:                               MatchingMap& matching) {
deba@2058:     MinCostMaxBipartiteMatching<BpUGraph, CostMap> 
deba@2058:       bpmatching(graph, cost);
deba@2058:     bpmatching.run();
deba@2058:     bpmatching.matching(matching);
deba@2058:     return bpmatching.matchingCost();
deba@2058:   }
deba@2058: 
deba@2040: }
deba@2040: 
deba@2040: #endif