# HG changeset patch
# User kpeter
# Date 1233957154 0
# Node ID e98bbe64cca4f68e74f2da559d49d45bceee69c4
# Parent  4f47c0f6be046a6005e7a4c94e8b34ac3e99971e
Rework Network Simplex
Use simpler and faster graph implementation instead of SmartGraph

diff -r 4f47c0f6be04 -r e98bbe64cca4 lemon/network_simplex.h
--- a/lemon/network_simplex.h	Wed Feb 04 14:42:31 2009 +0000
+++ b/lemon/network_simplex.h	Fri Feb 06 21:52:34 2009 +0000
@@ -28,9 +28,7 @@
 #include <limits>
 #include <algorithm>
 
-#include <lemon/graph_adaptor.h>
 #include <lemon/graph_utils.h>
-#include <lemon/smart_graph.h>
 #include <lemon/math.h>
 
 namespace lemon {
@@ -72,29 +70,21 @@
              typename SupplyMap = typename Graph::template NodeMap<int> >
   class NetworkSimplex
   {
+    GRAPH_TYPEDEFS(typename Graph);
+
     typedef typename CapacityMap::Value Capacity;
     typedef typename CostMap::Value Cost;
     typedef typename SupplyMap::Value Supply;
 
-    typedef SmartGraph SGraph;
-    GRAPH_TYPEDEFS(typename SGraph);
+    typedef std::vector<Edge> EdgeVector;
+    typedef std::vector<Node> NodeVector;
+    typedef std::vector<int> IntVector;
+    typedef std::vector<bool> BoolVector;
+    typedef std::vector<Capacity> CapacityVector;
+    typedef std::vector<Cost> CostVector;
+    typedef std::vector<Supply> SupplyVector;
 
-    typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
-    typedef typename SGraph::template EdgeMap<Cost> SCostMap;
-    typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
-    typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
-
-    typedef typename SGraph::template NodeMap<int> IntNodeMap;
-    typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
-    typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
-    typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
-    typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
-    typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;
-
-    typedef typename Graph::template NodeMap<Node> NodeRefMap;
-    typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
-
-    typedef std::vector<Edge> EdgeVector;
+    typedef typename Graph::template NodeMap<int> IntNodeMap;
 
   public:
 
@@ -116,35 +106,6 @@
 
   private:
 
-    /// \brief Map adaptor class for handling reduced edge costs.
-    ///
-    /// Map adaptor class for handling reduced edge costs.
-    class ReducedCostMap : public MapBase<Edge, Cost>
-    {
-    private:
-
-      const SGraph &_gr;
-      const SCostMap &_cost_map;
-      const SPotentialMap &_pot_map;
-
-    public:
-
-      ///\e
-      ReducedCostMap( const SGraph &gr,
-                      const SCostMap &cost_map,
-                      const SPotentialMap &pot_map ) :
-        _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
-
-      ///\e
-      Cost operator[](const Edge &e) const {
-        return _cost_map[e] + _pot_map[_gr.source(e)]
-                            - _pot_map[_gr.target(e)];
-      }
-
-    }; //class ReducedCostMap
-
-  private:
-
     /// \brief Implementation of the "First Eligible" pivot rule for the
     /// \ref NetworkSimplex "network simplex" algorithm.
     ///
@@ -157,40 +118,52 @@
     private:
 
       // References to the NetworkSimplex class
-      NetworkSimplex &_ns;
-      EdgeVector &_edges;
+      const EdgeVector &_edge;
+      const IntVector  &_source;
+      const IntVector  &_target;
+      const CostVector &_cost;
+      const IntVector  &_state;
+      const CostVector &_pi;
+      int &_in_edge;
+      int _edge_num;
 
+      // Pivot rule data
       int _next_edge;
 
     public:
 
       /// Constructor
-      FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
-        _ns(ns), _edges(edges), _next_edge(0) {}
+      FirstEligiblePivotRule(NetworkSimplex &ns) :
+        _edge(ns._edge), _source(ns._source), _target(ns._target),
+        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
+        _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
+      {}
 
       /// Find next entering edge
       bool findEnteringEdge() {
-        Edge e;
-        for (int i = _next_edge; i < int(_edges.size()); ++i) {
-          e = _edges[i];
-          if (_ns._state[e] * _ns._red_cost[e] < 0) {
-            _ns._in_edge = e;
-            _next_edge = i + 1;
+        Cost c;
+        for (int e = _next_edge; e < _edge_num; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < 0) {
+            _in_edge = e;
+            _next_edge = e + 1;
             return true;
           }
         }
-        for (int i = 0; i < _next_edge; ++i) {
-          e = _edges[i];
-          if (_ns._state[e] * _ns._red_cost[e] < 0) {
-            _ns._in_edge = e;
-            _next_edge = i + 1;
+        for (int e = 0; e < _next_edge; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < 0) {
+            _in_edge = e;
+            _next_edge = e + 1;
             return true;
           }
         }
         return false;
       }
+
     }; //class FirstEligiblePivotRule
 
+
     /// \brief Implementation of the "Best Eligible" pivot rule for the
     /// \ref NetworkSimplex "network simplex" algorithm.
     ///
@@ -203,30 +176,40 @@
     private:
 
       // References to the NetworkSimplex class
-      NetworkSimplex &_ns;
-      EdgeVector &_edges;
+      const EdgeVector &_edge;
+      const IntVector  &_source;
+      const IntVector  &_target;
+      const CostVector &_cost;
+      const IntVector  &_state;
+      const CostVector &_pi;
+      int &_in_edge;
+      int _edge_num;
 
     public:
 
       /// Constructor
-      BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
-        _ns(ns), _edges(edges) {}
+      BestEligiblePivotRule(NetworkSimplex &ns) :
+        _edge(ns._edge), _source(ns._source), _target(ns._target),
+        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
+        _in_edge(ns._in_edge), _edge_num(ns._edge_num)
+      {}
 
       /// Find next entering edge
       bool findEnteringEdge() {
-        Cost min = 0;
-        Edge e;
-        for (int i = 0; i < int(_edges.size()); ++i) {
-          e = _edges[i];
-          if (_ns._state[e] * _ns._red_cost[e] < min) {
-            min = _ns._state[e] * _ns._red_cost[e];
-            _ns._in_edge = e;
+        Cost c, min = 0;
+        for (int e = 0; e < _edge_num; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < min) {
+            min = c;
+            _in_edge = e;
           }
         }
         return min < 0;
       }
+
     }; //class BestEligiblePivotRule
 
+
     /// \brief Implementation of the "Block Search" pivot rule for the
     /// \ref NetworkSimplex "network simplex" algorithm.
     ///
@@ -239,37 +222,45 @@
     private:
 
       // References to the NetworkSimplex class
-      NetworkSimplex &_ns;
-      EdgeVector &_edges;
+      const EdgeVector &_edge;
+      const IntVector  &_source;
+      const IntVector  &_target;
+      const CostVector &_cost;
+      const IntVector  &_state;
+      const CostVector &_pi;
+      int &_in_edge;
+      int _edge_num;
 
+      // Pivot rule data
       int _block_size;
-      int _next_edge, _min_edge;
+      int _next_edge;
 
     public:
 
       /// Constructor
-      BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
-        _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
+      BlockSearchPivotRule(NetworkSimplex &ns) :
+        _edge(ns._edge), _source(ns._source), _target(ns._target),
+        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
+        _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
       {
         // The main parameters of the pivot rule
         const double BLOCK_SIZE_FACTOR = 2.0;
         const int MIN_BLOCK_SIZE = 10;
 
-        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
+        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
                                 MIN_BLOCK_SIZE );
       }
 
       /// Find next entering edge
       bool findEnteringEdge() {
-        Cost curr, min = 0;
-        Edge e;
+        Cost c, min = 0;
         int cnt = _block_size;
-        int i;
-        for (i = _next_edge; i < int(_edges.size()); ++i) {
-          e = _edges[i];
-          if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
-            min = curr;
-            _min_edge = i;
+        int e, min_edge = _next_edge;
+        for (e = _next_edge; e < _edge_num; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < min) {
+            min = c;
+            min_edge = e;
           }
           if (--cnt == 0) {
             if (min < 0) break;
@@ -277,11 +268,11 @@
           }
         }
         if (min == 0 || cnt > 0) {
-          for (i = 0; i < _next_edge; ++i) {
-            e = _edges[i];
-            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
-              min = curr;
-              _min_edge = i;
+          for (e = 0; e < _next_edge; ++e) {
+            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+            if (c < min) {
+              min = c;
+              min_edge = e;
             }
             if (--cnt == 0) {
               if (min < 0) break;
@@ -290,12 +281,14 @@
           }
         }
         if (min >= 0) return false;
-        _ns._in_edge = _edges[_min_edge];
-        _next_edge = i;
+        _in_edge = min_edge;
+        _next_edge = e;
         return true;
       }
+
     }; //class BlockSearchPivotRule
 
+
     /// \brief Implementation of the "Candidate List" pivot rule for the
     /// \ref NetworkSimplex "network simplex" algorithm.
     ///
@@ -308,19 +301,28 @@
     private:
 
       // References to the NetworkSimplex class
-      NetworkSimplex &_ns;
-      EdgeVector &_edges;
+      const EdgeVector &_edge;
+      const IntVector  &_source;
+      const IntVector  &_target;
+      const CostVector &_cost;
+      const IntVector  &_state;
+      const CostVector &_pi;
+      int &_in_edge;
+      int _edge_num;
 
-      EdgeVector _candidates;
+      // Pivot rule data
+      IntVector _candidates;
       int _list_length, _minor_limit;
       int _curr_length, _minor_count;
-      int _next_edge, _min_edge;
+      int _next_edge;
 
     public:
 
       /// Constructor
-      CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
-        _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
+      CandidateListPivotRule(NetworkSimplex &ns) :
+        _edge(ns._edge), _source(ns._source), _target(ns._target),
+        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
+        _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
       {
         // The main parameters of the pivot rule
         const double LIST_LENGTH_FACTOR = 1.0;
@@ -328,7 +330,7 @@
         const double MINOR_LIMIT_FACTOR = 0.1;
         const int MIN_MINOR_LIMIT = 3;
 
-        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())),
+        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge_num)),
                                  MIN_LIST_LENGTH );
         _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
                                  MIN_MINOR_LIMIT );
