deba@2040: /* -*- C++ -*-
deba@2040:  *
deba@2040:  * This file is a part of LEMON, a generic C++ optimization library
deba@2040:  *
deba@2040:  * Copyright (C) 2003-2006
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@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@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@2040:     MaxBipartiteMatching(const BpUGraph& _graph) 
deba@2040:       : anode_matching(_graph), bnode_matching(_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@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         anode_matching[it] = INVALID;
deba@2040:       }
deba@2040:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         bnode_matching[it] = INVALID;
deba@2040:       }
deba@2051:       matching_size = 0;
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@2051:       matching_size = 0;
deba@2040:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         bnode_matching[it] = INVALID;
deba@2040:       }
deba@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         anode_matching[it] = INVALID;
deba@2040:         for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
deba@2040:           if (bnode_matching[graph->bNode(jt)] == INVALID) {
deba@2040:             anode_matching[it] = jt;
deba@2040:             bnode_matching[graph->bNode(jt)] = jt;
deba@2051:             ++matching_size;
deba@2040:             break;
deba@2040:           }
deba@2040:         }
deba@2040:       }
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@2040:     void matchingInit(const MatchingMap& matching) {
deba@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         anode_matching[it] = INVALID;
deba@2040:       }
deba@2040:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         bnode_matching[it] = INVALID;
deba@2040:       }
deba@2051:       matching_size = 0;
deba@2040:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2040:         if (matching[it]) {
deba@2051:           ++matching_size;
deba@2040:           anode_matching[graph->aNode(it)] = it;
deba@2040:           bnode_matching[graph->bNode(it)] = it;
deba@2040:         }
deba@2040:       }
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@2040:     void checkedMatchingInit(const MatchingMap& matching) {
deba@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         anode_matching[it] = INVALID;
deba@2040:       }
deba@2040:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         bnode_matching[it] = INVALID;
deba@2040:       }
deba@2051:       matching_size = 0;
deba@2040:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2040:         if (matching[it]) {
deba@2051:           ++matching_size;
deba@2040:           if (anode_matching[graph->aNode(it)] != INVALID) {
deba@2040:             return false;
deba@2040:           }
deba@2040:           anode_matching[graph->aNode(it)] = it;
deba@2040:           if (bnode_matching[graph->aNode(it)] != INVALID) {
deba@2040:             return false;
deba@2040:           }
deba@2040:           bnode_matching[graph->bNode(it)] = it;
deba@2040:         }
deba@2040:       }
deba@2040:       return false;
deba@2040:     }
deba@2040: 
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@2040:     /// algorithm. The phase finds maximum count of edge disjoint
deba@2040:     /// augmenting paths and augments on these paths. The algorithm
deba@2040:     /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
deba@2040:     /// \f$ O(e) \f$ long.
deba@2040:     bool augment() {
deba@2040: 
deba@2040:       typename Graph::template ANodeMap<bool> areached(*graph, false);
deba@2040:       typename Graph::template BNodeMap<bool> breached(*graph, false);
deba@2040: 
deba@2040:       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
deba@2040: 
deba@2040:       std::vector<Node> queue, bqueue;
deba@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         if (anode_matching[it] == INVALID) {
deba@2040:           queue.push_back(it);
deba@2040:           areached[it] = true;
deba@2040:         }
deba@2040:       }
deba@2040: 
deba@2040:       bool success = false;
deba@2040: 
deba@2040:       while (!success && !queue.empty()) {
deba@2040:         std::vector<Node> newqueue;
deba@2040:         for (int i = 0; i < (int)queue.size(); ++i) {
deba@2040:           Node anode = queue[i];
deba@2040:           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
deba@2040:             Node bnode = graph->bNode(jt);
deba@2040:             if (breached[bnode]) continue;
deba@2040:             breached[bnode] = true;
deba@2040:             bpred[bnode] = jt;
deba@2040:             if (bnode_matching[bnode] == INVALID) {
deba@2040:               bqueue.push_back(bnode);
deba@2040:               success = true;
deba@2040:             } else {           
deba@2040:               Node newanode = graph->aNode(bnode_matching[bnode]);
deba@2040:               if (!areached[newanode]) {
deba@2040:                 areached[newanode] = true;
deba@2040:                 newqueue.push_back(newanode);
deba@2040:               }
deba@2040:             }
deba@2040:           }
deba@2040:         }
deba@2040:         queue.swap(newqueue);
deba@2040:       }
deba@2040: 
deba@2040:       if (success) {
deba@2040: 
deba@2040:         typename Graph::template ANodeMap<bool> aused(*graph, false);
deba@2040:         
deba@2040:         for (int i = 0; i < (int)bqueue.size(); ++i) {
deba@2040:           Node bnode = bqueue[i];
deba@2040: 
deba@2040:           bool used = false;
deba@2040: 
deba@2040:           while (bnode != INVALID) {
deba@2040:             UEdge uedge = bpred[bnode];
deba@2040:             Node anode = graph->aNode(uedge);
deba@2040:             
deba@2040:             if (aused[anode]) {
deba@2040:               used = true;
deba@2040:               break;
deba@2040:             }
deba@2040:             
deba@2040:             bnode = anode_matching[anode] != INVALID ? 
deba@2040:               graph->bNode(anode_matching[anode]) : INVALID;
deba@2040:             
deba@2040:           }
deba@2040:           
deba@2040:           if (used) continue;
deba@2040: 
deba@2040:           bnode = bqueue[i];
deba@2040:           while (bnode != INVALID) {
deba@2040:             UEdge uedge = bpred[bnode];
deba@2040:             Node anode = graph->aNode(uedge);
deba@2040:             
deba@2040:             bnode_matching[bnode] = uedge;
deba@2040: 
deba@2040:             bnode = anode_matching[anode] != INVALID ? 
deba@2040:               graph->bNode(anode_matching[anode]) : INVALID;
deba@2040:             
deba@2040:             anode_matching[anode] = uedge;
deba@2040: 
deba@2040:             aused[anode] = true;
deba@2040:           }
deba@2051:           ++matching_size;
deba@2040: 
deba@2040:         }
deba@2040:       }
deba@2040:       return success;
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief An augmenting phase of the Ford-Fulkerson algorithm
deba@2040:     ///
deba@2040:     /// It runs an augmenting phase of the Ford-Fulkerson
deba@2040:     /// algorithm. The phase finds only one augmenting path and 
deba@2040:     /// augments only on this paths. The algorithm consists at most 
deba@2040:     /// of \f$ O(n) \f$ simple phase and one phase is at most 
deba@2040:     /// \f$ O(e) \f$ long.
deba@2040:     bool simpleAugment() {
deba@2040: 
deba@2040:       typename Graph::template ANodeMap<bool> areached(*graph, false);
deba@2040:       typename Graph::template BNodeMap<bool> breached(*graph, false);
deba@2040: 
deba@2040:       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
deba@2040: 
deba@2040:       std::vector<Node> queue;
deba@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         if (anode_matching[it] == INVALID) {
deba@2040:           queue.push_back(it);
deba@2040:           areached[it] = true;
deba@2040:         }
deba@2040:       }
deba@2040: 
deba@2040:       while (!queue.empty()) {
deba@2040:         std::vector<Node> newqueue;
deba@2040:         for (int i = 0; i < (int)queue.size(); ++i) {
deba@2040:           Node anode = queue[i];
deba@2040:           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
deba@2040:             Node bnode = graph->bNode(jt);
deba@2040:             if (breached[bnode]) continue;
deba@2040:             breached[bnode] = true;
deba@2040:             bpred[bnode] = jt;
deba@2040:             if (bnode_matching[bnode] == INVALID) {
deba@2040:               while (bnode != INVALID) {
deba@2040:                 UEdge uedge = bpred[bnode];
deba@2040:                 anode = graph->aNode(uedge);
deba@2040:                 
deba@2040:                 bnode_matching[bnode] = uedge;
deba@2040:                 
deba@2040:                 bnode = anode_matching[anode] != INVALID ? 
deba@2040:                   graph->bNode(anode_matching[anode]) : INVALID;
deba@2040:                 
deba@2040:                 anode_matching[anode] = uedge;
deba@2040:                 
deba@2040:               }
deba@2051:               ++matching_size;
deba@2040:               return true;
deba@2040:             } else {           
deba@2040:               Node newanode = graph->aNode(bnode_matching[bnode]);
deba@2040:               if (!areached[newanode]) {
deba@2040:                 areached[newanode] = true;
deba@2040:                 newqueue.push_back(newanode);
deba@2040:               }
deba@2040:             }
deba@2040:           }
deba@2040:         }
deba@2040:         queue.swap(newqueue);
deba@2040:       }
deba@2040:       
deba@2040:       return false;
deba@2040:     }
deba@2040: 
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@2040:     /// \brief Returns an minimum covering of the nodes.
deba@2040:     ///
deba@2040:     /// The minimum covering set problem is the dual solution of the
deba@2040:     /// maximum bipartite matching. It provides an 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@2040:       typename Graph::template ANodeMap<bool> areached(*graph, false);
deba@2040:       typename Graph::template BNodeMap<bool> breached(*graph, false);
deba@2040:       
deba@2040:       std::vector<Node> queue;
deba@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         if (anode_matching[it] == INVALID) {
deba@2040:           queue.push_back(it);
deba@2040:         }
deba@2040:       }
deba@2040: 
deba@2040:       while (!queue.empty()) {
deba@2040:         std::vector<Node> newqueue;
deba@2040:         for (int i = 0; i < (int)queue.size(); ++i) {
deba@2040:           Node anode = queue[i];
deba@2040:           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
deba@2040:             Node bnode = graph->bNode(jt);
deba@2040:             if (breached[bnode]) continue;
deba@2040:             breached[bnode] = true;
deba@2040:             if (bnode_matching[bnode] != INVALID) {
deba@2040:               Node newanode = graph->aNode(bnode_matching[bnode]);
deba@2040:               if (!areached[newanode]) {
deba@2040:                 areached[newanode] = true;
deba@2040:                 newqueue.push_back(newanode);
deba@2040:               }
deba@2040:             }
deba@2040:           }
deba@2040:         }
deba@2040:         queue.swap(newqueue);
deba@2040:       }
deba@2040: 
deba@2040:       int size = 0;
deba@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         covering[it] = !areached[it] && anode_matching[it] != INVALID;
deba@2040:         if (!areached[it] && anode_matching[it] != INVALID) {
deba@2040:           ++size;
deba@2040:         }
deba@2040:       }
deba@2040:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         covering[it] = breached[it];
deba@2040:         if (breached[it]) {
deba@2040:           ++size;
deba@2040:         }
deba@2040:       }
deba@2040:       return size;
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief Set true all matching uedge in the map.
deba@2040:     /// 
deba@2040:     /// Set true all matching uedge in the map. It does not change the
deba@2040:     /// value mapped to the other uedges.
deba@2040:     /// \return The number of the matching edges.
deba@2040:     template <typename MatchingMap>
deba@2058:     int quickMatching(MatchingMap& matching) const {
deba@2040:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2040:         if (anode_matching[it] != INVALID) {
deba@2040:           matching[anode_matching[it]] = true;
deba@2040:         }
deba@2040:       }
deba@2051:       return matching_size;
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief Set true all matching uedge in the map and the others to false.
deba@2040:     /// 
deba@2040:     /// Set true all matching uedge in the map and the others to false.
deba@2040:     /// \return The number of the matching edges.
deba@2040:     template <typename MatchingMap>
deba@2058:     int matching(MatchingMap& matching) const {
deba@2040:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2040:         matching[it] = it == anode_matching[graph->aNode(it)];
deba@2040:       }
deba@2051:       return matching_size;
deba@2040:     }
deba@2040: 
deba@2040: 
deba@2040:     /// \brief Return true if the given uedge is in the matching.
deba@2040:     /// 
deba@2040:     /// It returns true if the given uedge is in the matching.
deba@2058:     bool matchingEdge(const UEdge& edge) const {
deba@2040:       return anode_matching[graph->aNode(edge)] == edge;
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief Returns the matching edge from the node.
deba@2040:     /// 
deba@2040:     /// Returns the matching edge from the node. If there is not such
deba@2040:     /// edge it gives back \c INVALID.
deba@2058:     UEdge matchingEdge(const Node& node) const {
deba@2040:       if (graph->aNode(node)) {
deba@2040:         return anode_matching[node];
deba@2040:       } else {
deba@2040:         return bnode_matching[node];
deba@2040:       }
deba@2040:     }
deba@2040: 
deba@2040:     /// \brief Gives back the number of the matching edges.
deba@2040:     ///
deba@2040:     /// Gives back the number of the matching edges.
deba@2051:     int matchingSize() const {
deba@2051:       return matching_size;
deba@2040:     }
deba@2040: 
deba@2040:     /// @}
deba@2040: 
deba@2040:   private:
deba@2040: 
deba@2040:     ANodeMatchingMap anode_matching;
deba@2040:     BNodeMatchingMap bnode_matching;
deba@2040:     const Graph *graph;
deba@2040: 
deba@2051:     int matching_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@2058:   /// \retval matching The undirected edge map which will be set to 
deba@2058:   /// the matching.
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@2058:     bpmatching.matching(matching);
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@2051:   struct WeightedBipartiteMatchingDefaultTraits {
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@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 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
deba@2051:     typedef BinHeap<typename BpUGraph::Node, 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@2051:             typename _Traits = WeightedBipartiteMatchingDefaultTraits<_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@2051:     MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
deba@2051:       if(local_heap_cross_ref) {
deba@2051: 	delete _heap_cross_ref;
deba@2051: 	local_heap_cross_ref = false;
deba@2051:       }
deba@2051:       _heap_cross_ref = &crossRef;
deba@2051:       if(local_heap) {
deba@2051: 	delete _heap;
deba@2051: 	local_heap = false;
deba@2051:       }
deba@2051:       _heap = &heap;
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 
deba@2051:     /// algorithm. The 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@2051:           if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
deba@2051:             bdist[bnode] = bvalue;
deba@2051:             bpred[bnode] = jt;
deba@2051:           }
deba@2051:           if (bvalue > bdistMax) {
deba@2051:             bdistMax = bvalue;
deba@2051:           }
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@2058:     void potential(PotentialMap& potential) const {
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         potential[it] = anode_potential[it];
deba@2051:       }
deba@2051:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         potential[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@2058:     int quickMatching(MatchingMap& matching) const {
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (anode_matching[it] != INVALID) {
deba@2051:           matching[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@2058:     int matching(MatchingMap& matching) const {
deba@2051:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2051:         matching[it] = it == anode_matching[graph->aNode(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@2051:     typedef BinHeap<typename BpUGraph::Node, 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@2051:             typename _Traits = 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@2051:     MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
deba@2051:       if(local_heap_cross_ref) {
deba@2051: 	delete _heap_cross_ref;
deba@2051: 	local_heap_cross_ref = false;
deba@2051:       }
deba@2051:       _heap_cross_ref = &crossRef;
deba@2051:       if(local_heap) {
deba@2051: 	delete _heap;
deba@2051: 	local_heap = false;
deba@2051:       }
deba@2051:       _heap = &heap;
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@2051:           if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
deba@2051:             bdist[bnode] = bvalue;
deba@2051:             bpred[bnode] = jt;
deba@2051:           }
deba@2136:           if (bvalue > bdistMax) {
deba@2136:             bdistMax = bvalue;
deba@2136:           }
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@2058:     /// Gives back the potential in the NodeMap. The potential is optimal with 
deba@2058:     /// the current number of edges if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
deba@2051:     /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
deba@2051:     /// for each edges.
deba@2051:     template <typename PotentialMap>
deba@2058:     void potential(PotentialMap& potential) const {
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         potential[it] = anode_potential[it];
deba@2051:       }
deba@2051:       for (BNodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         potential[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@2058:     int quickMatching(MatchingMap& matching) const {
deba@2051:       for (ANodeIt it(*graph); it != INVALID; ++it) {
deba@2051:         if (anode_matching[it] != INVALID) {
deba@2051:           matching[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@2058:     int matching(MatchingMap& matching) const {
deba@2051:       for (UEdgeIt it(*graph); it != INVALID; ++it) {
deba@2051:         matching[it] = it == anode_matching[graph->aNode(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 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@2058:   /// This function calculates the minimum cost matching of the maximum 
deba@2058:   /// cardinality matchings of a bipartite graph. It gives back the matching 
deba@2058:   /// in 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