RhodeCode

      return *_level;

    }

    /// \brief Sets the tolerance used by the algorithm.

    ///

    /// Sets the tolerance object used by the algorithm.

    /// \return <tt>(*this)</tt>

    Preflow& tolerance(const Tolerance& tolerance) {

      _tolerance = tolerance;

      return *this;

    }

    /// \brief Returns a const reference to the tolerance.

    ///

    /// Returns a const reference to the tolerance object used by

    /// the algorithm.

    const Tolerance& tolerance() const {

      return _tolerance;

    }

    /// \name Execution Control

    /// The simplest way to execute the preflow algorithm is to use

    /// \ref run() or \ref runMinCut().\n

    /// If you need better control on the initial solution or the execution,

    /// you have to call one of the \ref init() functions first, then

    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().

    ///@{

    /// \brief Initializes the internal data structures.

    ///

    /// Initializes the internal data structures and sets the initial

    /// flow to zero on each arc.

    void init() {

      createStructures();

      _phase = true;

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

        (*_excess)[n] = 0;

      }

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

        _flow->set(e, 0);

      }

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

      _level->initStart();

      _level->initAddItem(_target);

      std::vector<Node> queue;

      reached[_source] = true;

      queue.push_back(_target);

      reached[_target] = true;

      while (!queue.empty()) {

        _level->initNewLevel();

        std::vector<Node> nqueue;

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

          Node n = queue[i];

          for (InArcIt e(_graph, n); e != INVALID; ++e) {

            Node u = _graph.source(e);

            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {

              reached[u] = true;

              _level->initAddItem(u);

              nqueue.push_back(u);

            }

          }

        }

        queue.swap(nqueue);

      }

      _level->initFinish();

      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {

        if (_tolerance.positive((*_capacity)[e])) {

          Node u = _graph.target(e);

          if ((*_level)[u] == _level->maxLevel()) continue;

          _flow->set(e, (*_capacity)[e]);

          (*_excess)[u] += (*_capacity)[e];

          if (u != _target && !_level->active(u)) {

            _level->activate(u);

          }

        }

      }

    }

    /// \brief Initializes the internal data structures using the

    /// given flow map.

    ///

    /// Initializes the internal data structures and sets the initial

    /// flow to the given \c flowMap. The \c flowMap should contain a

    /// flow or at least a preflow, i.e. at each node excluding the

    /// source node the incoming flow should greater or equal to the

    /// outgoing flow.

    /// \return \c false if the given \c flowMap is not a preflow.

    template <typename FlowMap>

    bool init(const FlowMap& flowMap) {

      createStructures();

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

        _flow->set(e, flowMap[e]);

      }

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

        Value excess = 0;

        for (InArcIt e(_graph, n); e != INVALID; ++e) {

          excess += (*_flow)[e];

        }

        for (OutArcIt e(_graph, n); e != INVALID; ++e) {

          excess -= (*_flow)[e];

        }

        if (excess < 0 && n != _source) return false;

        (*_excess)[n] = excess;

      }

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

      _level->initStart();

      _level->initAddItem(_target);

      std::vector<Node> queue;

      reached[_source] = true;

      queue.push_back(_target);

      reached[_target] = true;

      while (!queue.empty()) {

        _level->initNewLevel();

        std::vector<Node> nqueue;

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

          Node n = queue[i];

          for (InArcIt e(_graph, n); e != INVALID; ++e) {

            Node u = _graph.source(e);

            if (!reached[u] &&

                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {

              reached[u] = true;

              _level->initAddItem(u);

              nqueue.push_back(u);

            }

          }

          for (OutArcIt e(_graph, n); e != INVALID; ++e) {

            Node v = _graph.target(e);

            if (!reached[v] && _tolerance.positive((*_flow)[e])) {

              reached[v] = true;

              _level->initAddItem(v);

              nqueue.push_back(v);

            }

          }

        }

        queue.swap(nqueue);

      }

      _level->initFinish();

      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {

        Value rem = (*_capacity)[e] - (*_flow)[e];

        if (_tolerance.positive(rem)) {

          Node u = _graph.target(e);

          if ((*_level)[u] == _level->maxLevel()) continue;

          _flow->set(e, (*_capacity)[e]);

          (*_excess)[u] += rem;

          if (u != _target && !_level->active(u)) {

            _level->activate(u);

          }

        }

      }

      for (InArcIt e(_graph, _source); e != INVALID; ++e) {

        Value rem = (*_flow)[e];

        if (_tolerance.positive(rem)) {

          Node v = _graph.source(e);

          if ((*_level)[v] == _level->maxLevel()) continue;

          _flow->set(e, 0);

          (*_excess)[v] += rem;

          if (v != _target && !_level->active(v)) {

            _level->activate(v);

          }

        }

      }

      return true;

    }

    /// \brief Starts the first phase of the preflow algorithm.

    ///

    /// The preflow algorithm consists of two phases, this method runs

    /// the first phase. After the first phase the maximum flow value

    /// and a minimum value cut can already be computed, although a

    /// maximum flow is not yet obtained. So after calling this method

    /// \ref flowValue() returns the value of a maximum flow and \ref

    /// minCut() returns a minimum cut.

    /// \pre One of the \ref init() functions must be called before

    /// using this function.

    void startFirstPhase() {

      _phase = true;

        int num = _node_num;

          Value excess = (*_excess)[n];

          int new_level = _level->maxLevel();

          for (OutArcIt e(_graph, n); e != INVALID; ++e) {

            Value rem = (*_capacity)[e] - (*_flow)[e];

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.target(e);

            if ((*_level)[v] < level) {

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

              if (!_tolerance.less(rem, excess)) {

                _flow->set(e, (*_flow)[e] + excess);

                (*_excess)[v] += excess;

                excess = 0;

                goto no_more_push_1;

              } else {

                excess -= rem;

                (*_excess)[v] += rem;

                _flow->set(e, (*_capacity)[e]);

              }

            } else if (new_level > (*_level)[v]) {

              new_level = (*_level)[v];

            }

          }

          for (InArcIt e(_graph, n); e != INVALID; ++e) {

            Value rem = (*_flow)[e];

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.source(e);

            if ((*_level)[v] < level) {

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

              if (!_tolerance.less(rem, excess)) {

                _flow->set(e, (*_flow)[e] - excess);

                (*_excess)[v] += excess;

                excess = 0;

                goto no_more_push_1;

              } else {

                excess -= rem;

                (*_excess)[v] += rem;

                _flow->set(e, 0);

              }

            } else if (new_level > (*_level)[v]) {

              new_level = (*_level)[v];

            }

          }

        no_more_push_1:

          (*_excess)[n] = excess;

          if (excess != 0) {

            if (new_level + 1 < _level->maxLevel()) {

              _level->liftHighestActive(new_level + 1);

            } else {

              _level->liftHighestActiveToTop();

            }

            if (_level->emptyLevel(level)) {

              _level->liftToTop(level);

            }

          } else {

            _level->deactivate(n);

          }

        }

        num = _node_num * 20;

          Value excess = (*_excess)[n];

          int new_level = _level->maxLevel();

          for (OutArcIt e(_graph, n); e != INVALID; ++e) {

            Value rem = (*_capacity)[e] - (*_flow)[e];

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.target(e);

            if ((*_level)[v] < level) {

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

              if (!_tolerance.less(rem, excess)) {

                _flow->set(e, (*_flow)[e] + excess);

                (*_excess)[v] += excess;

                excess = 0;

                goto no_more_push_2;

              } else {

                excess -= rem;

                (*_excess)[v] += rem;

                _flow->set(e, (*_capacity)[e]);

              }

            } else if (new_level > (*_level)[v]) {

              new_level = (*_level)[v];

            }

          }

          for (InArcIt e(_graph, n); e != INVALID; ++e) {

            Value rem = (*_flow)[e];

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.source(e);

            if ((*_level)[v] < level) {

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

              if (!_tolerance.less(rem, excess)) {

                _flow->set(e, (*_flow)[e] - excess);

                (*_excess)[v] += excess;

                excess = 0;

                goto no_more_push_2;

              } else {

                excess -= rem;

                (*_excess)[v] += rem;

                _flow->set(e, 0);

              }

            } else if (new_level > (*_level)[v]) {

              new_level = (*_level)[v];

            }

          }

        no_more_push_2:

          (*_excess)[n] = excess;

          if (excess != 0) {

            if (new_level + 1 < _level->maxLevel()) {

              _level->liftActiveOn(level, new_level + 1);

            } else {

              _level->liftActiveToTop(level);

            }

            if (_level->emptyLevel(level)) {

              _level->liftToTop(level);

            }

          } else {

            _level->deactivate(n);

          }

          }

          }

    }

    /// \brief Starts the second phase of the preflow algorithm.

    ///

    /// The preflow algorithm consists of two phases, this method runs

    /// the second phase. After calling one of the \ref init() functions

    /// and \ref startFirstPhase() and then \ref startSecondPhase(),

    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the

    /// value of a maximum flow, \ref minCut() returns a minimum cut

    /// \pre One of the \ref init() functions and \ref startFirstPhase()

    /// must be called before using this function.

    void startSecondPhase() {

      _phase = false;

      typename Digraph::template NodeMap<bool> reached(_graph);

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

        reached[n] = (*_level)[n] < _level->maxLevel();

      }

      _level->initStart();

      _level->initAddItem(_source);

      std::vector<Node> queue;

      queue.push_back(_source);

      reached[_source] = true;

      while (!queue.empty()) {

        _level->initNewLevel();

        std::vector<Node> nqueue;

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

          Node n = queue[i];

          for (OutArcIt e(_graph, n); e != INVALID; ++e) {

            Node v = _graph.target(e);

            if (!reached[v] && _tolerance.positive((*_flow)[e])) {

              reached[v] = true;

              _level->initAddItem(v);

              nqueue.push_back(v);

            }

          }

          for (InArcIt e(_graph, n); e != INVALID; ++e) {

            Node u = _graph.source(e);

            if (!reached[u] &&

                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {

              reached[u] = true;

              _level->initAddItem(u);

              nqueue.push_back(u);

            }

          }

        }

        queue.swap(nqueue);

      }

      _level->initFinish();

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

        if (!reached[n]) {

          _level->dirtyTopButOne(n);

        } else if ((*_excess)[n] > 0 && _target != n) {

          _level->activate(n);

        }

      }

      Node n;

      while ((n = _level->highestActive()) != INVALID) {

        Value excess = (*_excess)[n];

        int level = _level->highestActiveLevel();

        int new_level = _level->maxLevel();

        for (OutArcIt e(_graph, n); e != INVALID; ++e) {

          Value rem = (*_capacity)[e] - (*_flow)[e];

          if (!_tolerance.positive(rem)) continue;

          Node v = _graph.target(e);

          if ((*_level)[v] < level) {

            if (!_level->active(v) && v != _source) {

              _level->activate(v);

            }

            if (!_tolerance.less(rem, excess)) {

              _flow->set(e, (*_flow)[e] + excess);

              (*_excess)[v] += excess;

              excess = 0;

              goto no_more_push;

            } else {

              excess -= rem;

              (*_excess)[v] += rem;

              _flow->set(e, (*_capacity)[e]);

            }

          } else if (new_level > (*_level)[v]) {

            new_level = (*_level)[v];

          }

        }

        for (InArcIt e(_graph, n); e != INVALID; ++e) {

          Value rem = (*_flow)[e];

          if (!_tolerance.positive(rem)) continue;

          Node v = _graph.source(e);

          if ((*_level)[v] < level) {

            if (!_level->active(v) && v != _source) {

              _level->activate(v);

            }

            if (!_tolerance.less(rem, excess)) {

              _flow->set(e, (*_flow)[e] - excess);

              (*_excess)[v] += excess;

              excess = 0;

              goto no_more_push;

            } else {

              excess -= rem;

              (*_excess)[v] += rem;

              _flow->set(e, 0);

            }

          } else if (new_level > (*_level)[v]) {

            new_level = (*_level)[v];

          }

        }

      no_more_push:

        (*_excess)[n] = excess;

        if (excess != 0) {

          if (new_level + 1 < _level->maxLevel()) {

            _level->liftHighestActive(new_level + 1);

          } else {

            // Calculation error

            _level->liftHighestActiveToTop();

          }

          if (_level->emptyLevel(level)) {

            // Calculation error

            _level->liftToTop(level);

          }

        } else {

          _level->deactivate(n);

        }

      }

    }

    /// \brief Runs the preflow algorithm.

    ///

    /// Runs the preflow algorithm.

    /// \note pf.run() is just a shortcut of the following code.

    /// \code

    ///   pf.init();

    ///   pf.startFirstPhase();

    ///   pf.startSecondPhase();

    /// \endcode

    void run() {

      init();

      startFirstPhase();

      startSecondPhase();

    }

    /// \brief Runs the preflow algorithm to compute the minimum cut.

    ///

    /// Runs the preflow algorithm to compute the minimum cut.

    /// \note pf.runMinCut() is just a shortcut of the following code.

    /// \code

    ///   pf.init();

    ///   pf.startFirstPhase();

    /// \endcode

    void runMinCut() {

      init();

      startFirstPhase();

    }

    /// @}

    /// \name Query Functions

    /// The results of the preflow algorithm can be obtained using these

    /// functions.\n

    /// Either one of the \ref run() "run*()" functions or one of the

    /// \ref startFirstPhase() "start*()" functions should be called

    /// before using them.

    ///@{

    /// \brief Returns the value of the maximum flow.

    ///

    /// Returns the value of the maximum flow by returning the excess

    /// of the target node. This value equals to the value of

    /// the maximum flow already after the first phase of the algorithm.

    ///

    /// \pre Either \ref run() or \ref init() must be called before

    /// using this function.

    Value flowValue() const {

      return (*_excess)[_target];

    }

    /// \brief Returns the flow value on the given arc.

    ///

    /// Returns the flow value on the given arc. This method can

    /// be called after the second phase of the algorithm.

    ///

    /// \pre Either \ref run() or \ref init() must be called before