# HG changeset patch
# User deba
# Date 1145038053 0
# Node ID 08652c1763f6566080c711781141127d8ea1ea51
# Parent  d9a221218ea4a6f203e1250f1fe8c772de84a7ec
MaxWeightedBipartiteMatching
MinCostMaxBipartiteMatching

Both algorithms are based on successive shortest
path algorithm with dijkstra shortest path
finding

diff -r d9a221218ea4 -r 08652c1763f6 lemon/bipartite_matching.h
--- a/lemon/bipartite_matching.h	Fri Apr 14 18:05:02 2006 +0000
+++ b/lemon/bipartite_matching.h	Fri Apr 14 18:07:33 2006 +0000
@@ -19,8 +19,10 @@
 #ifndef LEMON_BIPARTITE_MATCHING
 #define LEMON_BIPARTITE_MATCHING
 
-#include <lemon/bpugraph_adaptor.h>
-#include <lemon/bfs.h>
+#include <functional>
+
+#include <lemon/bin_heap.h>
+#include <lemon/maps.h>
 
 #include <iostream>
 
@@ -35,7 +37,7 @@
   /// \brief Bipartite Max Cardinality Matching algorithm
   ///
   /// Bipartite Max Cardinality Matching algorithm. This class implements
-  /// the Hopcroft-Karp algorithm wich has \f$ O(e\sqrt{n}) \f$ time
+  /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
   /// complexity.
   template <typename BpUGraph>
   class MaxBipartiteMatching {
@@ -83,7 +85,7 @@
       for (BNodeIt it(*graph); it != INVALID; ++it) {
         bnode_matching[it] = INVALID;
       }
-      matching_value = 0;
+      matching_size = 0;
     }
 
     /// \brief Initalize the data structures.
@@ -92,7 +94,7 @@
     /// matching.  From this matching sometimes it is faster to get
     /// the matching than from the initial empty matching.
     void greedyInit() {
-      matching_value = 0;
+      matching_size = 0;
       for (BNodeIt it(*graph); it != INVALID; ++it) {
         bnode_matching[it] = INVALID;
       }
@@ -102,7 +104,7 @@
           if (bnode_matching[graph->bNode(jt)] == INVALID) {
             anode_matching[it] = jt;
             bnode_matching[graph->bNode(jt)] = jt;
-            ++matching_value;
+            ++matching_size;
             break;
           }
         }
@@ -120,10 +122,10 @@
       for (BNodeIt it(*graph); it != INVALID; ++it) {
         bnode_matching[it] = INVALID;
       }
-      matching_value = 0;
+      matching_size = 0;
       for (UEdgeIt it(*graph); it != INVALID; ++it) {
         if (matching[it]) {
-          ++matching_value;
+          ++matching_size;
           anode_matching[graph->aNode(it)] = it;
           bnode_matching[graph->bNode(it)] = it;
         }
@@ -142,10 +144,10 @@
       for (BNodeIt it(*graph); it != INVALID; ++it) {
         bnode_matching[it] = INVALID;
       }
-      matching_value = 0;
+      matching_size = 0;
       for (UEdgeIt it(*graph); it != INVALID; ++it) {
         if (matching[it]) {
-          ++matching_value;
+          ++matching_size;
           if (anode_matching[graph->aNode(it)] != INVALID) {
             return false;
           }
@@ -246,7 +248,7 @@
 
             aused[anode] = true;
           }
-          ++matching_value;
+          ++matching_size;
 
         }
       }
@@ -297,7 +299,7 @@
                 anode_matching[anode] = uedge;
                 
               }
-              ++matching_value;
+              ++matching_size;
               return true;
             } else {           
               Node newanode = graph->aNode(bnode_matching[bnode]);
@@ -407,7 +409,7 @@
           matching[anode_matching[it]] = true;
         }
       }
-      return matching_value;
+      return matching_size;
     }
 
     /// \brief Set true all matching uedge in the map and the others to false.
@@ -419,7 +421,7 @@
       for (UEdgeIt it(*graph); it != INVALID; ++it) {
         matching[it] = it == anode_matching[graph->aNode(it)];
       }
-      return matching_value;
+      return matching_size;
     }
 
 
@@ -445,8 +447,8 @@
     /// \brief Gives back the number of the matching edges.
     ///
     /// Gives back the number of the matching edges.