@@ -338,51 +340,52 @@
 
       /// Find next entering edge
       bool findEnteringEdge() {
-        Cost min, curr;
+        Cost min, c;
+        int e, min_edge = _next_edge;
         if (_curr_length > 0 && _minor_count < _minor_limit) {
           // Minor iteration: select the best eligible edge from the
           // current candidate list
           ++_minor_count;
-          Edge e;
           min = 0;
           for (int i = 0; i < _curr_length; ++i) {
             e = _candidates[i];
-            curr = _ns._state[e] * _ns._red_cost[e];
-            if (curr < min) {
-              min = curr;
-              _ns._in_edge = e;
+            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+            if (c < min) {
+              min = c;
+              min_edge = e;
             }
-            if (curr >= 0) {
+            if (c >= 0) {
               _candidates[i--] = _candidates[--_curr_length];
             }
           }
-          if (min < 0) return true;
+          if (min < 0) {
+            _in_edge = min_edge;
+            return true;
+          }
         }
 
         // Major iteration: build a new candidate list
-        Edge e;
         min = 0;
         _curr_length = 0;
-        int i;
-        for (i = _next_edge; i < int(_edges.size()); ++i) {
-          e = _edges[i];
-          if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+        for (e = _next_edge; e < _edge_num; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < 0) {
             _candidates[_curr_length++] = e;
-            if (curr < min) {
-              min = curr;
-              _min_edge = i;
+            if (c < min) {
+              min = c;
+              min_edge = e;
             }
             if (_curr_length == _list_length) break;
           }
         }
         if (_curr_length < _list_length) {
-          for (i = 0; i < _next_edge; ++i) {
-            e = _edges[i];
-            if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+          for (e = 0; e < _next_edge; ++e) {
+            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+            if (c < 0) {
               _candidates[_curr_length++] = e;
-              if (curr < min) {
-                min = curr;
-                _min_edge = i;
+              if (c < min) {
+                min = c;
+                min_edge = e;
               }
               if (_curr_length == _list_length) break;
             }
@@ -390,12 +393,14 @@
         }
         if (_curr_length == 0) return false;
         _minor_count = 1;
-        _ns._in_edge = _edges[_min_edge];
-        _next_edge = i;
+        _in_edge = min_edge;
+        _next_edge = e;
         return true;
       }
+
     }; //class CandidateListPivotRule
 
+
     /// \brief Implementation of the "Altering Candidate List" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
     ///
@@ -408,23 +413,29 @@
     private:
 
       // References to the NetworkSimplex class
-      NetworkSimplex &_ns;
-      EdgeVector &_edges;
+      const EdgeVector &_edge;
+      const IntVector  &_source;
+      const IntVector  &_target;
+      const CostVector &_cost;
+      const IntVector  &_state;
+      const CostVector &_pi;
+      int &_in_edge;
+      int _edge_num;
 
-      EdgeVector _candidates;
-      SCostMap _cand_cost;
       int _block_size, _head_length, _curr_length;
       int _next_edge;
+      IntVector _candidates;
+      CostVector _cand_cost;
 
       // Functor class to compare edges during sort of the candidate list
       class SortFunc
       {
       private:
-        const SCostMap &_map;
+        const CostVector &_map;
       public:
-        SortFunc(const SCostMap &map) : _map(map) {}
-        bool operator()(const Edge &e1, const Edge &e2) {
-          return _map[e1] > _map[e2];
+        SortFunc(const CostVector &map) : _map(map) {}
+        bool operator()(int left, int right) {
+          return _map[left] > _map[right];
         }
       };
 
@@ -433,9 +444,11 @@
     public:
 
       /// Constructor
-      AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
-        _ns(ns), _edges(edges), _cand_cost(_ns._graph),
-        _next_edge(0), _sort_func(_cand_cost)
+      AlteringListPivotRule(NetworkSimplex &ns) :
+        _edge(ns._edge), _source(ns._source), _target(ns._target),
+        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
+        _in_edge(ns._in_edge), _edge_num(ns._edge_num),
+         _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost)
       {
         // The main parameters of the pivot rule
         const double BLOCK_SIZE_FACTOR = 1.5;
@@ -443,7 +456,7 @@
         const double HEAD_LENGTH_FACTOR = 0.1;
         const int MIN_HEAD_LENGTH = 3;
 
-        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
+        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
                                 MIN_BLOCK_SIZE );
         _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
                                  MIN_HEAD_LENGTH );
@@ -454,11 +467,13 @@
       /// Find next entering edge
       bool findEnteringEdge() {
         // Check the current candidate list
-        Edge e;
-        for (int ix = 0; ix < _curr_length; ++ix) {
-          e = _candidates[ix];
-          if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
-            _candidates[ix--] = _candidates[--_curr_length];
+        int e;
+        for (int i = 0; i < _curr_length; ++i) {
+          e = _candidates[i];
+          _cand_cost[e] = _state[e] * 
+            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (_cand_cost[e] >= 0) {
+            _candidates[i--] = _candidates[--_curr_length];
           }
         }
 
@@ -466,11 +481,13 @@
         int cnt = _block_size;
         int last_edge = 0;
         int limit = _head_length;
-        for (int i = _next_edge; i < int(_edges.size()); ++i) {
-          e = _edges[i];
-          if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+        
+        for (int e = _next_edge; e < _edge_num; ++e) {
+          _cand_cost[e] = _state[e] *
+            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (_cand_cost[e] < 0) {
             _candidates[_curr_length++] = e;
-            last_edge = i;
+            last_edge = e;
           }
           if (--cnt == 0) {
             if (_curr_length > limit) break;
@@ -479,11 +496,12 @@
           }
         }
         if (_curr_length <= limit) {
-          for (int i = 0; i < _next_edge; ++i) {
-            e = _edges[i];
-            if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+          for (int e = 0; e < _next_edge; ++e) {
+            _cand_cost[e] = _state[e] *
+              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+            if (_cand_cost[e] < 0) {
               _candidates[_curr_length++] = e;
-              last_edge = i;
+              last_edge = e;
             }
             if (--cnt == 0) {
               if (_curr_length > limit) break;
@@ -500,12 +518,13 @@
                    _sort_func );
 
         // Pop the first element of the heap
-        _ns._in_edge = _candidates[0];
+        _in_edge = _candidates[0];
         pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
                   _sort_func );
         _curr_length = std::min(_head_length, _curr_length - 1);
         return true;
       }
+
     }; //class AlteringListPivotRule
 
   private:
@@ -519,66 +538,87 @@
 
   private:
 
-    // The directed graph the algorithm runs on
-    SGraph _graph;
     // The original graph
-    const Graph &_graph_ref;
+    const Graph &_orig_graph;
     // The original lower bound map
-    const LowerMap *_lower;
-    // The capacity map
-    SCapacityMap _capacity;
-    // The cost map
-    SCostMap _cost;
-    // The supply map
-    SSupplyMap _supply;
-    bool _valid_supply;
+    const LowerMap *_orig_lower;
+    // The original capacity map
+    const CapacityMap &_orig_cap;
+    // The original cost map
+    const CostMap &_orig_cost;
+    // The original supply map
+    const SupplyMap *_orig_supply;
+    // The source node (if no supply map was given)
+    Node _orig_source;
+    // The target node (if no supply map was given)
+    Node _orig_target;
+    // The flow value (if no supply map was given)
+    Capacity _orig_flow_value;
 
-    // Edge map of the current flow
-    SCapacityMap _flow;
-    // Node map of the current potentials
-    SPotentialMap _potential;
+    // The flow result map
+    FlowMap *_flow_result;
+    // The potential result map
+    PotentialMap *_potential_result;
+    // Indicate if the flow result map is local
+    bool _local_flow;
+    // Indicate if the potential result map is local
+    bool _local_potential;
 
-    // The depth node map of the spanning tree structure
-    IntNodeMap _depth;
-    // The parent node map of the spanning tree structure
-    NodeNodeMap _parent;
-    // The pred_edge node map of the spanning tree structure
-    EdgeNodeMap _pred_edge;
-    // The thread node map of the spanning tree structure
-    NodeNodeMap _thread;
-    // The forward node map of the spanning tree structure
-    BoolNodeMap _forward;
-    // The state edge map
-    IntEdgeMap _state;
-    // The artificial root node of the spanning tree structure
-    Node _root;
+    // The edge references
+    EdgeVector _edge;
+    // The node references
+    NodeVector _node;
+    // The node ids
+    IntNodeMap _node_id;
+    // The source nodes of the edges
+    IntVector _source;
+    // The target nodess of the edges
+    IntVector _target;
 
-    // The reduced cost map
-    ReducedCostMap _red_cost;
+    // The (modified) capacity vector
+    CapacityVector _cap;
+    // The cost vector
+    CostVector _cost;
+    // The (modified) supply vector
+    CostVector _supply;
+    // The current flow vector
+    CapacityVector _flow;
+    // The current potential vector
+    CostVector _pi;
 
-    // The non-artifical edges
-    EdgeVector _edges;
+    // The number of nodes in the original graph
+    int _node_num;
+    // The number of edges in the original graph
+    int _edge_num;
 
-    // Members for handling the original graph
-    FlowMap *_flow_result;
-    PotentialMap *_potential_result;
-    bool _local_flow;
-    bool _local_potential;
-    NodeRefMap _node_ref;
-    EdgeRefMap _edge_ref;
+    // The depth vector of the spanning tree structure
+    IntVector _depth;
+    // The parent vector of the spanning tree structure
+    IntVector _parent;
+    // The pred_edge vector of the spanning tree structure
+    IntVector _pred;
+    // The thread vector of the spanning tree structure
+    IntVector _thread;
+    // The forward vector of the spanning tree structure
+    BoolVector _forward;
+    // The state vector
+    IntVector _state;
+    // The root node
+    int _root;
 
     // The entering edge of the current pivot iteration
-    Edge _in_edge;
+    int _in_edge;
 
     // Temporary nodes used in the current pivot iteration
-    Node join, u_in, v_in, u_out, v_out;
-    Node right, first, second, last;
-    Node stem, par_stem, new_stem;
+    int join, u_in, v_in, u_out, v_out;
+    int right, first, second, last;
+    int stem, par_stem, new_stem;
+
     // The maximum augment amount along the found cycle in the current
     // pivot iteration
     Capacity delta;
 
-  public :
+  public:
 
     /// \brief General constructor (with lower bounds).
     ///
@@ -594,41 +634,12 @@
                     const CapacityMap &capacity,
                     const CostMap &cost,
                     const SupplyMap &supply ) :
-      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
-      _cost(_graph), _supply(_graph), _flow(_graph),
-      _potential(_graph), _depth(_graph), _parent(_graph),
-      _pred_edge(_graph), _thread(_graph), _forward(_graph),
-      _state(_graph), _red_cost(_graph, _cost, _potential),
+      _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
+      _orig_cost(cost), _orig_supply(&supply),
       _flow_result(NULL), _potential_result(NULL),
       _local_flow(false), _local_potential(false),
-      _node_ref(graph), _edge_ref(graph)
-    {
-      // Check the sum of supply values
-      Supply sum = 0;
-      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
-        sum += supply[n];
-      if (!(_valid_supply = sum == 0)) return;
-
-      // Copy _graph_ref to _graph
-      _graph.reserveNode(countNodes(_graph_ref) + 1);
-      _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
-      copyGraph(_graph, _graph_ref)
-        .edgeMap(_capacity, capacity)
-        .edgeMap(_cost, cost)
-        .nodeMap(_supply, supply)
-        .nodeRef(_node_ref)
-        .edgeRef(_edge_ref)
-        .run();
-
-      // Remove non-zero lower bounds
-      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
-        if (lower[e] != 0) {
-        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
-          _supply[_node_ref[_graph_ref.source(e)]] -= lower[e];
-          _supply[_node_ref[_graph_ref.target(e)]] += lower[e];
-      }
-      }
-    }
+      _node_id(graph)
+    {}
 
     /// \brief General constructor (without lower bounds).
     ///
@@ -642,32 +653,12 @@
                     const CapacityMap &capacity,
                     const CostMap &cost,
                     const SupplyMap &supply ) :
-      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
-      _cost(_graph), _supply(_graph), _flow(_graph),
-      _potential(_graph), _depth(_graph), _parent(_graph),
-      _pred_edge(_graph), _thread(_graph), _forward(_graph),
-      _state(_graph), _red_cost(_graph, _cost, _potential),
+      _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
+      _orig_cost(cost), _orig_supply(&supply),
       _flow_result(NULL), _potential_result(NULL),
       _local_flow(false), _local_potential(false),
-      _node_ref(graph), _edge_ref(graph)
-    {
-      // Check the sum of supply values
-      Supply sum = 0;
-      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
-        sum += supply[n];
-      if (!(_valid_supply = sum == 0)) return;
-
-      // Copy _graph_ref to _graph
-      _graph.reserveNode(countNodes(_graph_ref) + 1);
-      _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
-      copyGraph(_graph, _graph_ref)
-        .edgeMap(_capacity, capacity)
-        .edgeMap(_cost, cost)
-        .nodeMap(_supply, supply)
-        .nodeRef(_node_ref)
-        .edgeRef(_edge_ref)
-        .run();
-    }
+      _node_id(graph)
+    {}
 
     /// \brief Simple constructor (with lower bounds).
     ///
@@ -685,40 +676,15 @@
                     const LowerMap &lower,
                     const CapacityMap &capacity,
                     const CostMap &cost,
-                    typename Graph::Node s,
-                    typename Graph::Node t,
-                    typename SupplyMap::Value flow_value ) :
-      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
-      _cost(_graph), _supply(_graph, 0), _flow(_graph),
-      _potential(_graph), _depth(_graph), _parent(_graph),
-      _pred_edge(_graph), _thread(_graph), _forward(_graph),
-      _state(_graph), _red_cost(_graph, _cost, _potential),
+                    Node s, Node t,
+                    Capacity flow_value ) :
+      _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
+      _orig_cost(cost), _orig_supply(NULL),
+      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
       _flow_result(NULL), _potential_result(NULL),
       _local_flow(false), _local_potential(false),
-      _node_ref(graph), _edge_ref(graph)
-    {
-      // Copy _graph_ref to graph
-      _graph.reserveNode(countNodes(_graph_ref) + 1);
-      _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
-      copyGraph(_graph, _graph_ref)
-        .edgeMap(_capacity, capacity)
-        .edgeMap(_cost, cost)
-        .nodeRef(_node_ref)
-        .edgeRef(_edge_ref)
-        .run();
-
-      // Remove non-zero lower bounds
-      _supply[_node_ref[s]] =  flow_value;
-      _supply[_node_ref[t]] = -flow_value;
-      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
-        if (lower[e] != 0) {
-        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
-          _supply[_node_ref[_graph_ref.source(e)]] -= lower[e];
-          _supply[_node_ref[_graph_ref.target(e)]] += lower[e];
-      }
-      }
-      _valid_supply = true;
-    }
+      _node_id(graph)
+    {}
 
     /// \brief Simple constructor (without lower bounds).
     ///
@@ -734,31 +700,15 @@
     NetworkSimplex( const Graph &graph,
                     const CapacityMap &capacity,
                     const CostMap &cost,
-                    typename Graph::Node s,
-                    typename Graph::Node t,
-                    typename SupplyMap::Value flow_value ) :
-      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
-      _cost(_graph), _supply(_graph, 0), _flow(_graph),
-      _potential(_graph), _depth(_graph), _parent(_graph),
-      _pred_edge(_graph), _thread(_graph), _forward(_graph),
-      _state(_graph), _red_cost(_graph, _cost, _potential),
+                    Node s, Node t,
+                    Capacity flow_value ) :
+      _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
+      _orig_cost(cost), _orig_supply(NULL),
+      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
       _flow_result(NULL), _potential_result(NULL),
       _local_flow(false), _local_potential(false),
