RhodeCode

        for (NodeIt n(_graph); n != INVALID; ++n) {

          _excess[_node_id[n]] = sup[n];

        }

        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {

          int u = _target[a];

          int ra = _reverse[a];

          _res_cap[a] = -_sum_supply + 1;

          _res_cap[ra] = -_excess[u];

          _cost[a] = 0;

          _cost[ra] = 0;

          _excess[u] = 0;

        }

      } else {

        for (ArcIt a(_graph); a != INVALID; ++a) {

          Value fa = flow[a];

          _res_cap[_arc_idf[a]] = cap[a] - fa;

          _res_cap[_arc_idb[a]] = fa;

        }

        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {

          int ra = _reverse[a];

          _res_cap[a] = 0;

          _res_cap[ra] = 0;

          _cost[a] = 0;

          _cost[ra] = 0;

        }

      }

      return OPTIMAL;

    }

    // Execute the algorithm and transform the results

    void start(Method method) {

      // Maximum path length for partial augment

      const int MAX_PATH_LENGTH = 4;

      // Initialize data structures for buckets

      _max_rank = _alpha * _res_node_num;

      _buckets.resize(_max_rank);

      _bucket_next.resize(_res_node_num + 1);

      _bucket_prev.resize(_res_node_num + 1);

      _rank.resize(_res_node_num + 1);

      // Execute the algorithm

      switch (method) {

        case PUSH:

          startPush();

          break;

        case AUGMENT:

          break;

        case PARTIAL_AUGMENT:

          startAugment(MAX_PATH_LENGTH);

          break;

      }

      // Compute node potentials for the original costs

      _arc_vec.clear();

      _cost_vec.clear();

      for (int j = 0; j != _res_arc_num; ++j) {

        if (_res_cap[j] > 0) {

          _arc_vec.push_back(IntPair(_source[j], _target[j]));

          _cost_vec.push_back(_scost[j]);

        }

      }

      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());

      typename BellmanFord<StaticDigraph, LargeCostArcMap>

        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);

      bf.distMap(_pi_map);

      bf.init(0);

      bf.start();

      // Handle non-zero lower bounds

      if (_have_lower) {

        int limit = _first_out[_root];

        for (int j = 0; j != limit; ++j) {

          if (!_forward[j]) _res_cap[j] += _lower[j];

        }

      }

    }

    // Initialize a cost scaling phase

    void initPhase() {

      // Saturate arcs not satisfying the optimality condition

      for (int u = 0; u != _res_node_num; ++u) {

        int last_out = _first_out[u+1];

        LargeCost pi_u = _pi[u];

        for (int a = _first_out[u]; a != last_out; ++a) {

          int v = _target[a];

          if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {

            Value delta = _res_cap[a];

            _excess[u] -= delta;

            _excess[v] += delta;

            _res_cap[a] = 0;

            _res_cap[_reverse[a]] += delta;

          }

        }

                int new_rank_v = old_rank_v;

                if (nrc < LargeCost(_max_rank))

                  new_rank_v = r + 1 + int(nrc);

                // Change the rank of v

                if (new_rank_v < old_rank_v) {

                  _rank[v] = new_rank_v;

                  _next_out[v] = _first_out[v];

                  // Remove v from its old bucket

                  if (old_rank_v < _max_rank) {

                    if (_buckets[old_rank_v] == v) {

                      _buckets[old_rank_v] = _bucket_next[v];

                    } else {

                      _bucket_next[_bucket_prev[v]] = _bucket_next[v];

                      _bucket_prev[_bucket_next[v]] = _bucket_prev[v];

                    }

                  }

                  // Insert v to its new bucket

                  _bucket_next[v] = _buckets[new_rank_v];

                  _bucket_prev[_buckets[new_rank_v]] = v;

                  _buckets[new_rank_v] = v;

                }

              }

            }

          }

          // Finish search if there are no more active nodes

          if (_excess[u] > 0) {

            total_excess -= _excess[u];

            if (total_excess <= 0) break;

          }

        }

        if (total_excess <= 0) break;

      }

      // Relabel nodes

      for (int u = 0; u != _res_node_num; ++u) {

        int k = std::min(_rank[u], r);

        if (k > 0) {

          _pi[u] -= _epsilon * k;

          _next_out[u] = _first_out[u];

        }

      }

    }

    /// Execute the algorithm performing augment and relabel operations

      // Paramters for heuristics

      const int EARLY_TERM_EPSILON_LIMIT = 1000;

      const double GLOBAL_UPDATE_FACTOR = 3.0;

      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *

        (_res_node_num + _sup_node_num * _sup_node_num));

      int next_update_limit = global_update_freq;

      int relabel_cnt = 0;

      // Perform cost scaling phases

      std::vector<int> path;

      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?

                                        1 : _epsilon / _alpha )

      {

        // Early termination heuristic

        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {

          if (earlyTermination()) break;

        }

        // Initialize current phase

        initPhase();

        // Perform partial augment and relabel operations

        while (true) {

          // Select an active node (FIFO selection)

          while (_active_nodes.size() > 0 &&

                 _excess[_active_nodes.front()] <= 0) {

            _active_nodes.pop_front();

          }

          if (_active_nodes.size() == 0) break;

          int start = _active_nodes.front();

          // Find an augmenting path from the start node

          path.clear();

          int tip = start;

          while (_excess[tip] >= 0 && int(path.size()) < max_length) {

            int u;

            LargeCost min_red_cost, rc, pi_tip = _pi[tip];

            int last_out = _first_out[tip+1];

            for (int a = _next_out[tip]; a != last_out; ++a) {

              u = _target[a];

              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {

                path.push_back(a);

                _next_out[tip] = a;

                tip = u;

                goto next_step;

              }