COIN-OR::LEMON - Graph Library

Ticket #340: 340-mcf-new-heur-87ff083a3f6b.patch

File 340-mcf-new-heur-87ff083a3f6b.patch, 27.6 KB (added by Peter Kovacs, 10 years ago)
  • lemon/capacity_scaling.h

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1265875904 -3600
    # Node ID 87ff083a3f6bcd65e4de7647134ec64a9c9e2554
    # Parent  a7e93de12cbda2267756b130476b8e84572002bf
    New heuristics for MCF algorithms (#340)
    and some implementation improvements.
    
     - A simple heuristic is added to NetworkSimplex to make the
       initial pivots faster.
     - A powerful global update heuristic is added to CostScaling
       and the implementation is reworked with various improvements.
     - A small improvement is made in CapacityScaling for better
       delta computation.
    
    diff --git a/lemon/capacity_scaling.h b/lemon/capacity_scaling.h
    a b  
    764764      // Initialize delta value
    765765      if (_factor > 1) {
    766766        // With scaling
    767         Value max_sup = 0, max_dem = 0;
    768         for (int i = 0; i != _node_num; ++i) {
     767        Value max_sup = 0, max_dem = 0, max_cap = 0;
     768        for (int i = 0; i != _root; ++i) {
    769769          Value ex = _excess[i];
    770770          if ( ex > max_sup) max_sup =  ex;
    771771          if (-ex > max_dem) max_dem = -ex;
    772         }
    773         Value max_cap = 0;
    774         for (int j = 0; j != _res_arc_num; ++j) {
    775           if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
     772          int last_out = _first_out[i+1] - 1;
     773          for (int j = _first_out[i]; j != last_out; ++j) {
     774            if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
     775          }
    776776        }
    777777        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
    778778        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
  • lemon/cost_scaling.h

    diff --git a/lemon/cost_scaling.h b/lemon/cost_scaling.h
    a b  
    197197    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    198198
    199199    typedef std::vector<int> IntVector;
    200     typedef std::vector<char> BoolVector;
     200    typedef std::vector<char> CharVector;
    201201    typedef std::vector<Value> ValueVector;
    202202    typedef std::vector<Cost> CostVector;
    203203    typedef std::vector<LargeCost> LargeCostVector;
     
    244244    // Parameters of the problem
    245245    bool _have_lower;
    246246    Value _sum_supply;
     247    int _sup_node_num;
    247248
    248249    // Data structures for storing the digraph
    249250    IntNodeMap _node_id;
    250251    IntArcMap _arc_idf;
    251252    IntArcMap _arc_idb;
    252253    IntVector _first_out;
    253     BoolVector _forward;
     254    CharVector _forward;
    254255    IntVector _source;
    255256    IntVector _target;
    256257    IntVector _reverse;
     
    272273    LargeCost _epsilon;
    273274    int _alpha;
    274275
     276    IntVector _buckets;
     277    IntVector _bucket_next;
     278    IntVector _bucket_prev;
     279    IntVector _rank;
     280    int _max_rank;
     281 
    275282    // Data for a StaticDigraph structure
    276283    typedef std::pair<int, int> IntPair;
    277284    StaticDigraph _sgr;
     
    802809        }
    803810      }
    804811
     812      _sup_node_num = 0;
     813      for (NodeIt n(_graph); n != INVALID; ++n) {
     814        if (sup[n] > 0) ++_sup_node_num;
     815      }
     816
    805817      // Find a feasible flow using Circulation
    806818      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
    807819        circ(_graph, low, cap, sup);
     
    836848        }
    837849        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
    838850          int ra = _reverse[a];
    839           _res_cap[a] = 1;
     851          _res_cap[a] = 0;
    840852          _res_cap[ra] = 0;
    841853          _cost[a] = 0;
    842854          _cost[ra] = 0;
     
    850862    void start(Method method) {
    851863      // Maximum path length for partial augment
    852864      const int MAX_PATH_LENGTH = 4;
    853      
     865
     866      // Initialize data structures for buckets     
     867      _max_rank = _alpha * _res_node_num;
     868      _buckets.resize(_max_rank);
     869      _bucket_next.resize(_res_node_num + 1);
     870      _bucket_prev.resize(_res_node_num + 1);
     871      _rank.resize(_res_node_num + 1);
     872 
    854873      // Execute the algorithm
    855874      switch (method) {
    856875        case PUSH:
     
    889908        }
    890909      }
    891910    }
     911   
     912    // Initialize a cost scaling phase
     913    void initPhase() {
     914      // Saturate arcs not satisfying the optimality condition
     915      for (int u = 0; u != _res_node_num; ++u) {
     916        int last_out = _first_out[u+1];
     917        LargeCost pi_u = _pi[u];
     918        for (int a = _first_out[u]; a != last_out; ++a) {
     919          int v = _target[a];
     920          if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
     921            Value delta = _res_cap[a];
     922            _excess[u] -= delta;
     923            _excess[v] += delta;
     924            _res_cap[a] = 0;
     925            _res_cap[_reverse[a]] += delta;
     926          }
     927        }
     928      }
     929     
     930      // Find active nodes (i.e. nodes with positive excess)
     931      for (int u = 0; u != _res_node_num; ++u) {
     932        if (_excess[u] > 0) _active_nodes.push_back(u);
     933      }
     934
     935      // Initialize the next arcs
     936      for (int u = 0; u != _res_node_num; ++u) {
     937        _next_out[u] = _first_out[u];
     938      }
     939    }
     940   
     941    // Early termination heuristic
     942    bool earlyTermination() {
     943      const double EARLY_TERM_FACTOR = 3.0;
     944
     945      // Build a static residual graph
     946      _arc_vec.clear();
     947      _cost_vec.clear();
     948      for (int j = 0; j != _res_arc_num; ++j) {
     949        if (_res_cap[j] > 0) {
     950          _arc_vec.push_back(IntPair(_source[j], _target[j]));
     951          _cost_vec.push_back(_cost[j] + 1);
     952        }
     953      }
     954      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
     955
     956      // Run Bellman-Ford algorithm to check if the current flow is optimal
     957      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
     958      bf.init(0);
     959      bool done = false;
     960      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
     961      for (int i = 0; i < K && !done; ++i) {
     962        done = bf.processNextWeakRound();
     963      }
     964      return done;
     965    }
     966
     967    // Global potential update heuristic
     968    void globalUpdate() {
     969      int bucket_end = _root + 1;
     970   
     971      // Initialize buckets
     972      for (int r = 0; r != _max_rank; ++r) {
     973        _buckets[r] = bucket_end;
     974      }
     975      Value total_excess = 0;
     976      for (int i = 0; i != _res_node_num; ++i) {
     977        if (_excess[i] < 0) {
     978          _rank[i] = 0;
     979          _bucket_next[i] = _buckets[0];
     980          _bucket_prev[_buckets[0]] = i;
     981          _buckets[0] = i;
     982        } else {
     983          total_excess += _excess[i];
     984          _rank[i] = _max_rank;
     985        }
     986      }
     987      if (total_excess == 0) return;
     988
     989      // Search the buckets
     990      int r = 0;
     991      for ( ; r != _max_rank; ++r) {
     992        while (_buckets[r] != bucket_end) {
     993          // Remove the first node from the current bucket
     994          int u = _buckets[r];
     995          _buckets[r] = _bucket_next[u];
     996         
     997          // Search the incomming arcs of u
     998          LargeCost pi_u = _pi[u];
     999          int last_out = _first_out[u+1];
     1000          for (int a = _first_out[u]; a != last_out; ++a) {
     1001            int ra = _reverse[a];
     1002            if (_res_cap[ra] > 0) {
     1003              int v = _source[ra];
     1004              int old_rank_v = _rank[v];
     1005              if (r < old_rank_v) {
     1006                // Compute the new rank of v
     1007                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
     1008                int new_rank_v = old_rank_v;
     1009                if (nrc < LargeCost(_max_rank))
     1010                  new_rank_v = r + 1 + int(nrc);
     1011                 
     1012                // Change the rank of v
     1013                if (new_rank_v < old_rank_v) {
     1014                  _rank[v] = new_rank_v;
     1015                  _next_out[v] = _first_out[v];
     1016                 
     1017                  // Remove v from its old bucket
     1018                  if (old_rank_v < _max_rank) {
     1019                    if (_buckets[old_rank_v] == v) {
     1020                      _buckets[old_rank_v] = _bucket_next[v];
     1021                    } else {
     1022                      _bucket_next[_bucket_prev[v]] = _bucket_next[v];
     1023                      _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
     1024                    }
     1025                  }
     1026                 
     1027                  // Insert v to its new bucket
     1028                  _bucket_next[v] = _buckets[new_rank_v];
     1029                  _bucket_prev[_buckets[new_rank_v]] = v;
     1030                  _buckets[new_rank_v] = v;
     1031                }
     1032              }
     1033            }
     1034          }
     1035
     1036          // Finish search if there are no more active nodes
     1037          if (_excess[u] > 0) {
     1038            total_excess -= _excess[u];
     1039            if (total_excess <= 0) break;
     1040          }
     1041        }
     1042        if (total_excess <= 0) break;
     1043      }
     1044     
     1045      // Relabel nodes
     1046      for (int u = 0; u != _res_node_num; ++u) {
     1047        int k = std::min(_rank[u], r);
     1048        if (k > 0) {
     1049          _pi[u] -= _epsilon * k;
     1050          _next_out[u] = _first_out[u];
     1051        }
     1052      }
     1053    }
    8921054
    8931055    /// Execute the algorithm performing augment and relabel operations
    8941056    void startAugment(int max_length = std::numeric_limits<int>::max()) {
    8951057      // Paramters for heuristics
    896       const int BF_HEURISTIC_EPSILON_BOUND = 1000;
    897       const int BF_HEURISTIC_BOUND_FACTOR  = 3;
     1058      const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1059      const double GLOBAL_UPDATE_FACTOR = 3.0;
    8981060
     1061      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1062        (_res_node_num + _sup_node_num * _sup_node_num));
     1063      int next_update_limit = global_update_freq;
     1064     
     1065      int relabel_cnt = 0;
     1066     
    8991067      // Perform cost scaling phases
    900       IntVector pred_arc(_res_node_num);
    901       std::vector<int> path_nodes;
     1068      std::vector<int> path;
    9021069      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    9031070                                        1 : _epsilon / _alpha )
    9041071      {
    905         // "Early Termination" heuristic: use Bellman-Ford algorithm
    906         // to check if the current flow is optimal
    907         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
    908           _arc_vec.clear();
    909           _cost_vec.clear();
    910           for (int j = 0; j != _res_arc_num; ++j) {
    911             if (_res_cap[j] > 0) {
    912               _arc_vec.push_back(IntPair(_source[j], _target[j]));
    913               _cost_vec.push_back(_cost[j] + 1);
    914             }
    915           }
    916           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    917 
    918           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
    919           bf.init(0);
    920           bool done = false;
    921           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
    922           for (int i = 0; i < K && !done; ++i)
    923             done = bf.processNextWeakRound();
    924           if (done) break;
    925         }
    926 
    927         // Saturate arcs not satisfying the optimality condition
    928         for (int a = 0; a != _res_arc_num; ++a) {
    929           if (_res_cap[a] > 0 &&
    930               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
    931             Value delta = _res_cap[a];
    932             _excess[_source[a]] -= delta;
    933             _excess[_target[a]] += delta;
    934             _res_cap[a] = 0;
    935             _res_cap[_reverse[a]] += delta;
    936           }
     1072        // Early termination heuristic
     1073        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
     1074          if (earlyTermination()) break;
    9371075        }
    9381076       
    939         // Find active nodes (i.e. nodes with positive excess)
    940         for (int u = 0; u != _res_node_num; ++u) {
    941           if (_excess[u] > 0) _active_nodes.push_back(u);
    942         }
    943 
    944         // Initialize the next arcs
    945         for (int u = 0; u != _res_node_num; ++u) {
    946           _next_out[u] = _first_out[u];
    947         }
    948 
     1077        // Initialize current phase
     1078        initPhase();
     1079       
    9491080        // Perform partial augment and relabel operations
    9501081        while (true) {
    9511082          // Select an active node (FIFO selection)
     
    9551086          }
    9561087          if (_active_nodes.size() == 0) break;
    9571088          int start = _active_nodes.front();
    958           path_nodes.clear();
    959           path_nodes.push_back(start);
    9601089
    9611090          // Find an augmenting path from the start node
     1091          path.clear();
    9621092          int tip = start;
    963           while (_excess[tip] >= 0 &&
    964                  int(path_nodes.size()) <= max_length) {
     1093          while (_excess[tip] >= 0 && int(path.size()) < max_length) {
    9651094            int u;
    966             LargeCost min_red_cost, rc;
    967             int last_out = _sum_supply < 0 ?
    968               _first_out[tip+1] : _first_out[tip+1] - 1;
     1095            LargeCost min_red_cost, rc, pi_tip = _pi[tip];
     1096            int last_out = _first_out[tip+1];
    9691097            for (int a = _next_out[tip]; a != last_out; ++a) {
    970               if (_res_cap[a] > 0 &&
    971                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
    972                 u = _target[a];
    973                 pred_arc[u] = a;
     1098              u = _target[a];
     1099              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
     1100                path.push_back(a);
    9741101                _next_out[tip] = a;
    9751102                tip = u;
    976                 path_nodes.push_back(tip);
    9771103                goto next_step;
    9781104              }
    9791105            }
     
    9811107            // Relabel tip node
    9821108            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
    9831109            for (int a = _first_out[tip]; a != last_out; ++a) {
    984               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
     1110              rc = _cost[a] + pi_tip - _pi[_target[a]];
    9851111              if (_res_cap[a] > 0 && rc < min_red_cost) {
    9861112                min_red_cost = rc;
    9871113              }
    9881114            }
    9891115            _pi[tip] -= min_red_cost + _epsilon;
    990 
    991             // Reset the next arc of tip
    9921116            _next_out[tip] = _first_out[tip];
     1117            ++relabel_cnt;
    9931118
    9941119            // Step back
    9951120            if (tip != start) {
    996               path_nodes.pop_back();
    997               tip = path_nodes.back();
     1121              tip = _source[path.back()];
     1122              path.pop_back();
    9981123            }
    9991124
    10001125          next_step: ;
     
    10021127
    10031128          // Augment along the found path (as much flow as possible)
    10041129          Value delta;
    1005           int u, v = path_nodes.front(), pa;
    1006           for (int i = 1; i < int(path_nodes.size()); ++i) {
     1130          int pa, u, v = start;
     1131          for (int i = 0; i != int(path.size()); ++i) {
     1132            pa = path[i];
    10071133            u = v;
    1008             v = path_nodes[i];
    1009             pa = pred_arc[v];
     1134            v = _target[pa];
    10101135            delta = std::min(_res_cap[pa], _excess[u]);
    10111136            _res_cap[pa] -= delta;
    10121137            _res_cap[_reverse[pa]] += delta;
     
    10151140            if (_excess[v] > 0 && _excess[v] <= delta)
    10161141              _active_nodes.push_back(v);
    10171142          }
     1143
     1144          // Global update heuristic
     1145          if (relabel_cnt >= next_update_limit) {
     1146            globalUpdate();
     1147            next_update_limit += global_update_freq;
     1148          }
    10181149        }
    10191150      }
    10201151    }
     
    10221153    /// Execute the algorithm performing push and relabel operations
    10231154    void startPush() {
    10241155      // Paramters for heuristics
    1025       const int BF_HEURISTIC_EPSILON_BOUND = 1000;
    1026       const int BF_HEURISTIC_BOUND_FACTOR  = 3;
     1156      const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1157      const double GLOBAL_UPDATE_FACTOR = 2.0;
    10271158
     1159      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1160        (_res_node_num + _sup_node_num * _sup_node_num));
     1161      int next_update_limit = global_update_freq;
     1162
     1163      int relabel_cnt = 0;
     1164     
    10281165      // Perform cost scaling phases
    1029       BoolVector hyper(_res_node_num, false);
     1166      CharVector hyper(_res_node_num, false);
    10301167      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    10311168                                        1 : _epsilon / _alpha )
    10321169      {
    1033         // "Early Termination" heuristic: use Bellman-Ford algorithm
    1034         // to check if the current flow is optimal
    1035         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
    1036           _arc_vec.clear();
    1037           _cost_vec.clear();
    1038           for (int j = 0; j != _res_arc_num; ++j) {
    1039             if (_res_cap[j] > 0) {
    1040               _arc_vec.push_back(IntPair(_source[j], _target[j]));
    1041               _cost_vec.push_back(_cost[j] + 1);
    1042             }
    1043           }
    1044           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    1045 
    1046           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
    1047           bf.init(0);
    1048           bool done = false;
    1049           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
    1050           for (int i = 0; i < K && !done; ++i)
    1051             done = bf.processNextWeakRound();
    1052           if (done) break;
     1170        // Early termination heuristic
     1171        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
     1172          if (earlyTermination()) break;
    10531173        }
    1054 
    1055         // Saturate arcs not satisfying the optimality condition
    1056         for (int a = 0; a != _res_arc_num; ++a) {
    1057           if (_res_cap[a] > 0 &&
    1058               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
    1059             Value delta = _res_cap[a];
    1060             _excess[_source[a]] -= delta;
    1061             _excess[_target[a]] += delta;
    1062             _res_cap[a] = 0;
    1063             _res_cap[_reverse[a]] += delta;
    1064           }
    1065         }
    1066 
    1067         // Find active nodes (i.e. nodes with positive excess)
    1068         for (int u = 0; u != _res_node_num; ++u) {
    1069           if (_excess[u] > 0) _active_nodes.push_back(u);
    1070         }
    1071 
    1072         // Initialize the next arcs
    1073         for (int u = 0; u != _res_node_num; ++u) {
    1074           _next_out[u] = _first_out[u];
    1075         }
     1174       
     1175        // Initialize current phase
     1176        initPhase();
    10761177
    10771178        // Perform push and relabel operations
    10781179        while (_active_nodes.size() > 0) {
    1079           LargeCost min_red_cost, rc;
     1180          LargeCost min_red_cost, rc, pi_n;
    10801181          Value delta;
    10811182          int n, t, a, last_out = _res_arc_num;
    10821183
     1184        next_node:
    10831185          // Select an active node (FIFO selection)
    1084         next_node:
    10851186          n = _active_nodes.front();
    1086           last_out = _sum_supply < 0 ?
    1087             _first_out[n+1] : _first_out[n+1] - 1;
    1088 
     1187          last_out = _first_out[n+1];
     1188          pi_n = _pi[n];
     1189         
    10891190          // Perform push operations if there are admissible arcs
    10901191          if (_excess[n] > 0) {
    10911192            for (a = _next_out[n]; a != last_out; ++a) {
    10921193              if (_res_cap[a] > 0 &&
    1093                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
     1194                  _cost[a] + pi_n - _pi[_target[a]] < 0) {
    10941195                delta = std::min(_res_cap[a], _excess[n]);
    10951196                t = _target[a];
    10961197
    10971198                // Push-look-ahead heuristic
    10981199                Value ahead = -_excess[t];
    1099                 int last_out_t = _sum_supply < 0 ?
    1100                   _first_out[t+1] : _first_out[t+1] - 1;
     1200                int last_out_t = _first_out[t+1];
     1201                LargeCost pi_t = _pi[t];
    11011202                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
    11021203                  if (_res_cap[ta] > 0 &&
    1103                       _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
     1204                      _cost[ta] + pi_t - _pi[_target[ta]] < 0)
    11041205                    ahead += _res_cap[ta];
    11051206                  if (ahead >= delta) break;
    11061207                }
    11071208                if (ahead < 0) ahead = 0;
    11081209
    11091210                // Push flow along the arc
    1110                 if (ahead < delta) {
     1211                if (ahead < delta && !hyper[t]) {
    11111212                  _res_cap[a] -= ahead;
    11121213                  _res_cap[_reverse[a]] += ahead;
    11131214                  _excess[n] -= ahead;
     
    11381239          if (_excess[n] > 0 || hyper[n]) {
    11391240            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
    11401241            for (int a = _first_out[n]; a != last_out; ++a) {
    1141               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
     1242              rc = _cost[a] + pi_n - _pi[_target[a]];
    11421243              if (_res_cap[a] > 0 && rc < min_red_cost) {
    11431244                min_red_cost = rc;
    11441245              }
    11451246            }
    11461247            _pi[n] -= min_red_cost + _epsilon;
     1248            _next_out[n] = _first_out[n];
    11471249            hyper[n] = false;
    1148 
    1149             // Reset the next arc
    1150             _next_out[n] = _first_out[n];
     1250            ++relabel_cnt;
    11511251          }
    11521252       
    11531253          // Remove nodes that are not active nor hyper
     
    11571257                  !hyper[_active_nodes.front()] ) {
    11581258            _active_nodes.pop_front();
    11591259          }
     1260         
     1261          // Global update heuristic
     1262          if (relabel_cnt >= next_update_limit) {
     1263            globalUpdate();
     1264            for (int u = 0; u != _res_node_num; ++u)
     1265              hyper[u] = false;
     1266            next_update_limit += global_update_freq;
     1267          }
    11601268        }
    11611269      }
    11621270    }
  • lemon/network_simplex.h

    diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
    a b  
    265265      // Find next entering arc
    266266      bool findEnteringArc() {
    267267        Cost c;
    268         for (int e = _next_arc; e < _search_arc_num; ++e) {
     268        for (int e = _next_arc; e != _search_arc_num; ++e) {
    269269          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    270270          if (c < 0) {
    271271            _in_arc = e;
     
    273273            return true;
    274274          }
    275275        }
    276         for (int e = 0; e < _next_arc; ++e) {
     276        for (int e = 0; e != _next_arc; ++e) {
    277277          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    278278          if (c < 0) {
    279279            _in_arc = e;
     
    313313      // Find next entering arc
    314314      bool findEnteringArc() {
    315315        Cost c, min = 0;
    316         for (int e = 0; e < _search_arc_num; ++e) {
     316        for (int e = 0; e != _search_arc_num; ++e) {
    317317          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    318318          if (c < min) {
    319319            min = c;
     
    354354        _next_arc(0)
    355355      {
    356356        // The main parameters of the pivot rule
    357         const double BLOCK_SIZE_FACTOR = 0.5;
     357        const double BLOCK_SIZE_FACTOR = 1.0;
    358358        const int MIN_BLOCK_SIZE = 10;
    359359
    360360        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
     
    367367        Cost c, min = 0;
    368368        int cnt = _block_size;
    369369        int e;
    370         for (e = _next_arc; e < _search_arc_num; ++e) {
     370        for (e = _next_arc; e != _search_arc_num; ++e) {
    371371          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    372372          if (c < min) {
    373373            min = c;
     
    378378            cnt = _block_size;
    379379          }
    380380        }
    381         for (e = 0; e < _next_arc; ++e) {
     381        for (e = 0; e != _next_arc; ++e) {
    382382          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    383383          if (c < min) {
    384384            min = c;
     
    469469        // Major iteration: build a new candidate list
    470470        min = 0;
    471471        _curr_length = 0;
    472         for (e = _next_arc; e < _search_arc_num; ++e) {
     472        for (e = _next_arc; e != _search_arc_num; ++e) {
    473473          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    474474          if (c < 0) {
    475475            _candidates[_curr_length++] = e;
     
    480480            if (_curr_length == _list_length) goto search_end;
    481481          }
    482482        }
    483         for (e = 0; e < _next_arc; ++e) {
     483        for (e = 0; e != _next_arc; ++e) {
    484484          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    485485          if (c < 0) {
    486486            _candidates[_curr_length++] = e;
     
    564564      bool findEnteringArc() {
    565565        // Check the current candidate list
    566566        int e;
    567         for (int i = 0; i < _curr_length; ++i) {
     567        for (int i = 0; i != _curr_length; ++i) {
    568568          e = _candidates[i];
    569569          _cand_cost[e] = _state[e] *
    570570            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     
    577577        int cnt = _block_size;
    578578        int limit = _head_length;
    579579
    580         for (e = _next_arc; e < _search_arc_num; ++e) {
     580        for (e = _next_arc; e != _search_arc_num; ++e) {
    581581          _cand_cost[e] = _state[e] *
    582582            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    583583          if (_cand_cost[e] < 0) {
     
    589589            cnt = _block_size;
    590590          }
    591591        }
    592         for (e = 0; e < _next_arc; ++e) {
     592        for (e = 0; e != _next_arc; ++e) {
    593593          _cand_cost[e] = _state[e] *
    594594            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    595595          if (_cand_cost[e] < 0) {
     
    13281328      }
    13291329
    13301330      // Update _rev_thread using the new _thread values
    1331       for (int i = 0; i < int(_dirty_revs.size()); ++i) {
     1331      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
    13321332        u = _dirty_revs[i];
    13331333        _rev_thread[_thread[u]] = u;
    13341334      }
     
    14001400      }
    14011401    }
    14021402
     1403    // Heuristic initial pivots
     1404    bool initialPivots() {
     1405      Value curr, total = 0;
     1406      std::vector<Node> supply_nodes, demand_nodes;
     1407      for (NodeIt u(_graph); u != INVALID; ++u) {
     1408        curr = _supply[_node_id[u]];
     1409        if (curr > 0) {
     1410          total += curr;
     1411          supply_nodes.push_back(u);
     1412        }
     1413        else if (curr < 0) {
     1414          demand_nodes.push_back(u);
     1415        }
     1416      }
     1417      if (_sum_supply > 0) total -= _sum_supply;
     1418      if (total <= 0) return true;
     1419
     1420      IntVector arc_vector;
     1421      if (_sum_supply >= 0) {
     1422        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
     1423          // Perform a reverse graph search from the sink to the source
     1424          typename GR::template NodeMap<bool> reached(_graph, false);
     1425          Node s = supply_nodes[0], t = demand_nodes[0];
     1426          std::vector<Node> stack;
     1427          reached[t] = true;
     1428          stack.push_back(t);
     1429          while (!stack.empty()) {
     1430            Node u, v = stack.back();
     1431            stack.pop_back();
     1432            if (v == s) break;
     1433            for (InArcIt a(_graph, v); a != INVALID; ++a) {
     1434              if (reached[u = _graph.source(a)]) continue;
     1435              int j = _arc_id[a];
     1436              if (_cap[j] >= total) {
     1437                arc_vector.push_back(j);
     1438                reached[u] = true;
     1439                stack.push_back(u);
     1440              }
     1441            }
     1442          }
     1443        } else {
     1444          // Find the min. cost incomming arc for each demand node
     1445          for (int i = 0; i != int(demand_nodes.size()); ++i) {
     1446            Node v = demand_nodes[i];
     1447            Cost c, min_cost = std::numeric_limits<Cost>::max();
     1448            Arc min_arc = INVALID;
     1449            for (InArcIt a(_graph, v); a != INVALID; ++a) {
     1450              c = _cost[_arc_id[a]];
     1451              if (c < min_cost) {
     1452                min_cost = c;
     1453                min_arc = a;
     1454              }
     1455            }
     1456            if (min_arc != INVALID) {
     1457              arc_vector.push_back(_arc_id[min_arc]);
     1458            }
     1459          }
     1460        }
     1461      } else {
     1462        // Find the min. cost outgoing arc for each supply node
     1463        for (int i = 0; i != int(supply_nodes.size()); ++i) {
     1464          Node u = supply_nodes[i];
     1465          Cost c, min_cost = std::numeric_limits<Cost>::max();
     1466          Arc min_arc = INVALID;
     1467          for (OutArcIt a(_graph, u); a != INVALID; ++a) {
     1468            c = _cost[_arc_id[a]];
     1469            if (c < min_cost) {
     1470              min_cost = c;
     1471              min_arc = a;
     1472            }
     1473          }
     1474          if (min_arc != INVALID) {
     1475            arc_vector.push_back(_arc_id[min_arc]);
     1476          }
     1477        }
     1478      }
     1479
     1480      // Perform heuristic initial pivots
     1481      for (int i = 0; i != int(arc_vector.size()); ++i) {
     1482        in_arc = arc_vector[i];
     1483        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
     1484            _pi[_target[in_arc]]) >= 0) continue;
     1485        findJoinNode();
     1486        bool change = findLeavingArc();
     1487        if (delta >= MAX) return false;
     1488        changeFlow(change);
     1489        if (change) {
     1490          updateTreeStructure();
     1491          updatePotential();
     1492        }
     1493      }
     1494      return true;
     1495    }
     1496
    14031497    // Execute the algorithm
    14041498    ProblemType start(PivotRule pivot_rule) {
    14051499      // Select the pivot rule implementation
     
    14221516    ProblemType start() {
    14231517      PivotRuleImpl pivot(*this);
    14241518
     1519      // Perform heuristic initial pivots
     1520      if (!initialPivots()) return UNBOUNDED;
     1521
    14251522      // Execute the Network Simplex algorithm
    14261523      while (pivot.findEnteringArc()) {
    14271524        findJoinNode();