# HG changeset patch
# User Peter Kovacs <kpeter@inf.elte.hu>
# Date 1242122800 -7200
# Node ID e239cde7ae572839447f31364ff9e67ae86a059a
# Parent  ca92c2f936b03d0966865203a090abbd2e703cbc
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.

diff --git a/doc/Makefile.am b/doc/Makefile.am
--- a/doc/Makefile.am
+++ b/doc/Makefile.am
@@ -8,6 +8,7 @@
 	doc/license.dox \
 	doc/mainpage.dox \
 	doc/migration.dox \
+	doc/min_cost_flow.dox \
 	doc/named-param.dox \
 	doc/namespaces.dox \
 	doc/html \
diff --git a/doc/groups.dox b/doc/groups.dox
--- a/doc/groups.dox
+++ b/doc/groups.dox
@@ -345,91 +345,16 @@
 */
 
 /**
-@defgroup min_cost_flow Minimum Cost Flow Algorithms
+@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
 @ingroup algs
 
 \brief Algorithms for finding minimum cost flows and circulations.
 
 This group contains the algorithms for finding minimum cost flows and
-circulations.
+circulations. For more information about this problem and its dual
+solution see \ref min_cost_flow "Minimum Cost Flow Problem".
 
-The \e minimum \e cost \e flow \e problem is to find a feasible flow of
-minimum total cost from a set of supply nodes to a set of demand nodes
-in a network with capacity constraints (lower and upper bounds)
-and arc costs.
-Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
-\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
-upper bounds for the flow values on the arcs, for which
-\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
-\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
-on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
-signed supply values of the nodes.
-If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
-supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
-\f$-sup(u)\f$ demand.
-A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
-of the following optimization problem.
-
-\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]
-
-The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
-zero or negative in order to have a feasible solution (since the sum
-of the expressions on the left-hand side of the inequalities is zero).
-It means that the total demand must be greater or equal to the total
-supply and all the supplies have to be carried out from the supply nodes,
-but there could be demands that are not satisfied.
-If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
-constraints have to be satisfied with equality, i.e. all demands
-have to be satisfied and all supplies have to be used.
-
-If you need the opposite inequalities in the supply/demand constraints
-(i.e. the total demand is less than the total supply and all the demands
-have to be satisfied while there could be supplies that are not used),
-then you could easily transform the problem to the above form by reversing
-the direction of the arcs and taking the negative of the supply values
-(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
-However \ref NetworkSimplex algorithm also supports this form directly
-for the sake of convenience.
-
-A feasible solution for this problem can be found using \ref Circulation.
-
-Note that the above formulation is actually more general than the usual
-definition of the minimum cost flow problem, in which strict equalities
-are required in the supply/demand contraints, i.e.
-
-\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
-    sup(u) \quad \forall u\in V. \f]
-
-However if the sum of the supply values is zero, then these two problems
-are equivalent. So if you need the equality form, you have to ensure this
-additional contraint for the algorithms.
-
-The dual solution of the minimum cost flow problem is represented by node 
-potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
-An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
-is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
-node potentials the following \e complementary \e slackness optimality
-conditions hold.
-
- - For all \f$uv\in A\f$ arcs:
-   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
-   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
-   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
- - For all \f$u\in V\f$ nodes:
-   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
-     then \f$\pi(u)=0\f$.
- 
-Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
-\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
-\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
-
-All algorithms provide dual solution (node potentials) as well,
-if an optimal flow is found.
-
-LEMON contains several algorithms for solving minimum cost flow problems.
+LEMON contains several algorithms for this problem.
  - \ref NetworkSimplex Primal Network Simplex algorithm with various
    pivot strategies.
  - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
@@ -439,10 +364,6 @@
  - \ref CancelAndTighten The Cancel and Tighten algorithm.
  - \ref CycleCanceling Cycle-Canceling algorithms.
 
-Most of these implementations support the general inequality form of the
-minimum cost flow problem, but CancelAndTighten and CycleCanceling
-only support the equality form due to the primal method they use.
-
 In general NetworkSimplex is the most efficient implementation,
 but in special cases other algorithms could be faster.
 For example, if the total supply and/or capacities are rather small,
diff --git a/doc/min_cost_flow.dox b/doc/min_cost_flow.dox
new file mode 100644
--- /dev/null
+++ b/doc/min_cost_flow.dox
@@ -0,0 +1,153 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2009
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+namespace lemon {
+
+/**
+\page min_cost_flow Minimum Cost Flow Problem
+
+\section mcf_def Definition
+
+The \e minimum \e cost \e flow \e problem is to find a feasible flow of
+minimum total cost from a set of supply nodes to a set of demand nodes
+in a network with capacity constraints (lower and upper bounds)
+and arc costs.
+
+Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$,
+\f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and
+upper bounds for the flow values on the arcs, for which
+\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
+\f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow
+on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the
+signed supply values of the nodes.
+If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
+supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
+\f$-sup(u)\f$ demand.
+A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution
+of the following optimization problem.
+
+\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]
+
+The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
+zero or negative in order to have a feasible solution (since the sum
+of the expressions on the left-hand side of the inequalities is zero).
+It means that the total demand must be greater or equal to the total
+supply and all the supplies have to be carried out from the supply nodes,
+but there could be demands that are not satisfied.
+If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
+constraints have to be satisfied with equality, i.e. all demands
+have to be satisfied and all supplies have to be used.
+
+
+\section mcf_algs Algorithms
+
+LEMON contains several algorithms for solving this problem, for more
+information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms".
+
+A feasible solution for this problem can be found using \ref Circulation.
+
+
+\section mcf_dual Dual Solution
+
+The dual solution of the minimum cost flow problem is represented by
+node potentials \f$\pi: V\rightarrow\mathbf{R}\f$.
+An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal
+if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials
+the following \e complementary \e slackness optimality conditions hold.
+
+ - For all \f$uv\in A\f$ arcs:
+   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
+   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
+   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
+ - For all \f$u\in V\f$ nodes:
+   - \f$\pi(u)<=0\f$;
+   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
+     then \f$\pi(u)=0\f$.
+ 
+Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
+\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
+\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
+
+All algorithms provide dual solution (node potentials), as well,
+if an optimal flow is found.
+
+
+\section mcf_eq Equality Form
+
+The above \ref mcf_def "definition" is actually more general than the
+usual formulation of the minimum cost flow problem, in which strict
+equalities are required in the supply/demand contraints.
+
+\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) =
+    sup(u) \quad \forall u\in V \f]
+\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
+
+However if the sum of the supply values is zero, then these two problems
+are equivalent.
+The \ref min_cost_flow_algs "algorithms" in LEMON support the general
+form, so if you need the equality form, you have to ensure this additional
+contraint manually.
+
+
+\section mcf_leq Opposite Inequality Form
+
+Another possible definition of the minimum cost flow problem is
+when there are <em>"less or equal"</em> (LEQ) supply/demand constraints,
+instead of the <em>"greater or equal"</em> (GEQ) constraints.
+
+\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.
+The equality form is also a special case of this form, of course.
+
+You could easily transform this case to the \ref mcf_def "GEQ form"
+of the problem by reversing the direction of the arcs and taking the
+negative of the supply values (e.g. using \ref ReverseDigraph and
+\ref NegMap adaptors).
+However \ref NetworkSimplex algorithm also supports this form directly
+for the sake of convenience.
+
+Note that the optimality conditions for this supply constraint type are
+slightly differ from the conditions that are discussed for the GEQ form,
+namely the potentials have to be non-negative instead of non-positive.
+An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem
+is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$
+node potentials the following conditions hold.
+
+ - For all \f$uv\in A\f$ arcs:
+   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
+   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
+   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
+ - For all \f$u\in V\f$ nodes:
+   - \f$\pi(u)>=0\f$;
+   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
+     then \f$\pi(u)=0\f$.
+
+*/
+}
diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
--- a/lemon/network_simplex.h
+++ b/lemon/network_simplex.h
@@ -19,7 +19,7 @@
 #ifndef LEMON_NETWORK_SIMPLEX_H
 #define LEMON_NETWORK_SIMPLEX_H
 