-      _node_ref(graph), _edge_ref(graph)
-    {
-      // Copy _graph_ref to graph
-      _graph.reserveNode(countNodes(_graph_ref) + 1);
-      _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
-      copyGraph(_graph, _graph_ref)
-        .edgeMap(_capacity, capacity)
-        .edgeMap(_cost, cost)
-        .nodeRef(_node_ref)
-        .edgeRef(_edge_ref)
-        .run();
-      _supply[_node_ref[s]] =  flow_value;
-      _supply[_node_ref[t]] = -flow_value;
-      _valid_supply = true;
-    }
+      _node_id(graph)
+    {}
 
     /// Destructor.
     ~NetworkSimplex() {
@@ -894,8 +844,8 @@
     /// \pre \ref run() must be called before using this function.
     Cost totalCost() const {
       Cost c = 0;
-      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
-        c += (*_flow_result)[e] * _cost[_edge_ref[e]];
+      for (EdgeIt e(_orig_graph); e != INVALID; ++e)
+        c += (*_flow_result)[e] * _orig_cost[e];
       return c;
     }
 
@@ -903,153 +853,216 @@
 
   private:
 
-    // Extend the underlying graph and initialize all the node and edge maps
+    // Initialize internal data structures
     bool init() {
-      if (!_valid_supply) return false;
-
       // Initialize result maps
       if (!_flow_result) {
-        _flow_result = new FlowMap(_graph_ref);
+        _flow_result = new FlowMap(_orig_graph);
         _local_flow = true;
       }
       if (!_potential_result) {
-        _potential_result = new PotentialMap(_graph_ref);
+        _potential_result = new PotentialMap(_orig_graph);
         _local_potential = true;
       }
-
-      // Initialize the edge vector and the edge maps
-      _edges.reserve(countEdges(_graph));
-      for (EdgeIt e(_graph); e != INVALID; ++e) {
-        _edges.push_back(e);
-        _flow[e] = 0;
-        _state[e] = STATE_LOWER;
+      
+      // Initialize vectors
+      _node_num = countNodes(_orig_graph);
+      _edge_num = countEdges(_orig_graph);
+      int all_node_num = _node_num + 1;
+      int all_edge_num = _edge_num + _node_num;
+      
+      _edge.resize(_edge_num);
+      _node.reserve(_node_num);
+      _source.resize(all_edge_num);
+      _target.resize(all_edge_num);
+      
+      _cap.resize(all_edge_num);
+      _cost.resize(all_edge_num);
+      _supply.resize(all_node_num);
+      _flow.resize(all_edge_num, 0);
+      _pi.resize(all_node_num, 0);
+      
+      _depth.resize(all_node_num);
+      _parent.resize(all_node_num);
+      _pred.resize(all_node_num);
+      _thread.resize(all_node_num);
+      _forward.resize(all_node_num);
+      _state.resize(all_edge_num, STATE_LOWER);
+      
+      // Initialize node related data
+      bool valid_supply = true;
+      if (_orig_supply) {
+        Supply sum = 0;
+        int i = 0;
+        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
+          _node.push_back(n);
+          _node_id[n] = i;
+          _supply[i] = (*_orig_supply)[n];
+          sum += _supply[i];
+        }
+        valid_supply = (sum == 0);
+      } else {
+        int i = 0;
+        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
+          _node.push_back(n);
+          _node_id[n] = i;
+          _supply[i] = 0;
+        }
+        _supply[_node_id[_orig_source]] =  _orig_flow_value;
+        _supply[_node_id[_orig_target]] = -_orig_flow_value;
+      }
+      if (!valid_supply) return false;
+      
+      // Set data for the artificial root node
+      _root = _node_num;
+      _depth[_root] = 0;
+      _parent[_root] = -1;
+      _pred[_root] = -1;
+      _thread[_root] = 0;
+      _supply[_root] = 0;
+      _pi[_root] = 0;
+      
+      // Store the edges in a mixed order
+      int k = std::max(int(sqrt(_edge_num)), 10);
+      int i = 0;
+      for (EdgeIt e(_orig_graph); e != INVALID; ++e) {
+        _edge[i] = e;
+        if ((i += k) >= _edge_num) i = (i % k) + 1;
       }
 
-      // Add an artificial root node to the graph
-      _root = _graph.addNode();
-      _parent[_root] = INVALID;
-      _pred_edge[_root] = INVALID;
-      _depth[_root] = 0;
-      _supply[_root] = 0;
-      _potential[_root] = 0;
+      // Initialize edge maps
+      for (int i = 0; i != _edge_num; ++i) {
+        Edge e = _edge[i];
+        _source[i] = _node_id[_orig_graph.source(e)];
+        _target[i] = _node_id[_orig_graph.target(e)];
+        _cost[i] = _orig_cost[e];
+        _cap[i] = _orig_cap[e];
+      }
 
-      // Add artificial edges to the graph and initialize the node maps
-      // of the spanning tree data structure
-      Node last = _root;
-      Edge e;
+      // Remove non-zero lower bounds
+      if (_orig_lower) {
+        for (int i = 0; i != _edge_num; ++i) {
+          Capacity c = (*_orig_lower)[_edge[i]];
+          if (c != 0) {
+            _cap[i] -= c;
+            _supply[_source[i]] -= c;
+            _supply[_target[i]] += c;
+          }
+        }
+      }
+
+      // Add artificial edges and initialize the spanning tree data structure
       Cost max_cost = std::numeric_limits<Cost>::max() / 4;
-      for (NodeIt u(_graph); u != INVALID; ++u) {
-        if (u == _root) continue;
-        _thread[last] = u;
-        last = u;
+      Capacity max_cap = std::numeric_limits<Capacity>::max();
+      for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
+        _thread[u] = u + 1;
+        _depth[u] = 1;
         _parent[u] = _root;
-        _depth[u] = 1;
+        _pred[u] = e;
         if (_supply[u] >= 0) {
-          e = _graph.addEdge(u, _root);
           _flow[e] = _supply[u];
           _forward[u] = true;
-          _potential[u] = -max_cost;
+          _pi[u] = -max_cost;
         } else {
-          e = _graph.addEdge(_root, u);
           _flow[e] = -_supply[u];
           _forward[u] = false;
-          _potential[u] = max_cost;
+          _pi[u] = max_cost;
         }
         _cost[e] = max_cost;
-        _capacity[e] = std::numeric_limits<Capacity>::max();
+        _cap[e] = max_cap;
         _state[e] = STATE_TREE;
-        _pred_edge[u] = e;
       }
-      _thread[last] = _root;
 
       return true;
     }
 
     // Find the join node
     void findJoinNode() {
-      Node u = _graph.source(_in_edge);
-      Node v = _graph.target(_in_edge);
+      int u = _source[_in_edge];
+      int v = _target[_in_edge];
+      while (_depth[u] > _depth[v]) u = _parent[u];
+      while (_depth[v] > _depth[u]) v = _parent[v];
       while (u != v) {
-        if (_depth[u] == _depth[v]) {
-          u = _parent[u];
-          v = _parent[v];
-        }
-        else if (_depth[u] > _depth[v]) u = _parent[u];
-        else v = _parent[v];
+        u = _parent[u];
+        v = _parent[v];
       }
       join = u;
     }
 