-    int matchingValue() const {
-      return matching_value;
+    int matchingSize() const {
+      return matching_size;
     }
 
     /// @}
@@ -457,7 +459,1022 @@
     BNodeMatchingMap bnode_matching;
     const Graph *graph;
 
-    int matching_value;
+    int matching_size;
+  
+  };
+
+  /// \brief Default traits class for weighted bipartite matching algoritms.
+  ///
+  /// Default traits class for weighted bipartite matching algoritms.
+  /// \param _BpUGraph The bipartite undirected graph type.
+  /// \param _WeightMap Type of weight map.
+  template <typename _BpUGraph, typename _WeightMap>
+  struct WeightedBipartiteMatchingDefaultTraits {
+    /// \brief The type of the weight of the undirected edges.
+    typedef typename _WeightMap::Value Value;
+
+    /// The undirected bipartite graph type the algorithm runs on. 
+    typedef _BpUGraph BpUGraph;
+
+    /// The map of the edges weights
+    typedef _WeightMap WeightMap;
+
+    /// \brief The cross reference type used by heap.
+    ///
+    /// The cross reference type used by heap.
+    /// Usually it is \c Graph::NodeMap<int>.
+    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
+
+    /// \brief Instantiates a HeapCrossRef.
+    ///
+    /// This function instantiates a \ref HeapCrossRef. 
+    /// \param graph is the graph, to which we would like to define the 
+    /// HeapCrossRef.
+    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
+      return new HeapCrossRef(graph);
+    }
+    
+    /// \brief The heap type used by weighted matching algorithms.
+    ///
+    /// The heap type used by weighted matching algorithms. It should
+    /// minimize the priorities and the heap's key type is the graph's
+    /// anode graph's node.
+    ///
+    /// \sa BinHeap
+    typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
+    
+    /// \brief Instantiates a Heap.
+    ///
+    /// This function instantiates a \ref Heap. 
+    /// \param crossref The cross reference of the heap.
+    static Heap *createHeap(HeapCrossRef& crossref) {
+      return new Heap(crossref);
+    }
+
+  };
+
+
+  /// \ingroup matching
+  ///
+  /// \brief Bipartite Max Weighted Matching algorithm
+  ///
+  /// This class implements the bipartite Max Weighted Matching
+  /// algorithm.  It uses the successive shortest path algorithm to
+  /// calculate the maximum weighted matching in the bipartite
+  /// graph. The algorithm can be used also to calculate the maximum
+  /// cardinality maximum weighted matching. The time complexity
+  /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
+  /// heap implementation but this can be improved to 
+  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
+  ///
+  /// The algorithm also provides a potential function on the nodes
+  /// which a dual solution of the matching algorithm and it can be
+  /// used to proof the optimality of the given pimal solution.
+#ifdef DOXYGEN
+  template <typename _BpUGraph, typename _WeightMap, typename _Traits>
+#else
+  template <typename _BpUGraph, 
+            typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
+            typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
+#endif
+  class MaxWeightedBipartiteMatching {
+  public:
+
+    typedef _Traits Traits;
+    typedef typename Traits::BpUGraph BpUGraph;
+    typedef typename Traits::WeightMap WeightMap;
+    typedef typename Traits::Value Value;
+
+  protected:
+
+    typedef typename Traits::HeapCrossRef HeapCrossRef;
+    typedef typename Traits::Heap Heap; 
+
+    
+    typedef typename BpUGraph::Node Node;
+    typedef typename BpUGraph::ANodeIt ANodeIt;
+    typedef typename BpUGraph::BNodeIt BNodeIt;
+    typedef typename BpUGraph::UEdge UEdge;
+    typedef typename BpUGraph::UEdgeIt UEdgeIt;
+    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
+
+    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
+    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
+
+    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
+    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
+
+
+  public:
+
+    /// \brief \ref Exception for uninitialized parameters.
+    ///
+    /// This error represents problems in the initialization
+    /// of the parameters of the algorithms.
+    class UninitializedParameter : public lemon::UninitializedParameter {
+    public:
+      virtual const char* exceptionName() const {
+	return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
+      }
+    };
+
+    ///\name Named template parameters
+
+    ///@{
+
+    template <class H, class CR>
+    struct DefHeapTraits : public Traits {
+      typedef CR HeapCrossRef;
+      typedef H Heap;
+      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
+	throw UninitializedParameter();
+      }
+      static Heap *createHeap(HeapCrossRef &) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting heap 
+    /// and cross reference type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting heap and cross 
+    /// reference type
+    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
+    struct DefHeap
+      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
+                                            DefHeapTraits<H, CR> > { 
+      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
+                                           DefHeapTraits<H, CR> > Create;
+    };
+
+    template <class H, class CR>
+    struct DefStandardHeapTraits : public Traits {
+      typedef CR HeapCrossRef;
+      typedef H Heap;
+      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
+	return new HeapCrossRef(graph);
+      }
+      static Heap *createHeap(HeapCrossRef &crossref) {
+	return new Heap(crossref);
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting heap and 
+    /// cross reference type with automatic allocation
+    ///
+    /// \ref named-templ-param "Named parameter" for setting heap and cross 
+    /// reference type. It can allocate the heap and the cross reference 
+    /// object if the cross reference's constructor waits for the graph as 
+    /// parameter and the heap's constructor waits for the cross reference.
+    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
+    struct DefStandardHeap
+      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
+                                            DefStandardHeapTraits<H, CR> > { 
+      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
+                                           DefStandardHeapTraits<H, CR> > 
+      Create;
+    };
+
+    ///@}
+
+
+    /// \brief Constructor.
+    ///
+    /// Constructor of the algorithm. 
+    MaxWeightedBipartiteMatching(const BpUGraph& _graph, 
+                                 const WeightMap& _weight) 
+      : graph(&_graph), weight(&_weight),
+        anode_matching(_graph), bnode_matching(_graph),
+        anode_potential(_graph), bnode_potential(_graph),
+        _heap_cross_ref(0), local_heap_cross_ref(false),
+        _heap(0), local_heap(0) {}
+
+    /// \brief Destructor.
+    ///
+    /// Destructor of the algorithm.
+    ~MaxWeightedBipartiteMatching() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the heap and the cross reference used by algorithm.
+    ///
+    /// Sets the heap and the cross reference used by algorithm.
+    /// If you don't use this function before calling \ref run(),
+    /// it will allocate one. The destuctor deallocates this
+    /// automatically allocated map, of course.
+    /// \return \c (*this)
+    MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
+      if(local_heap_cross_ref) {
+	delete _heap_cross_ref;
+	local_heap_cross_ref = false;
+      }
+      _heap_cross_ref = &crossRef;
+      if(local_heap) {
+	delete _heap;
+	local_heap = false;
+      }
+      _heap = &heap;
+      return *this;
+    }
+
+    /// \name Execution control
+    /// The simplest way to execute the algorithm is to use
+    /// one of the member functions called \c run().
+    /// \n
+    /// If you need more control on the execution,
+    /// first you must call \ref init() or one alternative for it.
+    /// Finally \ref start() will perform the matching computation or
+    /// with step-by-step execution you can augment the solution.
+
+    /// @{
+
+    /// \brief Initalize the data structures.
+    ///
+    /// It initalizes the data structures and creates an empty matching.
+    void init() {
+      initStructures();
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        anode_matching[it] = INVALID;
+        anode_potential[it] = 0;
+      }
+      for (BNodeIt it(*graph); it != INVALID; ++it) {
+        bnode_matching[it] = INVALID;
+        bnode_potential[it] = 0;
+        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
+          if (-(*weight)[jt] <= bnode_potential[it]) {
+            bnode_potential[it] = -(*weight)[jt];
+          }
+        }
+      }
+      matching_value = 0;
+      matching_size = 0;
+    }
+
+
+    /// \brief An augmenting phase of the weighted matching algorithm
+    ///
+    /// It runs an augmenting phase of the weighted matching 
+    /// algorithm. The phase finds the best augmenting path and 
+    /// augments only on this paths. 
+    ///
+    /// The algorithm consists at most 
+    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 
+    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 
+    /// with binary heap.
+    /// \param decrease If the given parameter true the matching value
+    /// can be decreased in the augmenting phase. If we would like
+    /// to calculate the maximum cardinality maximum weighted matching
+    /// then we should let the algorithm to decrease the matching
+    /// value in order to increase the number of the matching edges.
+    bool augment(bool decrease = false) {
+
+      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
+      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
+
+      Node bestNode = INVALID;
+      Value bestValue = 0;
+
+      _heap->clear();
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
+      }
+
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        if (anode_matching[it] == INVALID) {
+          _heap->push(it, 0);
+        }
+      }
+
+      Value bdistMax = 0;
+      while (!_heap->empty()) {
+        Node anode = _heap->top();
+        Value avalue = _heap->prio();
+        _heap->pop();
+        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
+          if (jt == anode_matching[anode]) continue;
+          Node bnode = graph->bNode(jt);
+          Value bvalue = avalue + anode_potential[anode] - 
+            bnode_potential[bnode] - (*weight)[jt];
+          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
+            bdist[bnode] = bvalue;
+            bpred[bnode] = jt;
+          }
+          if (bvalue > bdistMax) {
+            bdistMax = bvalue;
+          }
+          if (bnode_matching[bnode] != INVALID) {
+            Node newanode = graph->aNode(bnode_matching[bnode]);
+            switch (_heap->state(newanode)) {
+            case Heap::PRE_HEAP:
+              _heap->push(newanode, bvalue);
+              break;
+            case Heap::IN_HEAP:
+              if (bvalue < (*_heap)[newanode]) {
+                _heap->decrease(newanode, bvalue);
+              }
+              break;
+            case Heap::POST_HEAP:
+              break;
+            }
+          } else {
+            if (bestNode == INVALID || 
+                - bvalue - bnode_potential[bnode] > bestValue) {
+              bestValue = - bvalue - bnode_potential[bnode];
+              bestNode = bnode;
+            }
+          }
+        }
+      }
+
+      if (bestNode == INVALID || (!decrease && bestValue < 0)) {
+        return false;
+      }
+
+      matching_value += bestValue;
+      ++matching_size;
+
+      for (BNodeIt it(*graph); it != INVALID; ++it) {
+        if (bpred[it] != INVALID) {
+          bnode_potential[it] += bdist[it];
+        } else {
+          bnode_potential[it] += bdistMax;
+        }
+      }
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        if (anode_matching[it] != INVALID) {
+          Node bnode = graph->bNode(anode_matching[it]);
+          if (bpred[bnode] != INVALID) {
+            anode_potential[it] += bdist[bnode];
+          } else {
+            anode_potential[it] += bdistMax;
+          }
+        }
+      }
+
+      while (bestNode != INVALID) {
+        UEdge uedge = bpred[bestNode];
+        Node anode = graph->aNode(uedge);
+        
+        bnode_matching[bestNode] = uedge;
+        if (anode_matching[anode] != INVALID) {
+          bestNode = graph->bNode(anode_matching[anode]);
+        } else {
+          bestNode = INVALID;
+        }
+        anode_matching[anode] = uedge;
+      }
+
+
+      return true;
+    }
+
+    /// \brief Starts the algorithm.
+    ///
+    /// Starts the algorithm. It runs augmenting phases until the
+    /// optimal solution reached.
+    ///
+    /// \param maxCardinality If the given value is true it will
+    /// calculate the maximum cardinality maximum matching instead of
+    /// the maximum matching.
+    void start(bool maxCardinality = false) {
+      while (augment(maxCardinality)) {}
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// It just initalize the algorithm and then start it.
+    ///
+    /// \param maxCardinality If the given value is true it will
+    /// calculate the maximum cardinality maximum matching instead of
+    /// the maximum matching.
+    void run(bool maxCardinality = false) {
+      init();
+      start(maxCardinality);
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %Matching algorithm can be obtained using these
+    /// functions.\n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+
+    /// \brief Gives back the potential in the NodeMap
+    ///
+    /// Gives back the potential in the NodeMap. The potential
+    /// is feasible if \f$ \pi(a) - \pi(b) - w(ab) = 0 \f$ for
+    /// each matching edges and \f$ \pi(a) - \pi(b) - w(ab) \ge 0 \f$
+    /// for each edges.
+    template <typename PotentialMap>
+    void potential(PotentialMap& potential) {
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        potential[it] = anode_potential[it];
+      }
+      for (BNodeIt it(*graph); it != INVALID; ++it) {
+        potential[it] = bnode_potential[it];
+      }
+    }
+
+    /// \brief Set true all matching uedge in the map.
+    /// 
+    /// Set true all matching uedge in the map. It does not change the
+    /// value mapped to the other uedges.
+    /// \return The number of the matching edges.
+    template <typename MatchingMap>
+    int quickMatching(MatchingMap& matching) {
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        if (anode_matching[it] != INVALID) {
+          matching[anode_matching[it]] = true;
+        }
+      }
+      return matching_size;
+    }
+
+    /// \brief Set true all matching uedge in the map and the others to false.
+    /// 
+    /// Set true all matching uedge in the map and the others to false.
+    /// \return The number of the matching edges.
+    template <typename MatchingMap>
+    int matching(MatchingMap& matching) {
+      for (UEdgeIt it(*graph); it != INVALID; ++it) {
+        matching[it] = it == anode_matching[graph->aNode(it)];
+      }
+      return matching_size;
+    }
+
+
+    /// \brief Return true if the given uedge is in the matching.
+    /// 
+    /// It returns true if the given uedge is in the matching.
+    bool matchingEdge(const UEdge& edge) {
+      return anode_matching[graph->aNode(edge)] == edge;
+    }
+
+    /// \brief Returns the matching edge from the node.
+    /// 
+    /// Returns the matching edge from the node. If there is not such
+    /// edge it gives back \c INVALID.
+    UEdge matchingEdge(const Node& node) {
+      if (graph->aNode(node)) {
+        return anode_matching[node];
+      } else {
+        return bnode_matching[node];
+      }
+    }
+
+    /// \brief Gives back the sum of weights of the matching edges.
+    ///
+    /// Gives back the sum of weights of the matching edges.
+    Value matchingValue() const {
+      return matching_value;
+    }
+
+    /// \brief Gives back the number of the matching edges.
+    ///
+    /// Gives back the number of the matching edges.
+    int matchingSize() const {
+      return matching_size;
+    }
+
+    /// @}
+
+  private:
+
+    void initStructures() {
+      if (!_heap_cross_ref) {
+	local_heap_cross_ref = true;
+	_heap_cross_ref = Traits::createHeapCrossRef(*graph);
+      }
+      if (!_heap) {
+	local_heap = true;
+	_heap = Traits::createHeap(*_heap_cross_ref);
+      }
+    }
+
+    void destroyStructures() {
+      if (local_heap_cross_ref) delete _heap_cross_ref;
+      if (local_heap) delete _heap;
+    }
+
+
+  private:
+    
+    const BpUGraph *graph;
+    const WeightMap* weight;
+
+    ANodeMatchingMap anode_matching;
+    BNodeMatchingMap bnode_matching;
+
+    ANodePotentialMap anode_potential;
+    BNodePotentialMap bnode_potential;
+
+    Value matching_value;
+    int matching_size;
+
+    HeapCrossRef *_heap_cross_ref;
+    bool local_heap_cross_ref;
+
+    Heap *_heap;
+    bool local_heap;
+  
+  };
+
+  /// \brief Default traits class for minimum cost bipartite matching
+  /// algoritms.
+  ///
+  /// Default traits class for minimum cost bipartite matching
+  /// algoritms.  
+  ///
+  /// \param _BpUGraph The bipartite undirected graph
+  /// type.  
+  ///
+  /// \param _CostMap Type of cost map.
+  template <typename _BpUGraph, typename _CostMap>
+  struct MinCostMaxBipartiteMatchingDefaultTraits {
+    /// \brief The type of the cost of the undirected edges.
+    typedef typename _CostMap::Value Value;
+
+    /// The undirected bipartite graph type the algorithm runs on. 
+    typedef _BpUGraph BpUGraph;
+
+    /// The map of the edges costs
+    typedef _CostMap CostMap;
+
+    /// \brief The cross reference type used by heap.
+    ///
+    /// The cross reference type used by heap.
+    /// Usually it is \c Graph::NodeMap<int>.
+    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
+
+    /// \brief Instantiates a HeapCrossRef.
+    ///
+    /// This function instantiates a \ref HeapCrossRef. 
+    /// \param graph is the graph, to which we would like to define the 
+    /// HeapCrossRef.
+    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
+      return new HeapCrossRef(graph);
+    }
+    
+    /// \brief The heap type used by costed matching algorithms.
+    ///
+    /// The heap type used by costed matching algorithms. It should
+    /// minimize the priorities and the heap's key type is the graph's
+    /// anode graph's node.
+    ///
+    /// \sa BinHeap
+    typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
+    
+    /// \brief Instantiates a Heap.
+    ///
+    /// This function instantiates a \ref Heap. 
+    /// \param crossref The cross reference of the heap.
+    static Heap *createHeap(HeapCrossRef& crossref) {
+      return new Heap(crossref);
+    }
+
+  };
+
+
+  /// \ingroup matching
+  ///
+  /// \brief Bipartite Min Cost Matching algorithm
+  ///
+  /// This class implements the bipartite Min Cost Matching algorithm.
+  /// It uses the successive shortest path algorithm to calculate the
+  /// minimum cost maximum matching in the bipartite graph. The time
+  /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
+  /// default binary heap implementation but this can be improved to
+  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
+  ///
+  /// The algorithm also provides a potential function on the nodes
+  /// which a dual solution of the matching algorithm and it can be
+  /// used to proof the optimality of the given pimal solution.
+#ifdef DOXYGEN
+  template <typename _BpUGraph, typename _CostMap, typename _Traits>
+#else
+  template <typename _BpUGraph, 
+            typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
+            typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
+#endif
+  class MinCostMaxBipartiteMatching {
+  public:
+
+    typedef _Traits Traits;
+    typedef typename Traits::BpUGraph BpUGraph;
+    typedef typename Traits::CostMap CostMap;
+    typedef typename Traits::Value Value;
+
+  protected:
+
+    typedef typename Traits::HeapCrossRef HeapCrossRef;
+    typedef typename Traits::Heap Heap; 
+
+    
+    typedef typename BpUGraph::Node Node;
+    typedef typename BpUGraph::ANodeIt ANodeIt;
+    typedef typename BpUGraph::BNodeIt BNodeIt;
+    typedef typename BpUGraph::UEdge UEdge;
+    typedef typename BpUGraph::UEdgeIt UEdgeIt;
+    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
+
+    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
+    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
+
+    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
+    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
+
+
+  public:
+
+    /// \brief \ref Exception for uninitialized parameters.
+    ///
+    /// This error represents problems in the initialization
+    /// of the parameters of the algorithms.
+    class UninitializedParameter : public lemon::UninitializedParameter {
+    public:
+      virtual const char* exceptionName() const {
+	return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
+      }
+    };
+
+    ///\name Named template parameters
+
+    ///@{
+
+    template <class H, class CR>
+    struct DefHeapTraits : public Traits {
+      typedef CR HeapCrossRef;
+      typedef H Heap;
+      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
+	throw UninitializedParameter();
+      }
+      static Heap *createHeap(HeapCrossRef &) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting heap 
+    /// and cross reference type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting heap and cross 
+    /// reference type
+    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
+    struct DefHeap
+      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
+                                            DefHeapTraits<H, CR> > { 
+      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
+                                           DefHeapTraits<H, CR> > Create;
+    };
+
+    template <class H, class CR>
+    struct DefStandardHeapTraits : public Traits {
+      typedef CR HeapCrossRef;
+      typedef H Heap;
+      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
+	return new HeapCrossRef(graph);
+      }
+      static Heap *createHeap(HeapCrossRef &crossref) {
+	return new Heap(crossref);
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting heap and 
+    /// cross reference type with automatic allocation
+    ///
+    /// \ref named-templ-param "Named parameter" for setting heap and cross 
+    /// reference type. It can allocate the heap and the cross reference 
+    /// object if the cross reference's constructor waits for the graph as 
+    /// parameter and the heap's constructor waits for the cross reference.
+    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
+    struct DefStandardHeap
+      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
+                                            DefStandardHeapTraits<H, CR> > { 
+      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
+                                           DefStandardHeapTraits<H, CR> > 
+      Create;
+    };
+
+    ///@}
+
+
+    /// \brief Constructor.
+    ///
+    /// Constructor of the algorithm. 
+    MinCostMaxBipartiteMatching(const BpUGraph& _graph, 
+                                 const CostMap& _cost) 
+      : graph(&_graph), cost(&_cost),
+        anode_matching(_graph), bnode_matching(_graph),
+        anode_potential(_graph), bnode_potential(_graph),
+        _heap_cross_ref(0), local_heap_cross_ref(false),
+        _heap(0), local_heap(0) {}
+
+    /// \brief Destructor.
+    ///
+    /// Destructor of the algorithm.
+    ~MinCostMaxBipartiteMatching() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the heap and the cross reference used by algorithm.
+    ///
+    /// Sets the heap and the cross reference used by algorithm.
+    /// If you don't use this function before calling \ref run(),
+    /// it will allocate one. The destuctor deallocates this
+    /// automatically allocated map, of course.
+    /// \return \c (*this)
+    MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
+      if(local_heap_cross_ref) {
+	delete _heap_cross_ref;
+	local_heap_cross_ref = false;
+      }
+      _heap_cross_ref = &crossRef;
+      if(local_heap) {
+	delete _heap;
+	local_heap = false;
+      }
+      _heap = &heap;
+      return *this;
+    }
+
+    /// \name Execution control
+    /// The simplest way to execute the algorithm is to use
+    /// one of the member functions called \c run().
+    /// \n
+    /// If you need more control on the execution,
+    /// first you must call \ref init() or one alternative for it.
+    /// Finally \ref start() will perform the matching computation or
+    /// with step-by-step execution you can augment the solution.
+
+    /// @{
+
+    /// \brief Initalize the data structures.
+    ///
+    /// It initalizes the data structures and creates an empty matching.
+    void init() {
+      initStructures();
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        anode_matching[it] = INVALID;
+        anode_potential[it] = 0;
+      }
+      for (BNodeIt it(*graph); it != INVALID; ++it) {
+        bnode_matching[it] = INVALID;
+        bnode_potential[it] = 0;
+      }
+      matching_cost = 0;
+      matching_size = 0;
+    }
+
+
+    /// \brief An augmenting phase of the costed matching algorithm
+    ///
+    /// It runs an augmenting phase of the matching algorithm. The
+    /// phase finds the best augmenting path and augments only on this
+    /// paths.
+    ///
+    /// The algorithm consists at most 
+    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 
+    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 
+    /// with binary heap.
+    bool augment() {
+
+      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
+      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
+
+      Node bestNode = INVALID;
+      Value bestValue = 0;
+
+      _heap->clear();
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
+      }
+
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        if (anode_matching[it] == INVALID) {
+          _heap->push(it, 0);
+        }
+      }
+
+      while (!_heap->empty()) {
+        Node anode = _heap->top();
+        Value avalue = _heap->prio();
+        _heap->pop();
+        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
+          if (jt == anode_matching[anode]) continue;
+          Node bnode = graph->bNode(jt);
+          Value bvalue = avalue + (*cost)[jt] + 
+            anode_potential[anode] - bnode_potential[bnode];
+          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
+            bdist[bnode] = bvalue;
+            bpred[bnode] = jt;
+          }
+          if (bnode_matching[bnode] != INVALID) {
+            Node newanode = graph->aNode(bnode_matching[bnode]);
+            switch (_heap->state(newanode)) {
+            case Heap::PRE_HEAP:
+              _heap->push(newanode, bvalue);
+              break;
+            case Heap::IN_HEAP:
+              if (bvalue < (*_heap)[newanode]) {
+                _heap->decrease(newanode, bvalue);
+              }
+              break;
+            case Heap::POST_HEAP:
+              break;
+            }
+          } else {
+            if (bestNode == INVALID || 
+                bvalue + bnode_potential[bnode] < bestValue) {
+              bestValue = bvalue + bnode_potential[bnode];
+              bestNode = bnode;
+            }
+          }
+        }
+      }
+
+      if (bestNode == INVALID) {
+        return false;
+      }
+
+      matching_cost += bestValue;
+      ++matching_size;
+
+      for (BNodeIt it(*graph); it != INVALID; ++it) {
+        if (bpred[it] != INVALID) {
+          bnode_potential[it] += bdist[it];
+        }
+      }
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        if (anode_matching[it] != INVALID) {
+          Node bnode = graph->bNode(anode_matching[it]);
+          if (bpred[bnode] != INVALID) {
+            anode_potential[it] += bdist[bnode];
+          }
+        }
+      }
+
+      while (bestNode != INVALID) {
+        UEdge uedge = bpred[bestNode];
+        Node anode = graph->aNode(uedge);
+        
+        bnode_matching[bestNode] = uedge;
+        if (anode_matching[anode] != INVALID) {
+          bestNode = graph->bNode(anode_matching[anode]);
+        } else {
+          bestNode = INVALID;
+        }
+        anode_matching[anode] = uedge;
+      }
+
+
+      return true;
+    }
+
+    /// \brief Starts the algorithm.
+    ///
+    /// Starts the algorithm. It runs augmenting phases until the
+    /// optimal solution reached.
+    void start() {
+      while (augment()) {}
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// It just initalize the algorithm and then start it.
+    void run() {
+      init();
+      start();
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %Matching algorithm can be obtained using these
+    /// functions.\n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+
+    /// \brief Gives back the potential in the NodeMap
+    ///
+    /// Gives back the potential in the NodeMap. The potential
+    /// is feasible if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
+    /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
+    /// for each edges.
+    template <typename PotentialMap>
+    void potential(PotentialMap& potential) {
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        potential[it] = anode_potential[it];
+      }
+      for (BNodeIt it(*graph); it != INVALID; ++it) {
+        potential[it] = bnode_potential[it];
+      }
+    }
+
+    /// \brief Set true all matching uedge in the map.
+    /// 
+    /// Set true all matching uedge in the map. It does not change the
+    /// value mapped to the other uedges.
+    /// \return The number of the matching edges.
+    template <typename MatchingMap>
+    int quickMatching(MatchingMap& matching) {
+      for (ANodeIt it(*graph); it != INVALID; ++it) {
+        if (anode_matching[it] != INVALID) {
+          matching[anode_matching[it]] = true;
+        }
+      }
+      return matching_size;
+    }
+
+    /// \brief Set true all matching uedge in the map and the others to false.
+    /// 
+    /// Set true all matching uedge in the map and the others to false.
+    /// \return The number of the matching edges.
+    template <typename MatchingMap>
+    int matching(MatchingMap& matching) {
+      for (UEdgeIt it(*graph); it != INVALID; ++it) {
+        matching[it] = it == anode_matching[graph->aNode(it)];
+      }
+      return matching_size;
+    }
+
+
+    /// \brief Return true if the given uedge is in the matching.
+    /// 
+    /// It returns true if the given uedge is in the matching.
+    bool matchingEdge(const UEdge& edge) {
+      return anode_matching[graph->aNode(edge)] == edge;
+    }
+
+    /// \brief Returns the matching edge from the node.
+    /// 
+    /// Returns the matching edge from the node. If there is not such
+    /// edge it gives back \c INVALID.
+    UEdge matchingEdge(const Node& node) {
+      if (graph->aNode(node)) {
+        return anode_matching[node];
+      } else {
+        return bnode_matching[node];
+      }
+    }
+
+    /// \brief Gives back the sum of costs of the matching edges.
+    ///
+    /// Gives back the sum of costs of the matching edges.
+    Value matchingCost() const {
+      return matching_cost;
+    }
+
+    /// \brief Gives back the number of the matching edges.
+    ///
+    /// Gives back the number of the matching edges.
+    int matchingSize() const {
+      return matching_size;
+    }
+
+    /// @}
+
+  private:
+
+    void initStructures() {
+      if (!_heap_cross_ref) {
+	local_heap_cross_ref = true;
+	_heap_cross_ref = Traits::createHeapCrossRef(*graph);
+      }
+      if (!_heap) {
+	local_heap = true;
+	_heap = Traits::createHeap(*_heap_cross_ref);
+      }
+    }
+
+    void destroyStructures() {
+      if (local_heap_cross_ref) delete _heap_cross_ref;
+      if (local_heap) delete _heap;
+    }
+
+
+  private:
+    
+    const BpUGraph *graph;
+    const CostMap* cost;
+
+    ANodeMatchingMap anode_matching;
+    BNodeMatchingMap bnode_matching;
+
+    ANodePotentialMap anode_potential;
+    BNodePotentialMap bnode_potential;
+
+    Value matching_cost;
+    int matching_size;
+
+    HeapCrossRef *_heap_cross_ref;
+    bool local_heap_cross_ref;
+
+    Heap *_heap;
+    bool local_heap;
   
   };
 
