Changeset 663:8b0df68370a4 in lemon-1.2 for lemon

Timestamp:

05/12/09 12:06:40 (14 years ago)

Author:

Peter Kovacs <kpeter@…>

Branch:

default

Phase:

public

Message:

Fix the GEQ/LEQ handling in NetworkSimplex? + improve doc (#291)

• Fix the optimality conditions for the GEQ/LEQ form.
• Fix the initialization of the algortihm. It ensures correct solutions and it is much faster for the inequality forms.
• Fix the pivot rules to search all the arcs that have to be allowed to get in the basis.
• Better block size for the Block Search pivot rule.
• Improve documentation of the problem and move it to a separate page.

File:

• lemon/network_simplex.h (modified) (27 diffs)

lemon/network_simplex.h

@@ -20 +20 @@
 #define LEMON_NETWORK_SIMPLEX_H
 
-/// \ingroup min_cost_flow
+/// \ingroup min_cost_flow_algs
 ///
 /// \file
@@ -34 +34 @@
 namespace lemon {
 
-  /// \addtogroup min_cost_flow
+  /// \addtogroup min_cost_flow_algs
   /// @{
 
@@ -103 +103 @@
     /// constraints of the \ref min_cost_flow "minimum cost flow problem".
     ///
-    /// The default supply type is \c GEQ, since this form is supported
-    /// by other minimum cost flow algorithms and the \ref Circulation
-    /// algorithm, as well.
-    /// The \c LEQ problem type can be selected using the \ref supplyType()
-    /// function.
-    ///
-    /// Note that the equality form is a special case of both supply types.
+    /// The default supply type is \c GEQ, the \c LEQ type can be
+    /// selected using \ref supplyType().
+    /// The equality form is a special case of both supply types.
     enum SupplyType {
-
       /// This option means that there are <em>"greater or equal"</em>
-      /// supply/demand constraints in the definition, i.e. the exact
-      /// formulation of the problem is the following.
-      /**
-          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
-          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
-              sup(u) \quad \forall u\in V \f]
-          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
-      */
-      /// It means that the total demand must be greater or equal to the 
-      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
-      /// negative) and all the supplies have to be carried out from 
-      /// the supply nodes, but there could be demands that are not 
-      /// satisfied.
+      /// supply/demand constraints in the definition of the problem.
       GEQ,
-      /// It is just an alias for the \c GEQ option.
-      CARRY_SUPPLIES = GEQ,
-
       /// This option means that there are <em>"less or equal"</em>
-      /// supply/demand constraints in the definition, i.e. the exact
-      /// formulation of the problem is the following.
-      /**
-          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
-          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
-              sup(u) \quad \forall u\in V \f]
-          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
-      */
-      /// It means that the total demand must be less or equal to the 
-      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
-      /// positive) and all the demands have to be satisfied, but there
-      /// could be supplies that are not carried out from the supply
-      /// nodes.
-      LEQ,
-      /// It is just an alias for the \c LEQ option.
-      SATISFY_DEMANDS = LEQ
+      /// supply/demand constraints in the definition of the problem.
+      LEQ
     };
     
@@ -216 +182 @@
     int _node_num;
     int _arc_num;
+    int _all_arc_num;
+    int _search_arc_num;
 
     // Parameters of the problem
@@ -278 +246 @@
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
       // Pivot rule data
@@ -289 +257 @@
         _source(ns._source), _target(ns._target),
         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
-        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
+        _next_arc(0)
       {}
 
@@ -295 +264 @@
       bool findEnteringArc() {
         Cost c;
-        for (int e = _next_arc; e < _arc_num; ++e) {
+        for (int e = _next_arc; e < _search_arc_num; ++e) {
           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
           if (c < 0) {
@@ -329 +298 @@
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
     public:
@@ -337 +306 @@
         _source(ns._source), _target(ns._target),
         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
-        _in_arc(ns.in_arc), _arc_num(ns._arc_num)
+        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
       {}
 
@@ -343 +312 @@
       bool findEnteringArc() {
         Cost c, min = 0;
-        for (int e = 0; e < _arc_num; ++e) {
+        for (int e = 0; e < _search_arc_num; ++e) {
           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
           if (c < min) {
@@ -368 +337 @@
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
       // Pivot rule data
@@ -380 +349 @@
         _source(ns._source), _target(ns._target),
         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
-        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
+        _next_arc(0)
       {
         // The main parameters of the pivot rule
-        const double BLOCK_SIZE_FACTOR = 2.0;
+        const double BLOCK_SIZE_FACTOR = 0.5;
         const int MIN_BLOCK_SIZE = 10;
 
         _block_size = std::max( int(BLOCK_SIZE_FACTOR *
-                                    std::sqrt(double(_arc_num))),
+                                    std::sqrt(double(_search_arc_num))),
                                 MIN_BLOCK_SIZE );
       }
@@ -396 +366 @@
         int cnt = _block_size;
         int e, min_arc = _next_arc;
-        for (e = _next_arc; e < _arc_num; ++e) {
+        for (e = _next_arc; e < _search_arc_num; ++e) {
           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
           if (c < min) {
@@ -441 +411 @@
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
       // Pivot rule data
@@ -455 +425 @@
         _source(ns._source), _target(ns._target),
         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
-        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
+        _next_arc(0)
       {
         // The main parameters of the pivot rule
@@ -464 +435 @@
 
         _list_length = std::max( int(LIST_LENGTH_FACTOR *
-                                     std::sqrt(double(_arc_num))),
+                                     std::sqrt(double(_search_arc_num))),
                                  MIN_LIST_LENGTH );
         _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
@@ -501 +472 @@
         min = 0;
         _curr_length = 0;
-        for (e = _next_arc; e < _arc_num; ++e) {
+        for (e = _next_arc; e < _search_arc_num; ++e) {
           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
           if (c < 0) {
@@ -547 +518 @@
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
       // Pivot rule data
@@ -575 +546 @@
         _source(ns._source), _target(ns._target),
         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
-        _in_arc(ns.in_arc), _arc_num(ns._arc_num),
-        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
+        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
+        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
       {
         // The main parameters of the pivot rule
@@ -585 +556 @@
 
         _block_size = std::max( int(BLOCK_SIZE_FACTOR *
-                                    std::sqrt(double(_arc_num))),
+                                    std::sqrt(double(_search_arc_num))),
                                 MIN_BLOCK_SIZE );
         _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
@@ -611 +582 @@
         int limit = _head_length;
 
-        for (int e = _next_arc; e < _arc_num; ++e) {
+        for (int e = _next_arc; e < _search_arc_num; ++e) {
           _cand_cost[e] = _state[e] *
             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
@@ -679 +650 @@
       _arc_num = countArcs(_graph);
       int all_node_num = _node_num + 1;
-      int all_arc_num = _arc_num + _node_num;
-
-      _source.resize(all_arc_num);
-      _target.resize(all_arc_num);
-
-      _lower.resize(all_arc_num);
-      _upper.resize(all_arc_num);
-      _cap.resize(all_arc_num);
-      _cost.resize(all_arc_num);
+      int max_arc_num = _arc_num + 2 * _node_num;
+
+      _source.resize(max_arc_num);
+      _target.resize(max_arc_num);
+
+      _lower.resize(_arc_num);
+      _upper.resize(_arc_num);
+      _cap.resize(max_arc_num);
+      _cost.resize(max_arc_num);
       _supply.resize(all_node_num);
-      _flow.resize(all_arc_num);
+      _flow.resize(max_arc_num);
       _pi.resize(all_node_num);
 
@@ -699 +670 @@
       _succ_num.resize(all_node_num);
       _last_succ.resize(all_node_num);
-      _state.resize(all_arc_num);
+      _state.resize(max_arc_num);
 
       // Copy the graph (store the arcs in a mixed order)
@@ -1070 +1041 @@
       Cost ART_COST;
       if (std::numeric_limits<Cost>::is_exact) {
-        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
+        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
       } else {
         ART_COST = std::numeric_limits<Cost>::min();
@@ -1094 +1065 @@
       _last_succ[_root] = _root - 1;
       _supply[_root] = -_sum_supply;
-      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
+      _pi[_root] = 0;
 
       // Add artificial arcs and initialize the spanning tree data structure
-      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
-        _parent[u] = _root;
-        _pred[u] = e;
-        _thread[u] = u + 1;
-        _rev_thread[u + 1] = u;
-        _succ_num[u] = 1;
-        _last_succ[u] = u;
-        _cost[e] = ART_COST;
-        _cap[e] = INF;
-        _state[e] = STATE_TREE;
-        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
-          _flow[e] = _supply[u];
-          _forward[u] = true;
-          _pi[u] = -ART_COST + _pi[_root];
-        } else {
-          _flow[e] = -_supply[u];
-          _forward[u] = false;
-          _pi[u] = ART_COST + _pi[_root];
-        }
+      if (_sum_supply == 0) {
+        // EQ supply constraints
+        _search_arc_num = _arc_num;
+        _all_arc_num = _arc_num + _node_num;
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          _parent[u] = _root;
+          _pred[u] = e;
+          _thread[u] = u + 1;
+          _rev_thread[u + 1] = u;
+          _succ_num[u] = 1;
+          _last_succ[u] = u;
+          _cap[e] = INF;
+          _state[e] = STATE_TREE;
+          if (_supply[u] >= 0) {
+            _forward[u] = true;
+            _pi[u] = 0;
+            _source[e] = u;
+            _target[e] = _root;
+            _flow[e] = _supply[u];
+            _cost[e] = 0;
+          } else {
+            _forward[u] = false;
+            _pi[u] = ART_COST;
+            _source[e] = _root;
+            _target[e] = u;
+            _flow[e] = -_supply[u];
+            _cost[e] = ART_COST;
+          }
+        }
+      }
+      else if (_sum_supply > 0) {
+        // LEQ supply constraints
+        _search_arc_num = _arc_num + _node_num;
+        int f = _arc_num + _node_num;
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          _parent[u] = _root;
+          _thread[u] = u + 1;
+          _rev_thread[u + 1] = u;
+          _succ_num[u] = 1;
+          _last_succ[u] = u;
+          if (_supply[u] >= 0) {
+            _forward[u] = true;
+            _pi[u] = 0;
+            _pred[u] = e;
+            _source[e] = u;
+            _target[e] = _root;
+            _cap[e] = INF;
+            _flow[e] = _supply[u];
+            _cost[e] = 0;
+            _state[e] = STATE_TREE;
+          } else {
+            _forward[u] = false;
+            _pi[u] = ART_COST;
+            _pred[u] = f;
+            _source[f] = _root;
+            _target[f] = u;
+            _cap[f] = INF;
+            _flow[f] = -_supply[u];
+            _cost[f] = ART_COST;
+            _state[f] = STATE_TREE;
+            _source[e] = u;
+            _target[e] = _root;
+            _cap[e] = INF;
+            _flow[e] = 0;
+            _cost[e] = 0;
+            _state[e] = STATE_LOWER;
+            ++f;
+          }
+        }
+        _all_arc_num = f;
+      }
+      else {
+        // GEQ supply constraints
+        _search_arc_num = _arc_num + _node_num;
+        int f = _arc_num + _node_num;
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          _parent[u] = _root;
+          _thread[u] = u + 1;
+          _rev_thread[u + 1] = u;
+          _succ_num[u] = 1;
+          _last_succ[u] = u;
+          if (_supply[u] <= 0) {
+            _forward[u] = false;
+            _pi[u] = 0;
+            _pred[u] = e;
+            _source[e] = _root;
+            _target[e] = u;
+            _cap[e] = INF;
+            _flow[e] = -_supply[u];
+            _cost[e] = 0;
+            _state[e] = STATE_TREE;
+          } else {
+            _forward[u] = true;
+            _pi[u] = -ART_COST;
+            _pred[u] = f;
+            _source[f] = u;
+            _target[f] = _root;
+            _cap[f] = INF;
+            _flow[f] = _supply[u];
+            _state[f] = STATE_TREE;
+            _cost[f] = ART_COST;
+            _source[e] = _root;
+            _target[e] = u;
+            _cap[e] = INF;
+            _flow[e] = 0;
+            _cost[e] = 0;
+            _state[e] = STATE_LOWER;
+            ++f;
+          }
+        }
+        _all_arc_num = f;
       }
 
@@ -1375 +1438 @@
       
       // Check feasibility
-      if (_sum_supply < 0) {
-        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
-          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
-        }
-      }
-      else if (_sum_supply > 0) {
-        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
-          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
-        }
-      }
-      else {
-        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
-          if (_flow[e] != 0) return INFEASIBLE;
-        }
+      for (int e = _search_arc_num; e != _all_arc_num; ++e) {
+        if (_flow[e] != 0) return INFEASIBLE;
       }
 
@@ -1402 +1453 @@
         }
       }
+      
+      // Shift potentials to meet the requirements of the GEQ/LEQ type
+      // optimality conditions
+      if (_sum_supply == 0) {
+        if (_stype == GEQ) {
+          Cost max_pot = std::numeric_limits<Cost>::min();
+          for (int i = 0; i != _node_num; ++i) {
+            if (_pi[i] > max_pot) max_pot = _pi[i];
+          }
+          if (max_pot > 0) {
+            for (int i = 0; i != _node_num; ++i)
+              _pi[i] -= max_pot;
+          }
+        } else {
+          Cost min_pot = std::numeric_limits<Cost>::max();
+          for (int i = 0; i != _node_num; ++i) {
+            if (_pi[i] < min_pot) min_pot = _pi[i];
+          }
+          if (min_pot < 0) {
+            for (int i = 0; i != _node_num; ++i)
+              _pi[i] -= min_pot;
+          }
+        }
+      }
 
       return OPTIMAL;