# HG changeset patch
# User deba
# Date 1145343715 0
# Node ID 0b1fc1566fdbf045449fe46829835f68ae6641f8
# Parent  6f84cdb6f4e88d9750e20e42a5f5aa29c2d9941d
Refinements in bipartite matching algorithms

diff -r 6f84cdb6f4e8 -r 0b1fc1566fdb lemon/bipartite_matching.h
--- a/lemon/bipartite_matching.h	Fri Apr 14 23:55:36 2006 +0000
+++ b/lemon/bipartite_matching.h	Tue Apr 18 07:01:55 2006 +0000
@@ -328,7 +328,7 @@
     ///
     /// It just initalize the algorithm and then start it.
     void run() {
-      init();
+      greedyInit();
       start();
     }
 
@@ -349,7 +349,7 @@
     /// problem what is proof of the optimality of the matching.
     /// \return The size of the cover set.
     template <typename CoverMap>
-    int coverSet(CoverMap& covering) {
+    int coverSet(CoverMap& covering) const {
 
       typename Graph::template ANodeMap<bool> areached(*graph, false);
       typename Graph::template BNodeMap<bool> breached(*graph, false);
@@ -403,7 +403,7 @@
     /// value mapped to the other uedges.
     /// \return The number of the matching edges.
     template <typename MatchingMap>
-    int quickMatching(MatchingMap& matching) {
+    int quickMatching(MatchingMap& matching) const {
       for (ANodeIt it(*graph); it != INVALID; ++it) {
         if (anode_matching[it] != INVALID) {
           matching[anode_matching[it]] = true;
@@ -417,7 +417,7 @@
     /// 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) {
+    int matching(MatchingMap& matching) const {
       for (UEdgeIt it(*graph); it != INVALID; ++it) {
         matching[it] = it == anode_matching[graph->aNode(it)];
       }
@@ -428,7 +428,7 @@
     /// \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) {
+    bool matchingEdge(const UEdge& edge) const {
       return anode_matching[graph->aNode(edge)] == edge;
     }
 
@@ -436,7 +436,7 @@
     /// 
     /// Returns the matching edge from the node. If there is not such
     /// edge it gives back \c INVALID.
-    UEdge matchingEdge(const Node& node) {
+    UEdge matchingEdge(const Node& node) const {
       if (graph->aNode(node)) {
         return anode_matching[node];
       } else {
@@ -463,6 +463,26 @@
   
   };
 
+  /// \ingroup matching
+  ///
+  /// \brief Maximum cardinality bipartite matching
+  ///
+  /// This function calculates the maximum cardinality matching
+  /// in a bipartite graph. It gives back the matching in an undirected
+  /// edge map.
+  ///
+  /// \param graph The bipartite graph.
+  /// \retval matching The undirected edge map which will be set to 
+  /// the matching.
+  /// \return The size of the matching.
+  template <typename BpUGraph, typename MatchingMap>
+  int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) {
+    MaxBipartiteMatching<BpUGraph> bpmatching(graph);
+    bpmatching.run();
+    bpmatching.matching(matching);
+    return bpmatching.matchingSize();
+  }
+
   /// \brief Default traits class for weighted bipartite matching algoritms.
   ///
   /// Default traits class for weighted bipartite matching algoritms.
@@ -701,8 +721,8 @@
         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];
+          if ((*weight)[jt] > bnode_potential[it]) {
+            bnode_potential[it] = (*weight)[jt];
           }
         }
       }
@@ -753,8 +773,8 @@
         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];
+          Value bvalue = avalue  - (*weight)[jt] +
+            anode_potential[anode] + bnode_potential[bnode];
           if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
             bdist[bnode] = bvalue;
             bpred[bnode] = jt;
@@ -778,8 +798,8 @@
             }
           } else {
             if (bestNode == INVALID || 
-                - bvalue - bnode_potential[bnode] > bestValue) {
-              bestValue = - bvalue - bnode_potential[bnode];
+                bnode_potential[bnode] - bvalue > bestValue) {
+              bestValue = bnode_potential[bnode] - bvalue;
               bestNode = bnode;
             }
           }
@@ -795,9 +815,9 @@
 
       for (BNodeIt it(*graph); it != INVALID; ++it) {
         if (bpred[it] != INVALID) {
-          bnode_potential[it] += bdist[it];
+          bnode_potential[it] -= bdist[it];
         } else {
-          bnode_potential[it] += bdistMax;
+          bnode_potential[it] -= bdistMax;
         }
       }
       for (ANodeIt it(*graph); it != INVALID; ++it) {
@@ -864,12 +884,12 @@
 
     /// \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.
+    /// Gives back the potential in the NodeMap. The matching is optimal
+    /// with the current number of edges 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) {
+    void potential(PotentialMap& potential) const {
       for (ANodeIt it(*graph); it != INVALID; ++it) {
         potential[it] = anode_potential[it];
       }
@@ -884,7 +904,7 @@
     /// value mapped to the other uedges.
     /// \return The number of the matching edges.
     template <typename MatchingMap>
-    int quickMatching(MatchingMap& matching) {
+    int quickMatching(MatchingMap& matching) const {
       for (ANodeIt it(*graph); it != INVALID; ++it) {
         if (anode_matching[it] != INVALID) {
           matching[anode_matching[it]] = true;
@@ -898,7 +918,7 @@
     /// 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) {
+    int matching(MatchingMap& matching) const {
       for (UEdgeIt it(*graph); it != INVALID; ++it) {
         matching[it] = it == anode_matching[graph->aNode(it)];
       }
@@ -909,7 +929,7 @@
     /// \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) {
+    bool matchingEdge(const UEdge& edge) const {
       return anode_matching[graph->aNode(edge)] == edge;
     }
 
@@ -917,7 +937,7 @@
     /// 
     /// Returns the matching edge from the node. If there is not such
     /// edge it gives back \c INVALID.
-    UEdge matchingEdge(const Node& node) {
+    UEdge matchingEdge(const Node& node) const {
       if (graph->aNode(node)) {
         return anode_matching[node];
       } else {
@@ -982,6 +1002,55 @@
   
   };
 
+  /// \ingroup matching
+  ///
+  /// \brief Maximum weighted bipartite matching
+  ///
+  /// This function calculates the maximum weighted matching
+  /// in a bipartite graph. It gives back the matching in an undirected
+  /// edge map.
+  ///
+  /// \param graph The bipartite graph.
+  /// \param weight The undirected edge map which contains the weights.
+  /// \retval matching The undirected edge map which will be set to 
+  /// the matching.
+  /// \return The value of the matching.
+  template <typename BpUGraph, typename WeightMap, typename MatchingMap>
+  typename WeightMap::Value 
+  maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,
+                               MatchingMap& matching) {
+    MaxWeightedBipartiteMatching<BpUGraph, WeightMap> 
+      bpmatching(graph, weight);
+    bpmatching.run();
+    bpmatching.matching(matching);
+    return bpmatching.matchingValue();
+  }
+
+  /// \ingroup matching
+  ///
+  /// \brief Maximum weighted maximum cardinality bipartite matching
+  ///
+  /// This function calculates the maximum weighted of the maximum cardinality
+  /// matchings of a bipartite graph. It gives back the matching in an 
+  /// undirected edge map.
+  ///
+  /// \param graph The bipartite graph.
+  /// \param weight The undirected edge map which contains the weights.
+  /// \retval matching The undirected edge map which will be set to 
+  /// the matching.
+  /// \return The value of the matching.
+  template <typename BpUGraph, typename WeightMap, typename MatchingMap>
+  typename WeightMap::Value 
+  maxWeightedMaxBipartiteMatching(const BpUGraph& graph, 
+                                  const WeightMap& weight,
+                                  MatchingMap& matching) {
+    MaxWeightedBipartiteMatching<BpUGraph, WeightMap> 
+      bpmatching(graph, weight);
+    bpmatching.run(true);
+    bpmatching.matching(matching);
+    return bpmatching.matchingValue();
+  }
+
   /// \brief Default traits class for minimum cost bipartite matching
   /// algoritms.
   ///
@@ -1360,12 +1429,12 @@
 
     /// \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
+    /// Gives back the potential in the NodeMap. The potential is optimal with 
+    /// the current number of edges 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) {
+    void potential(PotentialMap& potential) const {
       for (ANodeIt it(*graph); it != INVALID; ++it) {
         potential[it] = anode_potential[it];
       }
@@ -1380,7 +1449,7 @@
     /// value mapped to the other uedges.
     /// \return The number of the matching edges.
     template <typename MatchingMap>
-    int quickMatching(MatchingMap& matching) {
+    int quickMatching(MatchingMap& matching) const {
       for (ANodeIt it(*graph); it != INVALID; ++it) {
         if (anode_matching[it] != INVALID) {
           matching[anode_matching[it]] = true;
@@ -1394,7 +1463,7 @@
     /// 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) {
+    int matching(MatchingMap& matching) const {
       for (UEdgeIt it(*graph); it != INVALID; ++it) {
         matching[it] = it == anode_matching[graph->aNode(it)];
       }
@@ -1405,7 +1474,7 @@
     /// \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) {
+    bool matchingEdge(const UEdge& edge) const {
       return anode_matching[graph->aNode(edge)] == edge;
     }
 
@@ -1413,7 +1482,7 @@
     /// 
     /// Returns the matching edge from the node. If there is not such
     /// edge it gives back \c INVALID.
-    UEdge matchingEdge(const Node& node) {
+    UEdge matchingEdge(const Node& node) const {
       if (graph->aNode(node)) {
         return anode_matching[node];
       } else {
@@ -1478,6 +1547,31 @@
   
   };
 
+  /// \ingroup matching
+  ///
+  /// \brief Minimum cost maximum cardinality bipartite matching
+  ///
+  /// This function calculates the minimum cost matching of the maximum 
+  /// cardinality matchings of a bipartite graph. It gives back the matching 
+  /// in an undirected edge map.
+  ///
+  /// \param graph The bipartite graph.
+  /// \param cost The undirected edge map which contains the costs.
+  /// \retval matching The undirected edge map which will be set to 
+  /// the matching.
+  /// \return The cost of the matching.
+  template <typename BpUGraph, typename CostMap, typename MatchingMap>
+  typename CostMap::Value 
+  minCostMaxBipartiteMatching(const BpUGraph& graph, 
+                              const CostMap& cost,
+                              MatchingMap& matching) {
+    MinCostMaxBipartiteMatching<BpUGraph, CostMap> 
+      bpmatching(graph, cost);
+    bpmatching.run();
+    bpmatching.matching(matching);
+    return bpmatching.matchingCost();
+  }
+
 }
 
 #endif
diff -r 6f84cdb6f4e8 -r 0b1fc1566fdb test/bipartite_matching_test.cc
--- a/test/bipartite_matching_test.cc	Fri Apr 14 23:55:36 2006 +0000
+++ b/test/bipartite_matching_test.cc	Tue Apr 18 07:01:55 2006 +0000
@@ -23,8 +23,8 @@
   Graph graph;
   vector<Node> aNodes;
   vector<Node> bNodes;
-  int n = argc > 1 ? atoi(argv[1]) : 100;
-  int m = argc > 2 ? atoi(argv[2]) : 100;
+  int n = argc > 1 ? atoi(argv[1]) : 10;
+  int m = argc > 2 ? atoi(argv[2]) : 10;
   int e = argc > 3 ? atoi(argv[3]) : (int)((n+m) * log(n+m));
   int c = argc > 4 ? atoi(argv[4]) : 100;
 
@@ -33,6 +33,7 @@
   int max_cardinality;
   int max_weight;
   int max_cardinality_max_weight;
+  int min_cost_matching;
 
   for (int i = 0; i < n; ++i) {
     Node node = graph.addANode();
@@ -74,6 +75,21 @@
   }
 
   {
+    Graph::UEdgeMap<bool> mm(graph);
+
+    check(max_cardinality == maxBipartiteMatching(graph, mm),
+          "WRONG MATCHING");
+    
+    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");
+    }
+  }
+
+  {
     MaxBipartiteMatching<Graph> bpmatch(graph);
 
     bpmatch.greedyInit();
@@ -137,10 +153,10 @@
       
       for (UEdgeIt it(graph); it != INVALID; ++it) {
         if (mm[it]) {
-          check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+          check(pm[graph.aNode(it)] + pm[graph.bNode(it)] - weight[it] == 0,
                 "INVALID DUAL");
         } else {
-          check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+          check(pm[graph.aNode(it)] + pm[graph.bNode(it)] - weight[it] >= 0,
                 "INVALID DUAL");
         }
       }
@@ -159,6 +175,21 @@
   }
 
   {
+    Graph::UEdgeMap<bool> mm(graph);
+
+    check(max_weight == maxWeightedBipartiteMatching(graph, weight, mm),
+          "WRONG MATCHING");
+    
+    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");
+    }
+  }
+
+  {
     MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
 
     bpmatch.run();
@@ -171,10 +202,10 @@
     
     for (UEdgeIt it(graph); it != INVALID; ++it) {
       if (mm[it]) {
-        check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+        check(pm[graph.aNode(it)] + pm[graph.bNode(it)] - weight[it] == 0,
               "INVALID DUAL");
       } else {
-        check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+        check(pm[graph.aNode(it)] + pm[graph.bNode(it)] - weight[it] >= 0,
               "INVALID DUAL");
       }
     }
@@ -202,10 +233,10 @@
     
     for (UEdgeIt it(graph); it != INVALID; ++it) {
       if (mm[it]) {
-        check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+        check(pm[graph.aNode(it)] + pm[graph.bNode(it)] - weight[it] == 0,
               "INVALID DUAL");
       } else {
-        check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+        check(pm[graph.aNode(it)] + pm[graph.bNode(it)] - weight[it] >= 0,
               "INVALID DUAL");
       }
     }
@@ -223,9 +254,25 @@
   }
 
   {
-    Graph::UEdgeMap<int> cost(graph);
+    Graph::UEdgeMap<bool> mm(graph);
 
-    cost = subMap(constMap<UEdge>(c), weight);
+    check(max_cardinality_max_weight == 
+          maxWeightedMaxBipartiteMatching(graph, weight, mm),
+          "WRONG MATCHING");
+    
+    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");
+    }
+  }
+
+  Graph::UEdgeMap<int> cost(graph);
+  cost = subMap(constMap<UEdge>(c), weight);
+
+  {
 
     MinCostMaxBipartiteMatching<Graph> bpmatch(graph, cost);
 
@@ -255,11 +302,28 @@
       check(num <= 1, "INVALID PRIMAL");
     }
 
+    min_cost_matching = bpmatch.matchingCost();
     check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE");
     check(max_cardinality * c - max_cardinality_max_weight 
           == bpmatch.matchingCost(), "WRONG SIZE");
 
   }
 
+  {
+    Graph::UEdgeMap<bool> mm(graph);
+
+    check(min_cost_matching == 
+          minCostMaxBipartiteMatching(graph, cost, mm),
+          "WRONG MATCHING");
+    
+    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");
+    }
+  }
+
   return 0;
 }