diff -r d9a221218ea4 -r 08652c1763f6 test/bipartite_matching_test.cc
--- a/test/bipartite_matching_test.cc	Fri Apr 14 18:05:02 2006 +0000
+++ b/test/bipartite_matching_test.cc	Fri Apr 14 18:07:33 2006 +0000
@@ -9,6 +9,8 @@
 
 #include <lemon/graph_utils.h>
 
+#include <lemon/maps.h>
+
 #include "test_tools.h"
 
 using namespace std;
@@ -24,6 +26,13 @@
   int n = argc > 1 ? atoi(argv[1]) : 100;
   int m = argc > 2 ? atoi(argv[2]) : 100;
   int e = argc > 3 ? atoi(argv[3]) : (int)((n+m) * log(n+m));
+  int c = argc > 4 ? atoi(argv[4]) : 100;
+
+  Graph::UEdgeMap<int> weight(graph);
+
+  int max_cardinality;
+  int max_weight;
+  int max_cardinality_max_weight;
 
   for (int i = 0; i < n; ++i) {
     Node node = graph.addANode();
@@ -36,7 +45,8 @@
   for (int i = 0; i < e; ++i) {
     Node aNode = aNodes[urandom(n)];
     Node bNode = bNodes[urandom(m)];
-    graph.addEdge(aNode, bNode);
+    UEdge uedge = graph.addEdge(aNode, bNode);
+    weight[uedge] = urandom(c);
   }
 
   {
@@ -57,9 +67,10 @@
       int num = 0;
       for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
         if (mm[jt]) ++num;
-    }
+      }
       check(num <= 1, "INVALID PRIMAL");
     }
