RhodeCode

    // Process all rounds of computing path data for the current component.

    // _data[v][k] is the length of a shortest directed walk from the root

    // node to node v containing exactly k arcs.

    void processRounds() {

      Node start = (*_nodes)[0];

      _data[start][0] = PathData(0);

      _process.clear();

      _process.push_back(start);

      int k, n = _nodes->size();

      int next_check = 4;

      bool terminate = false;

      for (k = 1; k <= n && int(_process.size()) < n && !terminate; ++k) {

        processNextBuildRound(k);

        if (k == next_check || k == n) {

          terminate = checkTermination(k);

          next_check = next_check * 3 / 2;

        }

      }

      for ( ; k <= n && !terminate; ++k) {

        processNextFullRound(k);

        if (k == next_check || k == n) {

          terminate = checkTermination(k);

          next_check = next_check * 3 / 2;

        }

      }

    }

    // Process one round and rebuild _process

    void processNextBuildRound(int k) {

      std::vector<Node> next;

      Node u, v;

      Arc e;

      LargeValue d;

      for (int i = 0; i < int(_process.size()); ++i) {

        u = _process[i];

        for (int j = 0; j < int(_out_arcs[u].size()); ++j) {

          e = _out_arcs[u][j];

          v = _gr.target(e);

          d = _data[u][k-1].dist + _length[e];

          if (_tolerance.less(d, _data[v][k].dist)) {

            if (_data[v][k].dist == INF) next.push_back(v);

            _data[v][k] = PathData(d, e);

          }

        }

      }

      _process.swap(next);

    }

    // Process one round using _nodes instead of _process

    void processNextFullRound(int k) {

      Node u, v;

      Arc e;

      LargeValue d;

      for (int i = 0; i < int(_nodes->size()); ++i) {

        u = (*_nodes)[i];

        for (int j = 0; j < int(_out_arcs[u].size()); ++j) {

          e = _out_arcs[u][j];

          v = _gr.target(e);

          d = _data[u][k-1].dist + _length[e];

          if (_tolerance.less(d, _data[v][k].dist)) {

            _data[v][k] = PathData(d, e);

          }

        }

      }

    }

    // Check early termination

    bool checkTermination(int k) {

      typedef std::pair<int, int> Pair;

      typename GR::template NodeMap<Pair> level(_gr, Pair(-1, 0));

      typename GR::template NodeMap<LargeValue> pi(_gr);

      int n = _nodes->size();

      LargeValue length;

      int size;

      Node u;

      // Search for cycles that are already found

      _curr_found = false;

      for (int i = 0; i < n; ++i) {

        u = (*_nodes)[i];

        if (_data[u][k].dist == INF) continue;

        for (int j = k; j >= 0; --j) {

          if (level[u].first == i && level[u].second > 0) {

            // A cycle is found

            length = _data[u][level[u].second].dist - _data[u][j].dist;

            size = level[u].second - j;

            if (!_curr_found || length * _curr_size < _curr_length * size) {

              _curr_length = length;

              _curr_size = size;

              _curr_node = u;

              _curr_level = level[u].second;

              _curr_found = true;

            }

          }

          level[u] = Pair(i, j);

        }

      }

      // If at least one cycle is found, check the optimality condition

      LargeValue d;

      if (_curr_found && k < n) {

        // Find node potentials

        for (int i = 0; i < n; ++i) {

          u = (*_nodes)[i];

          pi[u] = INF;

          for (int j = 0; j <= k; ++j) {

            if (_data[u][j].dist < INF) {

              d = _data[u][j].dist * _curr_size - j * _curr_length;

              if (_tolerance.less(d, pi[u])) pi[u] = d;

            }

          }

        }

        // Check the optimality condition for all arcs

        bool done = true;

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

          if (_tolerance.less(_length[a] * _curr_size - _curr_length,

                              pi[_gr.target(a)] - pi[_gr.source(a)]) ) {

            done = false;

            break;

          }

        }

        return done;

      }

      return (k == n);

    }

  }; //class HartmannOrlin

  ///@}

} //namespace lemon

#endif //LEMON_HARTMANN_ORLIN_H