-    // Find the leaving edge of the cycle and returns true if the 
+    // Find the leaving edge of the cycle and returns true if the
     // leaving edge is not the same as the entering edge
     bool findLeavingEdge() {
       // Initialize first and second nodes according to the direction
       // of the cycle
       if (_state[_in_edge] == STATE_LOWER) {
-        first  = _graph.source(_in_edge);
-        second = _graph.target(_in_edge);
+        first  = _source[_in_edge];
+        second = _target[_in_edge];
       } else {
-        first  = _graph.target(_in_edge);
-        second = _graph.source(_in_edge);
+        first  = _target[_in_edge];
+        second = _source[_in_edge];
       }
-      delta = _capacity[_in_edge];
-      bool result = false;
+      delta = _cap[_in_edge];
+      int result = 0;
       Capacity d;
-      Edge e;
+      int e;
 
       // Search the cycle along the path form the first node to the root
-      for (Node u = first; u != join; u = _parent[u]) {
-        e = _pred_edge[u];
-        d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
+      for (int u = first; u != join; u = _parent[u]) {
+        e = _pred[u];
+        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
         if (d < delta) {
           delta = d;
           u_out = u;
-          u_in = first;
-          v_in = second;
-          result = true;
+          result = 1;
         }
       }
       // Search the cycle along the path form the second node to the root
-      for (Node u = second; u != join; u = _parent[u]) {
-        e = _pred_edge[u];
-        d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
+      for (int u = second; u != join; u = _parent[u]) {
+        e = _pred[u];
+        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
         if (d <= delta) {
           delta = d;
           u_out = u;
-          u_in = second;
-          v_in = first;
-          result = true;
+          result = 2;
         }
       }
-      return result;
+
+      if (result == 1) {
+        u_in = first;
+        v_in = second;
+      } else {
+        u_in = second;
+        v_in = first;
+      }
+      return result != 0;
     }
 
-    // Change _flow and _state edge maps
-    void changeFlows(bool change) {
+    // Change _flow and _state vectors
+    void changeFlow(bool change) {
       // Augment along the cycle
       if (delta > 0) {
         Capacity val = _state[_in_edge] * delta;
         _flow[_in_edge] += val;
-        for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
-          _flow[_pred_edge[u]] += _forward[u] ? -val : val;
+        for (int u = _source[_in_edge]; u != join; u = _parent[u]) {
+          _flow[_pred[u]] += _forward[u] ? -val : val;
         }
-        for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
-          _flow[_pred_edge[u]] += _forward[u] ? val : -val;
+        for (int u = _target[_in_edge]; u != join; u = _parent[u]) {
+          _flow[_pred[u]] += _forward[u] ? val : -val;
         }
       }
       // Update the state of the entering and leaving edges
       if (change) {
         _state[_in_edge] = STATE_TREE;
-        _state[_pred_edge[u_out]] =
-          (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
+        _state[_pred[u_out]] =
+          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
       } else {
         _state[_in_edge] = -_state[_in_edge];
       }
     }
 
-    // Update _thread and _parent node maps
+    // Update _thread and _parent vectors
     void updateThreadParent() {
-      Node u;
+      int u;
       v_out = _parent[u_out];
 
       // Handle the case when join and v_out coincide
@@ -1105,30 +1118,30 @@
       }
     }
 
-    // Update _pred_edge and _forward node maps.
+    // Update _pred and _forward vectors
     void updatePredEdge() {
-      Node u = u_out, v;
+      int u = u_out, v;
       while (u != u_in) {
         v = _parent[u];
-        _pred_edge[u] = _pred_edge[v];
+        _pred[u] = _pred[v];
         _forward[u] = !_forward[v];
         u = v;
       }
-      _pred_edge[u_in] = _in_edge;
-      _forward[u_in] = (u_in == _graph.source(_in_edge));
+      _pred[u_in] = _in_edge;
+      _forward[u_in] = (u_in == _source[_in_edge]);
     }
 
-    // Update _depth and _potential node maps
+    // Update _depth and _potential vectors
     void updateDepthPotential() {
       _depth[u_in] = _depth[v_in] + 1;
       Cost sigma = _forward[u_in] ?
-        _potential[v_in] - _potential[u_in] - _cost[_pred_edge[u_in]] :
-        _potential[v_in] - _potential[u_in] + _cost[_pred_edge[u_in]];
-      _potential[u_in] += sigma;
-      for(Node u = _thread[u_in]; _parent[u] != INVALID; u = _thread[u]) {
+        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
+        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
+      _pi[u_in] += sigma;
+      for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
         _depth[u] = _depth[_parent[u]] + 1;
         if (_depth[u] <= _depth[u_in]) break;
-        _potential[u] += sigma;
+        _pi[u] += sigma;
       }
     }
 
@@ -1152,13 +1165,13 @@
 
     template<class PivotRuleImplementation>
     bool start() {
-      PivotRuleImplementation pivot(*this, _edges);
+      PivotRuleImplementation pivot(*this);
 
       // Execute the network simplex algorithm
       while (pivot.findEnteringEdge()) {
         findJoinNode();
         bool change = findLeavingEdge();
-        changeFlows(change);
+        changeFlow(change);
         if (change) {
           updateThreadParent();
           updatePredEdge();
@@ -1167,22 +1180,25 @@
       }
 
       // Check if the flow amount equals zero on all the artificial edges
-      for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
+      for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
         if (_flow[e] > 0) return false;
-      for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
-        if (_flow[e] > 0) return false;
+      }
 
       // Copy flow values to _flow_result
-      if (_lower) {
-        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
-          (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
+      if (_orig_lower) {
+        for (int i = 0; i != _edge_num; ++i) {
+          Edge e = _edge[i];
+          (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
+        }
       } else {
-        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
-          (*_flow_result)[e] = _flow[_edge_ref[e]];
+        for (int i = 0; i != _edge_num; ++i) {
+          (*_flow_result)[_edge[i]] = _flow[i];
+        }
       }
       // Copy potential values to _potential_result
-      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
-        (*_potential_result)[n] = _potential[_node_ref[n]];
+      for (int i = 0; i != _node_num; ++i) {
+        (*_potential_result)[_node[i]] = _pi[i];
+      }
 
       return true;
     }