RhodeCode

    typedef std::vector<char> BoolVector;

    // Note: vector<char> is used instead of vector<bool> for efficiency reasons

        Value max_sup = 0, max_dem = 0, max_cap = 0;

        for (int i = 0; i != _root; ++i) {

          int last_out = _first_out[i+1] - 1;

          for (int j = _first_out[i]; j != last_out; ++j) {

            if (_res_cap[j] > max_cap) max_cap = _res_cap[j];

          }

    typedef std::vector<char> BoolVector;

    // Note: vector<char> is used instead of vector<bool> for efficiency reasons

    int _sup_node_num;

    IntVector _buckets;

    IntVector _bucket_next;

    IntVector _bucket_prev;

    IntVector _rank;

    int _max_rank;

      _sup_node_num = 0;

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

        if (sup[n] > 0) ++_sup_node_num;

      }

          _res_cap[a] = 0;

      // 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);

    // 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;

          }

        }

      }

      // Find active nodes (i.e. nodes with positive excess)

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

        if (_excess[u] > 0) _active_nodes.push_back(u);

      }

      // Initialize the next arcs

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

        _next_out[u] = _first_out[u];

      }

    }

    // Early termination heuristic

    bool earlyTermination() {

      const double EARLY_TERM_FACTOR = 3.0;

      // Build a static residual graph

      _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(_cost[j] + 1);

        }

      }

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

      // Run Bellman-Ford algorithm to check if the current flow is optimal

      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);

      bf.init(0);

      bool done = false;

      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));

      for (int i = 0; i < K && !done; ++i) {

        done = bf.processNextWeakRound();

      }

      return done;

    }

    // Global potential update heuristic

    void globalUpdate() {

      int bucket_end = _root + 1;

      // Initialize buckets

      for (int r = 0; r != _max_rank; ++r) {

        _buckets[r] = bucket_end;

      }

      Value total_excess = 0;

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

        if (_excess[i] < 0) {

          _rank[i] = 0;

          _bucket_next[i] = _buckets[0];

          _bucket_prev[_buckets[0]] = i;

          _buckets[0] = i;

        } else {

          total_excess += _excess[i];

          _rank[i] = _max_rank;

        }

      }

      if (total_excess == 0) return;

      // Search the buckets

      int r = 0;

      for ( ; r != _max_rank; ++r) {

        while (_buckets[r] != bucket_end) {

          // Remove the first node from the current bucket

          int u = _buckets[r];

          _buckets[r] = _bucket_next[u];

          // Search the incomming arcs of u

          LargeCost pi_u = _pi[u];

          int last_out = _first_out[u+1];

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

            int ra = _reverse[a];

            if (_res_cap[ra] > 0) {

              int v = _source[ra];

              int old_rank_v = _rank[v];

              if (r < old_rank_v) {

                // Compute the new rank of v

                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;

                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];

        }

      }

    }

      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;

      std::vector<int> path;

        // Early termination heuristic

        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {

          if (earlyTermination()) break;

        // Initialize current phase

        initPhase();

          path.clear();

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

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

            int last_out = _first_out[tip+1];

              u = _target[a];

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

                path.push_back(a);

            min_red_cost = std::numeric_limits<LargeCost>::max();

            if (tip != start) {

              int ra = _reverse[path.back()];

              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];

            }

              rc = _cost[a] + pi_tip - _pi[_target[a]];

            ++relabel_cnt;

              tip = _source[path.back()];

              path.pop_back();

          int pa, u, v = start;

          for (int i = 0; i != int(path.size()); ++i) {

            pa = path[i];

            v = _target[pa];

          // Global update heuristic

          if (relabel_cnt >= next_update_limit) {

            globalUpdate();

            next_update_limit += global_update_freq;

          }

      const int EARLY_TERM_EPSILON_LIMIT = 1000;

      const double GLOBAL_UPDATE_FACTOR = 2.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;

      LargeCostVector hyper_cost(_res_node_num);

        // Early termination heuristic

        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {

          if (earlyTermination()) break;

        // Initialize current phase

        initPhase();

          LargeCost min_red_cost, rc, pi_n;

        next_node:

          last_out = _first_out[n+1];

          pi_n = _pi[n];

                  _cost[a] + pi_n - _pi[_target[a]] < 0) {

                int last_out_t = _first_out[t+1];

                LargeCost pi_t = _pi[t];

                      _cost[ta] + pi_t - _pi[_target[ta]] < 0)

                if (ahead < delta && !hyper[t]) {

                  hyper_cost[t] = _cost[a] + pi_n - pi_t;

             min_red_cost = hyper[n] ? -hyper_cost[n] :

               std::numeric_limits<LargeCost>::max();

              rc = _cost[a] + pi_n - _pi[_target[a]];

            _next_out[n] = _first_out[n];

            ++relabel_cnt;

          // Global update heuristic

          if (relabel_cnt >= next_update_limit) {

            globalUpdate();

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

              hyper[u] = false;

            next_update_limit += global_update_freq;

          }

    typedef std::vector<char> BoolVector;

    // Note: vector<char> is used instead of vector<bool> for efficiency reasons

    BoolVector _forward;

      BoolVector reached(_res_node_num);

      BoolVector processed(_res_node_num);

    typedef std::vector<char> BoolVector;

    // Note: vector<char> is used instead of vector<bool> for efficiency reasons

    BoolVector _forward;

    BoolVector _state;

      const BoolVector &_state;

        for (int e = _next_arc; e != _search_arc_num; ++e) {

        for (int e = 0; e != _next_arc; ++e) {

      const BoolVector &_state;

        for (int e = 0; e != _search_arc_num; ++e) {

      const BoolVector &_state;

        const double BLOCK_SIZE_FACTOR = 1.0;

        for (e = _next_arc; e != _search_arc_num; ++e) {

        for (e = 0; e != _next_arc; ++e) {

      const BoolVector &_state;

        for (e = _next_arc; e != _search_arc_num; ++e) {

        for (e = 0; e != _next_arc; ++e) {

      const BoolVector &_state;

        for (int i = 0; i != _curr_length; ++i) {

        for (e = _next_arc; e != _search_arc_num; ++e) {

        for (e = 0; e != _next_arc; ++e) {

      for (int i = 0; i != int(_dirty_revs.size()); ++i) {

    // Heuristic initial pivots

    bool initialPivots() {

      Value curr, total = 0;

      std::vector<Node> supply_nodes, demand_nodes;

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

        curr = _supply[_node_id[u]];

        if (curr > 0) {

          total += curr;

          supply_nodes.push_back(u);

        }

        else if (curr < 0) {

          demand_nodes.push_back(u);

        }

      }

      if (_sum_supply > 0) total -= _sum_supply;

      if (total <= 0) return true;

      IntVector arc_vector;

      if (_sum_supply >= 0) {

        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {

          // Perform a reverse graph search from the sink to the source

          typename GR::template NodeMap<bool> reached(_graph, false);

          Node s = supply_nodes[0], t = demand_nodes[0];

          std::vector<Node> stack;

          reached[t] = true;

          stack.push_back(t);

          while (!stack.empty()) {

            Node u, v = stack.back();

            stack.pop_back();

            if (v == s) break;

            for (InArcIt a(_graph, v); a != INVALID; ++a) {

              if (reached[u = _graph.source(a)]) continue;

              int j = _arc_id[a];

              if (_cap[j] >= total) {

                arc_vector.push_back(j);

                reached[u] = true;

                stack.push_back(u);

              }

            }

          }

        } else {

          // Find the min. cost incomming arc for each demand node

          for (int i = 0; i != int(demand_nodes.size()); ++i) {

            Node v = demand_nodes[i];

            Cost c, min_cost = std::numeric_limits<Cost>::max();

            Arc min_arc = INVALID;

            for (InArcIt a(_graph, v); a != INVALID; ++a) {

              c = _cost[_arc_id[a]];

              if (c < min_cost) {

                min_cost = c;

                min_arc = a;

              }

            }

            if (min_arc != INVALID) {

              arc_vector.push_back(_arc_id[min_arc]);

            }

          }

        }

      } else {

        // Find the min. cost outgoing arc for each supply node

        for (int i = 0; i != int(supply_nodes.size()); ++i) {

          Node u = supply_nodes[i];

          Cost c, min_cost = std::numeric_limits<Cost>::max();

          Arc min_arc = INVALID;

          for (OutArcIt a(_graph, u); a != INVALID; ++a) {

            c = _cost[_arc_id[a]];

            if (c < min_cost) {

              min_cost = c;

              min_arc = a;

            }

          }

          if (min_arc != INVALID) {

            arc_vector.push_back(_arc_id[min_arc]);

          }

        }

      }

      // Perform heuristic initial pivots

      for (int i = 0; i != int(arc_vector.size()); ++i) {

        in_arc = arc_vector[i];

        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -

            _pi[_target[in_arc]]) >= 0) continue;

        findJoinNode();

        bool change = findLeavingArc();

        if (delta >= MAX) return false;

        changeFlow(change);

        if (change) {

          updateTreeStructure();

          updatePotential();

        }

      }

      return true;

    }

      // Perform heuristic initial pivots

      if (!initialPivots()) return UNBOUNDED;