+    max_cardinality = bpmatch.matchingSize();
   }
 
   {
@@ -111,20 +122,63 @@
   }
 
   {
-    SwapBpUGraphAdaptor<Graph> swapgraph(graph);
-    MaxBipartiteMatching<SwapBpUGraphAdaptor<Graph> > bpmatch(swapgraph);
+    MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
+
+    bpmatch.init();
+
+    max_weight = 0;
+    while (bpmatch.augment(true)) {
+    
+      Graph::UEdgeMap<bool> mm(graph);
+      Graph::NodeMap<int> pm(graph);
+      
+      bpmatch.matching(mm);
+      bpmatch.potential(pm);
+      
+      for (UEdgeIt it(graph); it != INVALID; ++it) {
+        if (mm[it]) {
+          check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+                "INVALID DUAL");
+        } else {
+          check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+                "INVALID DUAL");
+        }
+      }
+
+      for (ANodeIt it(graph); it != INVALID; ++it) {
+        int num = 0;
+        for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
+          if (mm[jt]) ++num;
+        }
+        check(num <= 1, "INVALID PRIMAL");
+      }
+      if (bpmatch.matchingValue() > max_weight) {
+        max_weight = bpmatch.matchingValue();
+      }
+    }
+  }
+
+  {
+    MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
 
     bpmatch.run();
+    
+    Graph::UEdgeMap<bool> mm(graph);
+    Graph::NodeMap<int> pm(graph);
 
-    Graph::UEdgeMap<bool> mm(graph);
-    Graph::NodeMap<bool> cs(graph);
-    
-    check(bpmatch.coverSet(cs) == bpmatch.matching(mm), "INVALID PRIMAL-DUAL");
+    bpmatch.matching(mm);
+    bpmatch.potential(pm);
     
     for (UEdgeIt it(graph); it != INVALID; ++it) {
-      check(cs[graph.aNode(it)] || cs[graph.bNode(it)], "INVALID DUAL");
+      if (mm[it]) {
+        check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+              "INVALID DUAL");
+      } else {
+        check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+              "INVALID DUAL");
+      }
     }