-/// \ingroup min_cost_flow
+/// \ingroup min_cost_flow_algs
 ///
 /// \file
 /// \brief Network Simplex algorithm for finding a minimum cost flow.
@@ -33,7 +33,7 @@
 
 namespace lemon {
 
-  /// \addtogroup min_cost_flow
+  /// \addtogroup min_cost_flow_algs
   /// @{
 
   /// \brief Implementation of the primal Network Simplex algorithm
@@ -102,50 +102,16 @@
     /// i.e. the direction of the inequalities in the supply/demand
     /// 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
     };
     
     /// \brief Constants for selecting the pivot rule.
@@ -215,6 +181,8 @@
     const GR &_graph;
     int _node_num;
     int _arc_num;
+    int _all_arc_num;
+    int _search_arc_num;
 
     // Parameters of the problem
     bool _have_lower;
@@ -277,7 +245,7 @@
       const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
       // Pivot rule data
       int _next_arc;
@@ -288,13 +256,14 @@
       FirstEligiblePivotRule(NetworkSimplex &ns) :
         _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)
       {}
 
       // Find next entering arc
       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) {
             _in_arc = e;
@@ -328,7 +297,7 @@
       const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
     public:
 
@@ -336,13 +305,13 @@
       BestEligiblePivotRule(NetworkSimplex &ns) :
         _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)
       {}
 
       // Find next entering arc
       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) {
             min = c;
@@ -367,7 +336,7 @@
       const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
       // Pivot rule data
       int _block_size;
@@ -379,14 +348,15 @@
       BlockSearchPivotRule(NetworkSimplex &ns) :
         _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 );
       }
 
