Ticket #180: 180-cas4-08f7479bbc45.patch

lemon/capacity_scaling.h

@@ -55 +55 @@
   ///
   /// \warning Both value types must be signed and all input data must
   /// be integer.
-  /// \warning This implementation can handle only non-negative arc costs.
+  /// \warning This algorithm does not support negative costs for such
+  /// arcs that have infinite upper bound.
   template <typename GR, typename V = int, typename C = V>
   class CapacityScaling
   {
@@ -79 +80 @@
       /// bounded), and the algorithm has found optimal flow and node
       /// potentials (primal and dual solutions).
       OPTIMAL,
-      /// The objective function of the problem is unbounded, i.e.
-      /// there is a directed cycle having negative total cost and
-      /// infinite upper bound.
+      /// The digraph contains an arc of negative cost and infinite
+      /// upper bound. It means that the objective function is unbounded
+      /// on that arc, however note that it could actually be bounded
+      /// over the feasible flows, but this algroithm cannot handle
+      /// these cases.
       UNBOUNDED
     };
   
@@ -459 +462 @@
     /// \n \c OPTIMAL if the problem has optimal solution
     /// (i.e. it is feasible and bounded), and the algorithm has found
     /// optimal flow and node potentials (primal and dual solutions),
-    /// \n \c UNBOUNDED if the objective function of the problem is
-    /// unbounded, i.e. there is a directed cycle having negative total
-    /// cost and infinite upper bound.
+    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
+    /// and infinite upper bound. It means that the objective function
+    /// is unbounded on that arc, however note that it could actually be
+    /// bounded over the feasible flows, but this algroithm cannot handle
+    /// these cases.
     ///
     /// \see ProblemType
     ProblemType run(bool scaling = true) {
-      if (!init(scaling)) return INFEASIBLE;
+      ProblemType pt = init(scaling);
+      if (pt != OPTIMAL) return pt;
       return start();
     }
 
@@ -609 +615 @@
   private:
 
     // Initialize the algorithm
-    bool init(bool scaling) {
-      if (_node_num == 0) return false;
+    ProblemType init(bool scaling) {
+      if (_node_num == 0) return INFEASIBLE;
 
       // Check the sum of supply values
       _sum_supply = 0;
       for (int i = 0; i != _root; ++i) {
         _sum_supply += _supply[i];
       }
-      if (_sum_supply > 0) return false;
+      if (_sum_supply > 0) return INFEASIBLE;
       
       // Initialize maps
       for (int i = 0; i != _root; ++i) {
@@ -648 +654 @@
           _res_cap[j] = _forward[j] ? _upper[j] : 0;
         }
       }
+
+      // Handle negative costs
+      for (int u = 0; u != _root; ++u) {
+        for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
+          Value rc = _res_cap[a];
+          if (_cost[a] < 0 && rc > 0) {
+            if (rc == INF) return UNBOUNDED;
+            _excess[u] -= rc;
+            _excess[_target[a]] += rc;
+            _res_cap[a] = 0;
+            _res_cap[_reverse[a]] += rc;
+          }
+        }
+      }
       
       // Handle GEQ supply type
       if (_sum_supply < 0) {
@@ -696 +716 @@
         _delta = 1;
       }
 
-      return true;
+      return OPTIMAL;
     }
 
     ProblemType start() {

test/min_cost_flow_test.cc

@@ -501 +501 @@
     mcf.supplyMap(s6);
     checkMcf(mcf, mcf.run(),
              gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#C15", GEQ);
+
+    // Check negative costs
+    CapacityScaling<Digraph> neg_mcf(neg_gr);
+    neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
+    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
+      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#C16");
+    neg_mcf.upperMap(neg_u2);
+    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
+      neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#C17");
+    neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
+    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,
+      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#C18");
+
+    CapacityScaling<Digraph> negs_mcf(negs_gr);
+    negs_mcf.costMap(negs_c).supplyMap(negs_s);
+    checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
+      negs_c, negs_s, negs_mcf.UNBOUNDED, false, 0, "#C19", GEQ);
+    negs_mcf.upperMap(negs_u);
+    checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
+      negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#C20", GEQ);
   }
 
   return 0;