-    
+
     for (ANodeIt it(graph); it != INVALID; ++it) {
       int num = 0;
       for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
@@ -132,6 +186,79 @@
     }
       check(num <= 1, "INVALID PRIMAL");
     }
+    check(max_weight == bpmatch.matchingValue(), "WRONG WEIGHT");
+  }
+
+  {
+    MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
+
+    bpmatch.run(true);
+    
+    Graph::UEdgeMap<bool> mm(graph);
+    Graph::NodeMap<int> pm(graph);
+
+    bpmatch.matching(mm);
+    bpmatch.potential(pm);
+    
+    for (UEdgeIt it(graph); it != INVALID; ++it) {
+      if (mm[it]) {
+        check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+              "INVALID DUAL");
+      } else {
+        check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+              "INVALID DUAL");
+      }
+    }
+
+    for (ANodeIt it(graph); it != INVALID; ++it) {
+      int num = 0;
+      for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
+        if (mm[jt]) ++num;
+    }
+      check(num <= 1, "INVALID PRIMAL");
+    }
+    check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE");
+    max_cardinality_max_weight = bpmatch.matchingValue();
+
+  }
+
+  {
+    Graph::UEdgeMap<int> cost(graph);
+
+    cost = subMap(constMap<UEdge>(c), weight);
+
+    MinCostMaxBipartiteMatching<Graph> bpmatch(graph, cost);
+
+    bpmatch.run();
+    
+    Graph::UEdgeMap<bool> mm(graph);
+    Graph::NodeMap<int> pm(graph);
+
+    bpmatch.matching(mm);
+    bpmatch.potential(pm);
+    
+    for (UEdgeIt it(graph); it != INVALID; ++it) {
+      if (mm[it]) {
+        check(pm[graph.aNode(it)] + cost[it] - pm[graph.bNode(it)] == 0,
+              "INVALID DUAL");
+      } else {
+        check(pm[graph.aNode(it)] + cost[it] - pm[graph.bNode(it)] >= 0,
+              "INVALID DUAL");
+      }
+    }
+
+    for (ANodeIt it(graph); it != INVALID; ++it) {
+      int num = 0;
+      for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
+        if (mm[jt]) ++num;
+    }
+      check(num <= 1, "INVALID PRIMAL");
+    }
+
+    check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE");
+    check(max_cardinality * c - max_cardinality_max_weight 
+          == bpmatch.matchingCost(), "WRONG SIZE");
+
   }
 
   return 0;