@@ -395,7 +365,7 @@
         Cost c, min = 0;
         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) {
             min = c;
@@ -440,7 +410,7 @@
       const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
       // Pivot rule data
       IntVector _candidates;
@@ -454,7 +424,8 @@
       CandidateListPivotRule(NetworkSimplex &ns) :
         _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 LIST_LENGTH_FACTOR = 1.0;
@@ -463,7 +434,7 @@
         const int MIN_MINOR_LIMIT = 3;
 
         _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),
                                  MIN_MINOR_LIMIT );
@@ -500,7 +471,7 @@
         // Major iteration: build a new candidate list
         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) {
             _candidates[_curr_length++] = e;
@@ -546,7 +517,7 @@
       const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
-      int _arc_num;
+      int _search_arc_num;
 
       // Pivot rule data
       int _block_size, _head_length, _curr_length;
@@ -574,8 +545,8 @@
       AlteringListPivotRule(NetworkSimplex &ns) :
         _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
         const double BLOCK_SIZE_FACTOR = 1.5;
@@ -584,7 +555,7 @@
         const int MIN_HEAD_LENGTH = 3;
 
         _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),
                                  MIN_HEAD_LENGTH );
@@ -610,7 +581,7 @@
         int last_arc = 0;
         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]]);
           if (_cand_cost[e] < 0) {
@@ -678,17 +649,17 @@
       _node_num = countNodes(_graph);
       _arc_num = countArcs(_graph);
       int all_node_num = _node_num + 1;
-      int all_arc_num = _arc_num + _node_num;
+      int max_arc_num = _arc_num + 2 * _node_num;
 
-      _source.resize(all_arc_num);
-      _target.resize(all_arc_num);
+      _source.resize(max_arc_num);
+      _target.resize(max_arc_num);
 
-      _lower.resize(all_arc_num);
-      _upper.resize(all_arc_num);
-      _cap.resize(all_arc_num);
-      _cost.resize(all_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);
 
       _parent.resize(all_node_num);
@@ -698,7 +669,7 @@
       _rev_thread.resize(all_node_num);
       _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)
       int i = 0;
@@ -1069,7 +1040,7 @@
       // Initialize artifical cost
       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();
         for (int i = 0; i != _arc_num; ++i) {
@@ -1093,29 +1064,121 @@
       _succ_num[_root] = _node_num + 1;
       _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;
+      }
 
       return true;
     }
@@ -1374,20 +1437,8 @@
       }
       
       // 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;
       }
 
       // Transform the solution and the supply map to the original form
@@ -1401,6 +1452,30 @@
           }
         }
       }
+      
+      // 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;
     }