1.1 --- a/lemon/preflow.h Wed Nov 14 17:53:08 2007 +0000
1.2 +++ b/lemon/preflow.h Sat Nov 17 20:58:11 2007 +0000
1.3 @@ -19,897 +19,897 @@
1.4 #ifndef LEMON_PREFLOW_H
1.5 #define LEMON_PREFLOW_H
1.6
1.7 -#include <vector>
1.8 -#include <queue>
1.9 -
1.10 #include <lemon/error.h>
1.11 #include <lemon/bits/invalid.h>
1.12 #include <lemon/tolerance.h>
1.13 #include <lemon/maps.h>
1.14 #include <lemon/graph_utils.h>
1.15 +#include <lemon/elevator.h>
1.16
1.17 /// \file
1.18 /// \ingroup max_flow
1.19 /// \brief Implementation of the preflow algorithm.
1.20
1.21 -namespace lemon {
1.22 +namespace lemon {
1.23 +
1.24 + /// \brief Default traits class of Preflow class.
1.25 + ///
1.26 + /// Default traits class of Preflow class.
1.27 + /// \param _Graph Graph type.
1.28 + /// \param _CapacityMap Type of capacity map.
1.29 + template <typename _Graph, typename _CapacityMap>
1.30 + struct PreflowDefaultTraits {
1.31
1.32 - ///\ingroup max_flow
1.33 - ///\brief %Preflow algorithms class.
1.34 + /// \brief The graph type the algorithm runs on.
1.35 + typedef _Graph Graph;
1.36 +
1.37 + /// \brief The type of the map that stores the edge capacities.
1.38 + ///
1.39 + /// The type of the map that stores the edge capacities.
1.40 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.41 + typedef _CapacityMap CapacityMap;
1.42 +
1.43 + /// \brief The type of the length of the edges.
1.44 + typedef typename CapacityMap::Value Value;
1.45 +
1.46 + /// \brief The map type that stores the flow values.
1.47 + ///
1.48 + /// The map type that stores the flow values.
1.49 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1.50 + typedef typename Graph::template EdgeMap<Value> FlowMap;
1.51 +
1.52 + /// \brief Instantiates a FlowMap.
1.53 + ///
1.54 + /// This function instantiates a \ref FlowMap.
1.55 + /// \param graph The graph, to which we would like to define the flow map.
1.56 + static FlowMap* createFlowMap(const Graph& graph) {
1.57 + return new FlowMap(graph);
1.58 + }
1.59 +
1.60 + /// \brief The eleavator type used by Preflow algorithm.
1.61 + ///
1.62 + /// The elevator type used by Preflow algorithm.
1.63 + ///
1.64 + /// \sa Elevator
1.65 + /// \sa LinkedElevator
1.66 + typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
1.67 +
1.68 + /// \brief Instantiates an Elevator.
1.69 + ///
1.70 + /// This function instantiates a \ref Elevator.
1.71 + /// \param graph The graph, to which we would like to define the elevator.
1.72 + /// \param max_level The maximum level of the elevator.
1.73 + static Elevator* createElevator(const Graph& graph, int max_level) {
1.74 + return new Elevator(graph, max_level);
1.75 + }
1.76 +
1.77 + /// \brief The tolerance used by the algorithm
1.78 + ///
1.79 + /// The tolerance used by the algorithm to handle inexact computation.
1.80 + typedef Tolerance<Value> Tolerance;
1.81 +
1.82 + };
1.83 +
1.84 +
1.85 + /// \ingroup max_flow
1.86 ///
1.87 - ///This class provides an implementation of the \e preflow \e
1.88 - ///algorithm producing a flow of maximum value in a directed
1.89 - ///graph. The preflow algorithms are the fastest known max flow algorithms.
1.90 - ///The \e source node, the \e target node, the \e
1.91 - ///capacity of the edges and the \e starting \e flow value of the
1.92 - ///edges should be passed to the algorithm through the
1.93 - ///constructor. It is possible to change these quantities using the
1.94 - ///functions \ref source, \ref target, \ref capacityMap and \ref
1.95 - ///flowMap.
1.96 + /// \brief %Preflow algorithms class.
1.97 ///
1.98 - ///After running \ref lemon::Preflow::phase1() "phase1()"
1.99 - ///or \ref lemon::Preflow::run() "run()", the maximal flow
1.100 - ///value can be obtained by calling \ref flowValue(). The minimum
1.101 - ///value cut can be written into a <tt>bool</tt> node map by
1.102 - ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
1.103 - ///the inclusionwise minimum and maximum of the minimum value cuts,
1.104 - ///resp.)
1.105 + /// This class provides an implementation of the Goldberg's \e
1.106 + /// preflow \e algorithm producing a flow of maximum value in a
1.107 + /// directed graph. The preflow algorithms are the fastest known max
1.108 + /// flow algorithms. The current implementation use a mixture of the
1.109 + /// \e "highest label" and the \e "bound decrease" heuristics.
1.110 + /// The worst case time complexity of the algorithm is \f$O(n^3)\f$.
1.111 ///
1.112 - ///\param Graph The directed graph type the algorithm runs on.
1.113 - ///\param Num The number type of the capacities and the flow values.
1.114 - ///\param CapacityMap The capacity map type.
1.115 - ///\param FlowMap The flow map type.
1.116 - ///\param Tol The tolerance type.
1.117 + /// The algorithm consists from two phases. After the first phase
1.118 + /// the maximal flow value and the minimum cut can be obtained. The
1.119 + /// second phase constructs the feasible maximum flow on each edge.
1.120 ///
1.121 - ///\author Jacint Szabo
1.122 - ///\todo Second template parameter is superfluous
1.123 - template <typename Graph, typename Num,
1.124 - typename CapacityMap=typename Graph::template EdgeMap<Num>,
1.125 - typename FlowMap=typename Graph::template EdgeMap<Num>,
1.126 - typename Tol=Tolerance<Num> >
1.127 + /// \param _Graph The directed graph type the algorithm runs on.
1.128 + /// \param _CapacityMap The flow map type.
1.129 + /// \param _Traits Traits class to set various data types used by
1.130 + /// the algorithm. The default traits class is \ref
1.131 + /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the
1.132 + /// documentation of a %Preflow traits class.
1.133 + ///
1.134 + ///\author Jacint Szabo and Balazs Dezso
1.135 +#ifdef DOXYGEN
1.136 + template <typename _Graph, typename _CapacityMap, typename _Traits>
1.137 +#else
1.138 + template <typename _Graph,
1.139 + typename _CapacityMap = typename _Graph::template EdgeMap<int>,
1.140 + typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
1.141 +#endif
1.142 class Preflow {
1.143 - protected:
1.144 - typedef typename Graph::Node Node;
1.145 - typedef typename Graph::NodeIt NodeIt;
1.146 - typedef typename Graph::EdgeIt EdgeIt;
1.147 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.148 - typedef typename Graph::InEdgeIt InEdgeIt;
1.149 + public:
1.150 +
1.151 + typedef _Traits Traits;
1.152 + typedef typename Traits::Graph Graph;
1.153 + typedef typename Traits::CapacityMap CapacityMap;
1.154 + typedef typename Traits::Value Value;
1.155
1.156 - typedef typename Graph::template NodeMap<Node> NNMap;
1.157 - typedef typename std::vector<Node> VecNode;
1.158 + typedef typename Traits::FlowMap FlowMap;
1.159 + typedef typename Traits::Elevator Elevator;
1.160 + typedef typename Traits::Tolerance Tolerance;
1.161
1.162 - const Graph* _g;
1.163 - Node _source;
1.164 - Node _target;
1.165 + /// \brief \ref Exception for uninitialized parameters.
1.166 + ///
1.167 + /// This error represents problems in the initialization
1.168 + /// of the parameters of the algorithms.
1.169 + class UninitializedParameter : public lemon::UninitializedParameter {
1.170 + public:
1.171 + virtual const char* what() const throw() {
1.172 + return "lemon::Preflow::UninitializedParameter";
1.173 + }
1.174 + };
1.175 +
1.176 + private:
1.177 +
1.178 + GRAPH_TYPEDEFS(typename Graph);
1.179 +
1.180 + const Graph& _graph;
1.181 const CapacityMap* _capacity;
1.182 +
1.183 + int _node_num;
1.184 +
1.185 + Node _source, _target;
1.186 +
1.187 FlowMap* _flow;
1.188 + bool _local_flow;
1.189
1.190 - Tol _surely;
1.191 -
1.192 - int _node_num; //the number of nodes of G
1.193 -
1.194 - typename Graph::template NodeMap<int> level;
1.195 - typename Graph::template NodeMap<Num> excess;
1.196 + Elevator* _level;
1.197 + bool _local_level;
1.198
1.199 - // constants used for heuristics
1.200 - static const int H0=20;
1.201 - static const int H1=1;
1.202 + typedef typename Graph::template NodeMap<Value> ExcessMap;
1.203 + ExcessMap* _excess;
1.204 +
1.205 + Tolerance _tolerance;
1.206 +
1.207 + bool _phase;
1.208 +
1.209 + void createStructures() {
1.210 + _node_num = countNodes(_graph);
1.211 +
1.212 + if (!_flow) {
1.213 + _flow = Traits::createFlowMap(_graph);
1.214 + _local_flow = true;
1.215 + }
1.216 + if (!_level) {
1.217 + _level = Traits::createElevator(_graph, _node_num);
1.218 + _local_level = true;
1.219 + }
1.220 + if (!_excess) {
1.221 + _excess = new ExcessMap(_graph);
1.222 + }
1.223 + }
1.224 +
1.225 + void destroyStructures() {
1.226 + if (_local_flow) {
1.227 + delete _flow;
1.228 + }
1.229 + if (_local_level) {
1.230 + delete _level;
1.231 + }
1.232 + if (_excess) {
1.233 + delete _excess;
1.234 + }
1.235 + }
1.236
1.237 public:
1.238
1.239 - ///\ref Exception for the case when s=t.
1.240 + typedef Preflow Create;
1.241
1.242 - ///\ref Exception for the case when the source equals the target.
1.243 - class InvalidArgument : public lemon::LogicError {
1.244 - public:
1.245 - virtual const char* what() const throw() {
1.246 - return "lemon::Preflow::InvalidArgument";
1.247 + ///\name Named template parameters
1.248 +
1.249 + ///@{
1.250 +
1.251 + template <typename _FlowMap>
1.252 + struct DefFlowMapTraits : public Traits {
1.253 + typedef _FlowMap FlowMap;
1.254 + static FlowMap *createFlowMap(const Graph&) {
1.255 + throw UninitializedParameter();
1.256 }
1.257 };
1.258 -
1.259 -
1.260 - ///Indicates the property of the starting flow map.
1.261 -
1.262 - ///Indicates the property of the starting flow map.
1.263 +
1.264 + /// \brief \ref named-templ-param "Named parameter" for setting
1.265 + /// FlowMap type
1.266 ///
1.267 - enum FlowEnum{
1.268 - ///indicates an unspecified edge map. \c flow will be
1.269 - ///set to the constant zero flow in the beginning of
1.270 - ///the algorithm in this case.
1.271 - NO_FLOW,
1.272 - ///constant zero flow
1.273 - ZERO_FLOW,
1.274 - ///any flow, i.e. the sum of the in-flows equals to
1.275 - ///the sum of the out-flows in every node except the \c source and
1.276 - ///the \c target.
1.277 - GEN_FLOW,
1.278 - ///any preflow, i.e. the sum of the in-flows is at
1.279 - ///least the sum of the out-flows in every node except the \c source.
1.280 - PRE_FLOW
1.281 + /// \ref named-templ-param "Named parameter" for setting FlowMap
1.282 + /// type
1.283 + template <typename _FlowMap>
1.284 + struct DefFlowMap
1.285 + : public Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
1.286 + typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create;
1.287 };
1.288
1.289 - ///Indicates the state of the preflow algorithm.
1.290 + template <typename _Elevator>
1.291 + struct DefElevatorTraits : public Traits {
1.292 + typedef _Elevator Elevator;
1.293 + static Elevator *createElevator(const Graph&, int) {
1.294 + throw UninitializedParameter();
1.295 + }
1.296 + };
1.297
1.298 - ///Indicates the state of the preflow algorithm.
1.299 + /// \brief \ref named-templ-param "Named parameter" for setting
1.300 + /// Elevator type
1.301 ///
1.302 - enum StatusEnum {
1.303 - ///before running the algorithm or
1.304 - ///at an unspecified state.
1.305 - AFTER_NOTHING,
1.306 - ///right after running \ref phase1()
1.307 - AFTER_PREFLOW_PHASE_1,
1.308 - ///after running \ref phase2()
1.309 - AFTER_PREFLOW_PHASE_2
1.310 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.311 + /// type
1.312 + template <typename _Elevator>
1.313 + struct DefElevator
1.314 + : public Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > {
1.315 + typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create;
1.316 };
1.317 -
1.318 - protected:
1.319 - FlowEnum flow_prop;
1.320 - StatusEnum status; // Do not needle this flag only if necessary.
1.321 -
1.322 - public:
1.323 - ///The constructor of the class.
1.324
1.325 - ///The constructor of the class.
1.326 - ///\param _gr The directed graph the algorithm runs on.
1.327 - ///\param _s The source node.
1.328 - ///\param _t The target node.
1.329 - ///\param _cap The capacity of the edges.
1.330 - ///\param _f The flow of the edges.
1.331 - ///\param _sr Tol class.
1.332 - ///Except the graph, all of these parameters can be reset by
1.333 - ///calling \ref source, \ref target, \ref capacityMap and \ref
1.334 - ///flowMap, resp.
1.335 - Preflow(const Graph& _gr, Node _s, Node _t,
1.336 - const CapacityMap& _cap, FlowMap& _f,
1.337 - const Tol &_sr=Tol()) :
1.338 - _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
1.339 - _flow(&_f), _surely(_sr),
1.340 - _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
1.341 - flow_prop(NO_FLOW), status(AFTER_NOTHING) {
1.342 - if ( _source==_target )
1.343 - throw InvalidArgument();
1.344 - }
1.345 -
1.346 - ///Give a reference to the tolerance handler class
1.347 + template <typename _Elevator>
1.348 + struct DefStandardElevatorTraits : public Traits {
1.349 + typedef _Elevator Elevator;
1.350 + static Elevator *createElevator(const Graph& graph, int max_level) {
1.351 + return new Elevator(graph, max_level);
1.352 + }
1.353 + };
1.354
1.355 - ///Give a reference to the tolerance handler class
1.356 - ///\sa Tolerance
1.357 - Tol &tolerance() { return _surely; }
1.358 + /// \brief \ref named-templ-param "Named parameter" for setting
1.359 + /// Elevator type
1.360 + ///
1.361 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.362 + /// type. The Elevator should be standard constructor interface, ie.
1.363 + /// the graph and the maximum level should be passed to it.
1.364 + template <typename _Elevator>
1.365 + struct DefStandardElevator
1.366 + : public Preflow<Graph, CapacityMap,
1.367 + DefStandardElevatorTraits<_Elevator> > {
1.368 + typedef Preflow<Graph, CapacityMap,
1.369 + DefStandardElevatorTraits<_Elevator> > Create;
1.370 + };
1.371
1.372 - ///Runs the preflow algorithm.
1.373 + /// @}
1.374
1.375 - ///Runs the preflow algorithm.
1.376 + /// \brief The constructor of the class.
1.377 ///
1.378 - void run() {
1.379 - phase1(flow_prop);
1.380 - phase2();
1.381 - }
1.382 -
1.383 - ///Runs the preflow algorithm.
1.384 -
1.385 - ///Runs the preflow algorithm.
1.386 - ///\pre The starting flow map must be
1.387 - /// - a constant zero flow if \c fp is \c ZERO_FLOW,
1.388 - /// - an arbitrary flow if \c fp is \c GEN_FLOW,
1.389 - /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
1.390 - /// - any map if \c fp is NO_FLOW.
1.391 - ///If the starting flow map is a flow or a preflow then
1.392 - ///the algorithm terminates faster.
1.393 - void run(FlowEnum fp) {
1.394 - flow_prop=fp;
1.395 - run();
1.396 - }
1.397 -
1.398 - ///Runs the first phase of the preflow algorithm.
1.399 + /// The constructor of the class.
1.400 + /// \param graph The directed graph the algorithm runs on.
1.401 + /// \param capacity The capacity of the edges.
1.402 + /// \param source The source node.
1.403 + /// \param target The target node.
1.404 + Preflow(const Graph& graph, const CapacityMap& capacity,
1.405 + Node source, Node target)
1.406 + : _graph(graph), _capacity(&capacity),
1.407 + _node_num(0), _source(source), _target(target),
1.408 + _flow(0), _local_flow(false),
1.409 + _level(0), _local_level(false),
1.410 + _excess(0), _tolerance(), _phase() {}
1.411
1.412 - ///The preflow algorithm consists of two phases, this method runs
1.413 - ///the first phase. After the first phase the maximum flow value
1.414 - ///and a minimum value cut can already be computed, although a
1.415 - ///maximum flow is not yet obtained. So after calling this method
1.416 - ///\ref flowValue returns the value of a maximum flow and \ref
1.417 - ///minCut returns a minimum cut.
1.418 - ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
1.419 - ///value cuts unless calling \ref phase2.
1.420 - ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
1.421 - ///is needed for this phase.
1.422 - ///\pre The starting flow must be
1.423 - ///- a constant zero flow if \c fp is \c ZERO_FLOW,
1.424 - ///- an arbitary flow if \c fp is \c GEN_FLOW,
1.425 - ///- an arbitary preflow if \c fp is \c PRE_FLOW,
1.426 - ///- any map if \c fp is NO_FLOW.
1.427 - void phase1(FlowEnum fp)
1.428 - {
1.429 - flow_prop=fp;
1.430 - phase1();
1.431 + /// \brief Destrcutor.
1.432 + ///
1.433 + /// Destructor.
1.434 + ~Preflow() {
1.435 + destroyStructures();
1.436 }
1.437
1.438 + /// \brief Sets the capacity map.
1.439 + ///
1.440 + /// Sets the capacity map.
1.441 + /// \return \c (*this)
1.442 + Preflow& capacityMap(const CapacityMap& map) {
1.443 + _capacity = ↦
1.444 + return *this;
1.445 + }
1.446 +
1.447 + /// \brief Sets the flow map.
1.448 + ///
1.449 + /// Sets the flow map.
1.450 + /// \return \c (*this)
1.451 + Preflow& flowMap(FlowMap& map) {
1.452 + if (_local_flow) {
1.453 + delete _flow;
1.454 + _local_flow = false;
1.455 + }
1.456 + _flow = ↦
1.457 + return *this;
1.458 + }
1.459 +
1.460 + /// \brief Returns the flow map.
1.461 + ///
1.462 + /// \return The flow map.
1.463 + const FlowMap& flowMap() {
1.464 + return *_flow;
1.465 + }
1.466 +
1.467 + /// \brief Sets the elevator.
1.468 + ///
1.469 + /// Sets the elevator.
1.470 + /// \return \c (*this)
1.471 + Preflow& elevator(Elevator& elevator) {
1.472 + if (_local_level) {
1.473 + delete _level;
1.474 + _local_level = false;
1.475 + }
1.476 + _level = &elevator;
1.477 + return *this;
1.478 + }
1.479 +
1.480 + /// \brief Returns the elevator.
1.481 + ///
1.482 + /// \return The elevator.
1.483 + const Elevator& elevator() {
1.484 + return *_level;
1.485 + }
1.486 +
1.487 + /// \brief Sets the source node.
1.488 + ///
1.489 + /// Sets the source node.
1.490 + /// \return \c (*this)
1.491 + Preflow& source(const Node& node) {
1.492 + _source = node;
1.493 + return *this;
1.494 + }
1.495 +
1.496 + /// \brief Sets the target node.
1.497 + ///
1.498 + /// Sets the target node.
1.499 + /// \return \c (*this)
1.500 + Preflow& target(const Node& node) {
1.501 + _target = node;
1.502 + return *this;
1.503 + }
1.504 +
1.505 + /// \brief Sets the tolerance used by algorithm.
1.506 + ///
1.507 + /// Sets the tolerance used by algorithm.
1.508 + Preflow& tolerance(const Tolerance& tolerance) const {
1.509 + _tolerance = tolerance;
1.510 + return *this;
1.511 + }
1.512 +
1.513 + /// \brief Returns the tolerance used by algorithm.
1.514 + ///
1.515 + /// Returns the tolerance used by algorithm.
1.516 + const Tolerance& tolerance() const {
1.517 + return tolerance;
1.518 + }
1.519 +
1.520 + /// \name Execution control The simplest way to execute the
1.521 + /// algorithm is to use one of the member functions called \c
1.522 + /// run().
1.523 + /// \n
1.524 + /// If you need more control on initial solution or
1.525 + /// execution then you have to call one \ref init() function and then
1.526 + /// the startFirstPhase() and if you need the startSecondPhase().
1.527
1.528 - ///Runs the first phase of the preflow algorithm.
1.529 + ///@{
1.530
1.531 - ///The preflow algorithm consists of two phases, this method runs
1.532 - ///the first phase. After the first phase the maximum flow value
1.533 - ///and a minimum value cut can already be computed, although a
1.534 - ///maximum flow is not yet obtained. So after calling this method
1.535 - ///\ref flowValue returns the value of a maximum flow and \ref
1.536 - ///minCut returns a minimum cut.
1.537 - ///\warning \ref minMinCut() and \ref maxMinCut() do not
1.538 - ///give minimum value cuts unless calling \ref phase2().
1.539 - ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
1.540 - ///is needed for this phase.
1.541 - void phase1()
1.542 - {
1.543 - int heur0=int(H0*_node_num); //time while running 'bound decrease'
1.544 - int heur1=int(H1*_node_num); //time while running 'highest label'
1.545 - int heur=heur1; //starting time interval (#of relabels)
1.546 - int numrelabel=0;
1.547 + /// \brief Initializes the internal data structures.
1.548 + ///
1.549 + /// Initializes the internal data structures.
1.550 + ///
1.551 + void init() {
1.552 + createStructures();
1.553
1.554 - bool what_heur=1;
1.555 - //It is 0 in case 'bound decrease' and 1 in case 'highest label'
1.556 + _phase = true;
1.557 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.558 + _excess->set(n, 0);
1.559 + }
1.560
1.561 - bool end=false;
1.562 - //Needed for 'bound decrease', true means no active
1.563 - //nodes are above bound b.
1.564 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.565 + _flow->set(e, 0);
1.566 + }
1.567
1.568 - int k=_node_num-2; //bound on the highest level under n containing a node
1.569 - int b=k; //bound on the highest level under n containing an active node
1.570 + typename Graph::template NodeMap<bool> reached(_graph, false);
1.571
1.572 - VecNode first(_node_num, INVALID);
1.573 - NNMap next(*_g, INVALID);
1.574 + _level->initStart();
1.575 + _level->initAddItem(_target);
1.576
1.577 - NNMap left(*_g, INVALID);
1.578 - NNMap right(*_g, INVALID);
1.579 - VecNode level_list(_node_num,INVALID);
1.580 - //List of the nodes in level i<n, set to n.
1.581 + std::vector<Node> queue;
1.582 + reached.set(_source, true);
1.583
1.584 - preflowPreproc(first, next, level_list, left, right);
1.585 -
1.586 - //Push/relabel on the highest level active nodes.
1.587 - while ( true ) {
1.588 - if ( b == 0 ) {
1.589 - if ( !what_heur && !end && k > 0 ) {
1.590 - b=k;
1.591 - end=true;
1.592 - } else break;
1.593 - }
1.594 -
1.595 - if ( first[b]==INVALID ) --b;
1.596 - else {
1.597 - end=false;
1.598 - Node w=first[b];
1.599 - first[b]=next[w];
1.600 - int newlevel=push(w, next, first);
1.601 - if ( excess[w] != 0 ) {
1.602 - relabel(w, newlevel, first, next, level_list,
1.603 - left, right, b, k, what_heur);
1.604 - }
1.605 -
1.606 - ++numrelabel;
1.607 - if ( numrelabel >= heur ) {
1.608 - numrelabel=0;
1.609 - if ( what_heur ) {
1.610 - what_heur=0;
1.611 - heur=heur0;
1.612 - end=false;
1.613 - } else {
1.614 - what_heur=1;
1.615 - heur=heur1;
1.616 - b=k;
1.617 + queue.push_back(_target);
1.618 + reached.set(_target, true);
1.619 + while (!queue.empty()) {
1.620 + std::vector<Node> nqueue;
1.621 + for (int i = 0; i < int(queue.size()); ++i) {
1.622 + Node n = queue[i];
1.623 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.624 + Node u = _graph.source(e);
1.625 + if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
1.626 + reached.set(u, true);
1.627 + _level->initAddItem(u);
1.628 + nqueue.push_back(u);
1.629 }
1.630 }
1.631 }
1.632 + queue.swap(nqueue);
1.633 }
1.634 - flow_prop=PRE_FLOW;
1.635 - status=AFTER_PREFLOW_PHASE_1;
1.636 - }
1.637 - // Heuristics:
1.638 - // 2 phase
1.639 - // gap
1.640 - // list 'level_list' on the nodes on level i implemented by hand
1.641 - // stack 'active' on the active nodes on level i
1.642 - // runs heuristic 'highest label' for H1*n relabels
1.643 - // runs heuristic 'bound decrease' for H0*n relabels,
1.644 - // starts with 'highest label'
1.645 - // Parameters H0 and H1 are initialized to 20 and 1.
1.646 + _level->initFinish();
1.647
1.648 -
1.649 - ///Runs the second phase of the preflow algorithm.
1.650 -
1.651 - ///The preflow algorithm consists of two phases, this method runs
1.652 - ///the second phase. After calling \ref phase1() and then
1.653 - ///\ref phase2(),
1.654 - /// \ref flowMap() return a maximum flow, \ref flowValue
1.655 - ///returns the value of a maximum flow, \ref minCut returns a
1.656 - ///minimum cut, while the methods \ref minMinCut and \ref
1.657 - ///maxMinCut return the inclusionwise minimum and maximum cuts of
1.658 - ///minimum value, resp. \pre \ref phase1 must be called before.
1.659 - ///
1.660 - /// \todo The inexact computation can cause positive excess on a set of
1.661 - /// unpushable nodes. We may have to watch the empty level in this case
1.662 - /// due to avoid the terrible long running time.
1.663 - void phase2()
1.664 - {
1.665 -
1.666 - int k=_node_num-2; //bound on the highest level under n containing a node
1.667 - int b=k; //bound on the highest level under n of an active node
1.668 -
1.669 -
1.670 - VecNode first(_node_num, INVALID);
1.671 - NNMap next(*_g, INVALID);
1.672 - level.set(_source,0);
1.673 - std::queue<Node> bfs_queue;
1.674 - bfs_queue.push(_source);
1.675 -
1.676 - while ( !bfs_queue.empty() ) {
1.677 -
1.678 - Node v=bfs_queue.front();
1.679 - bfs_queue.pop();
1.680 - int l=level[v]+1;
1.681 -
1.682 - for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
1.683 - if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
1.684 - Node u=_g->source(e);
1.685 - if ( level[u] >= _node_num ) {
1.686 - bfs_queue.push(u);
1.687 - level.set(u, l);
1.688 - if ( excess[u] != 0 ) {
1.689 - next.set(u,first[l]);
1.690 - first[l]=u;
1.691 - }
1.692 - }
1.693 - }
1.694 -
1.695 - for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
1.696 - if ( !_surely.positive((*_flow)[e]) ) continue;
1.697 - Node u=_g->target(e);
1.698 - if ( level[u] >= _node_num ) {
1.699 - bfs_queue.push(u);
1.700 - level.set(u, l);
1.701 - if ( excess[u] != 0 ) {
1.702 - next.set(u,first[l]);
1.703 - first[l]=u;
1.704 - }
1.705 - }
1.706 - }
1.707 - }
1.708 - b=_node_num-2;
1.709 -
1.710 - while ( true ) {
1.711 -
1.712 - if ( b == 0 ) break;
1.713 - if ( first[b]==INVALID ) --b;
1.714 - else {
1.715 - Node w=first[b];
1.716 - first[b]=next[w];
1.717 - int newlevel=push(w,next, first);
1.718 -
1.719 - if ( newlevel == _node_num) {
1.720 - excess.set(w, 0);
1.721 - level.set(w,_node_num);
1.722 - }
1.723 - //relabel
1.724 - if ( excess[w] != 0 ) {
1.725 - level.set(w,++newlevel);
1.726 - next.set(w,first[newlevel]);
1.727 - first[newlevel]=w;
1.728 - b=newlevel;
1.729 - }
1.730 - }
1.731 - } // while(true)
1.732 - flow_prop=GEN_FLOW;
1.733 - status=AFTER_PREFLOW_PHASE_2;
1.734 - }
1.735 -
1.736 - /// Returns the value of the maximum flow.
1.737 -
1.738 - /// Returns the value of the maximum flow by returning the excess
1.739 - /// of the target node \c t. This value equals to the value of
1.740 - /// the maximum flow already after running \ref phase1.
1.741 - Num flowValue() const {
1.742 - return excess[_target];
1.743 - }
1.744 -
1.745 -
1.746 - ///Returns a minimum value cut.
1.747 -
1.748 - ///Sets \c M to the characteristic vector of a minimum value
1.749 - ///cut. This method can be called both after running \ref
1.750 - ///phase1 and \ref phase2. It is much faster after
1.751 - ///\ref phase1. \pre M should be a bool-valued node-map. \pre
1.752 - ///If \ref minCut() is called after \ref phase2() then M should
1.753 - ///be initialized to false.
1.754 - template<typename _CutMap>
1.755 - void minCut(_CutMap& M) const {
1.756 - switch ( status ) {
1.757 - case AFTER_PREFLOW_PHASE_1:
1.758 - for(NodeIt v(*_g); v!=INVALID; ++v) {
1.759 - if (level[v] < _node_num) {
1.760 - M.set(v, false);
1.761 - } else {
1.762 - M.set(v, true);
1.763 - }
1.764 - }
1.765 - break;
1.766 - case AFTER_PREFLOW_PHASE_2:
1.767 - minMinCut(M);
1.768 - break;
1.769 - case AFTER_NOTHING:
1.770 - break;
1.771 - }
1.772 - }
1.773 -
1.774 - ///Returns the inclusionwise minimum of the minimum value cuts.
1.775 -
1.776 - ///Sets \c M to the characteristic vector of the minimum value cut
1.777 - ///which is inclusionwise minimum. It is computed by processing a
1.778 - ///bfs from the source node \c s in the residual graph. \pre M
1.779 - ///should be a node map of bools initialized to false. \pre \ref
1.780 - ///phase2 should already be run.
1.781 - template<typename _CutMap>
1.782 - void minMinCut(_CutMap& M) const {
1.783 -
1.784 - std::queue<Node> queue;
1.785 - M.set(_source,true);
1.786 - queue.push(_source);
1.787 -
1.788 - while (!queue.empty()) {
1.789 - Node w=queue.front();
1.790 - queue.pop();
1.791 -
1.792 - for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
1.793 - Node v=_g->target(e);
1.794 - if (!M[v] && _surely.positive((*_capacity)[e] -(*_flow)[e]) ) {
1.795 - queue.push(v);
1.796 - M.set(v, true);
1.797 - }
1.798 - }
1.799 -
1.800 - for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
1.801 - Node v=_g->source(e);
1.802 - if (!M[v] && _surely.positive((*_flow)[e]) ) {
1.803 - queue.push(v);
1.804 - M.set(v, true);
1.805 - }
1.806 - }
1.807 - }
1.808 - }
1.809 -
1.810 - ///Returns the inclusionwise maximum of the minimum value cuts.
1.811 -
1.812 - ///Sets \c M to the characteristic vector of the minimum value cut
1.813 - ///which is inclusionwise maximum. It is computed by processing a
1.814 - ///backward bfs from the target node \c t in the residual graph.
1.815 - ///\pre \ref phase2() or run() should already be run.
1.816 - template<typename _CutMap>
1.817 - void maxMinCut(_CutMap& M) const {
1.818 -
1.819 - for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
1.820 -
1.821 - std::queue<Node> queue;
1.822 -
1.823 - M.set(_target,false);
1.824 - queue.push(_target);
1.825 -
1.826 - while (!queue.empty()) {
1.827 - Node w=queue.front();
1.828 - queue.pop();
1.829 -
1.830 - for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
1.831 - Node v=_g->source(e);
1.832 - if (M[v] && _surely.positive((*_capacity)[e] - (*_flow)[e]) ) {
1.833 - queue.push(v);
1.834 - M.set(v, false);
1.835 - }
1.836 - }
1.837 -
1.838 - for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
1.839 - Node v=_g->target(e);
1.840 - if (M[v] && _surely.positive((*_flow)[e]) ) {
1.841 - queue.push(v);
1.842 - M.set(v, false);
1.843 + for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
1.844 + if (_tolerance.positive((*_capacity)[e])) {
1.845 + Node u = _graph.target(e);
1.846 + if ((*_level)[u] == _level->maxLevel()) continue;
1.847 + _flow->set(e, (*_capacity)[e]);
1.848 + _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
1.849 + if (u != _target && !_level->active(u)) {
1.850 + _level->activate(u);
1.851 }
1.852 }
1.853 }
1.854 }
1.855
1.856 - ///Sets the source node to \c _s.
1.857 + /// \brief Initializes the internal data structures.
1.858 + ///
1.859 + /// Initializes the internal data structures and sets the initial
1.860 + /// flow to the given \c flowMap. The \c flowMap should contain a
1.861 + /// flow or at least a preflow, ie. in each node excluding the
1.862 + /// target the incoming flow should greater or equal to the
1.863 + /// outgoing flow.
1.864 + /// \return %False when the given \c flowMap is not a preflow.
1.865 + template <typename FlowMap>
1.866 + bool flowInit(const FlowMap& flowMap) {
1.867 + createStructures();
1.868
1.869 - ///Sets the source node to \c _s.
1.870 - ///
1.871 - void source(Node _s) {
1.872 - _source=_s;
1.873 - if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
1.874 - status=AFTER_NOTHING;
1.875 - }
1.876 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.877 + _flow->set(e, flowMap[e]);
1.878 + }
1.879
1.880 - ///Returns the source node.
1.881 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.882 + Value excess = 0;
1.883 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.884 + excess += (*_flow)[e];
1.885 + }
1.886 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.887 + excess -= (*_flow)[e];
1.888 + }
1.889 + if (excess < 0 && n != _source) return false;
1.890 + _excess->set(n, excess);
1.891 + }
1.892
1.893 - ///Returns the source node.
1.894 - ///
1.895 - Node source() const {
1.896 - return _source;
1.897 - }
1.898 + typename Graph::template NodeMap<bool> reached(_graph, false);
1.899
1.900 - ///Sets the target node to \c _t.
1.901 + _level->initStart();
1.902 + _level->initAddItem(_target);
1.903
1.904 - ///Sets the target node to \c _t.
1.905 - ///
1.906 - void target(Node _t) {
1.907 - _target=_t;
1.908 - if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
1.909 - status=AFTER_NOTHING;
1.910 - }
1.911 + std::vector<Node> queue;
1.912 + reached.set(_source, true);
1.913
1.914 - ///Returns the target node.
1.915 -
1.916 - ///Returns the target node.
1.917 - ///
1.918 - Node target() const {
1.919 - return _target;
1.920 - }
1.921 -
1.922 - /// Sets the edge map of the capacities to _cap.
1.923 -
1.924 - /// Sets the edge map of the capacities to _cap.
1.925 - ///
1.926 - void capacityMap(const CapacityMap& _cap) {
1.927 - _capacity=&_cap;
1.928 - status=AFTER_NOTHING;
1.929 - }
1.930 - /// Returns a reference to capacity map.
1.931 -
1.932 - /// Returns a reference to capacity map.
1.933 - ///
1.934 - const CapacityMap &capacityMap() const {
1.935 - return *_capacity;
1.936 - }
1.937 -
1.938 - /// Sets the edge map of the flows to _flow.
1.939 -
1.940 - /// Sets the edge map of the flows to _flow.
1.941 - ///
1.942 - void flowMap(FlowMap& _f) {
1.943 - _flow=&_f;
1.944 - flow_prop=NO_FLOW;
1.945 - status=AFTER_NOTHING;
1.946 - }
1.947 -
1.948 - /// Returns a reference to flow map.
1.949 -
1.950 - /// Returns a reference to flow map.
1.951 - ///
1.952 - const FlowMap &flowMap() const {
1.953 - return *_flow;
1.954 - }
1.955 -
1.956 - private:
1.957 -
1.958 - int push(Node w, NNMap& next, VecNode& first) {
1.959 -
1.960 - int lev=level[w];
1.961 - Num exc=excess[w];
1.962 - int newlevel=_node_num; //bound on the next level of w
1.963 -
1.964 - for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
1.965 - if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
1.966 - Node v=_g->target(e);
1.967 -
1.968 - if( lev > level[v] ) { //Push is allowed now
1.969 -
1.970 - if ( excess[v] == 0 && v!=_target && v!=_source ) {
1.971 - next.set(v,first[level[v]]);
1.972 - first[level[v]]=v;
1.973 - }
1.974 -
1.975 - Num cap=(*_capacity)[e];
1.976 - Num flo=(*_flow)[e];
1.977 - Num remcap=cap-flo;
1.978 -
1.979 - if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push.
1.980 -
1.981 - _flow->set(e, flo+exc);
1.982 - excess.set(v, excess[v]+exc);
1.983 - exc=0;
1.984 - break;
1.985 -
1.986 - } else { //A saturating push.
1.987 - _flow->set(e, cap);
1.988 - excess.set(v, excess[v]+remcap);
1.989 - exc-=remcap;
1.990 - }
1.991 - } else if ( newlevel > level[v] ) newlevel = level[v];
1.992 - } //for out edges wv
1.993 -
1.994 - if ( exc != 0 ) {
1.995 - for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
1.996 -
1.997 - if ( !_surely.positive((*_flow)[e]) ) continue;
1.998 - Node v=_g->source(e);
1.999 -
1.1000 - if( lev > level[v] ) { //Push is allowed now
1.1001 -
1.1002 - if ( excess[v] == 0 && v!=_target && v!=_source ) {
1.1003 - next.set(v,first[level[v]]);
1.1004 - first[level[v]]=v;
1.1005 - }
1.1006 -
1.1007 - Num flo=(*_flow)[e];
1.1008 -
1.1009 - if ( !_surely.less(flo, exc) ) { //A nonsaturating push.
1.1010 -
1.1011 - _flow->set(e, flo-exc);
1.1012 - excess.set(v, excess[v]+exc);
1.1013 - exc=0;
1.1014 - break;
1.1015 - } else { //A saturating push.
1.1016 -
1.1017 - excess.set(v, excess[v]+flo);
1.1018 - exc-=flo;
1.1019 - _flow->set(e,0);
1.1020 - }
1.1021 - } else if ( newlevel > level[v] ) newlevel = level[v];
1.1022 - } //for in edges vw
1.1023 -
1.1024 - } // if w still has excess after the out edge for cycle
1.1025 -
1.1026 - excess.set(w, exc);
1.1027 -
1.1028 - return newlevel;
1.1029 - }
1.1030 -
1.1031 -
1.1032 -
1.1033 - void preflowPreproc(VecNode& first, NNMap& next,
1.1034 - VecNode& level_list, NNMap& left, NNMap& right)
1.1035 - {
1.1036 - for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
1.1037 - std::queue<Node> bfs_queue;
1.1038 -
1.1039 - if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
1.1040 - //Reverse_bfs from t in the residual graph,
1.1041 - //to find the starting level.
1.1042 - level.set(_target,0);
1.1043 - bfs_queue.push(_target);
1.1044 -
1.1045 - while ( !bfs_queue.empty() ) {
1.1046 -
1.1047 - Node v=bfs_queue.front();
1.1048 - bfs_queue.pop();
1.1049 - int l=level[v]+1;
1.1050 -
1.1051 - for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
1.1052 - if ( !_surely.positive((*_capacity)[e] - (*_flow)[e] )) continue;
1.1053 - Node w=_g->source(e);
1.1054 - if ( level[w] == _node_num && w != _source ) {
1.1055 - bfs_queue.push(w);
1.1056 - Node z=level_list[l];
1.1057 - if ( z!=INVALID ) left.set(z,w);
1.1058 - right.set(w,z);
1.1059 - level_list[l]=w;
1.1060 - level.set(w, l);
1.1061 + queue.push_back(_target);
1.1062 + reached.set(_target, true);
1.1063 + while (!queue.empty()) {
1.1064 + _level->initNewLevel();
1.1065 + std::vector<Node> nqueue;
1.1066 + for (int i = 0; i < int(queue.size()); ++i) {
1.1067 + Node n = queue[i];
1.1068 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1069 + Node u = _graph.source(e);
1.1070 + if (!reached[u] &&
1.1071 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
1.1072 + reached.set(u, true);
1.1073 + _level->initAddItem(u);
1.1074 + nqueue.push_back(u);
1.1075 }
1.1076 }
1.1077 -
1.1078 - for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
1.1079 - if ( !_surely.positive((*_flow)[e]) ) continue;
1.1080 - Node w=_g->target(e);
1.1081 - if ( level[w] == _node_num && w != _source ) {
1.1082 - bfs_queue.push(w);
1.1083 - Node z=level_list[l];
1.1084 - if ( z!=INVALID ) left.set(z,w);
1.1085 - right.set(w,z);
1.1086 - level_list[l]=w;
1.1087 - level.set(w, l);
1.1088 - }
1.1089 - }
1.1090 - } //while
1.1091 - } //if
1.1092 -
1.1093 -
1.1094 - switch (flow_prop) {
1.1095 - case NO_FLOW:
1.1096 - for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
1.1097 - case ZERO_FLOW:
1.1098 - for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
1.1099 -
1.1100 - //Reverse_bfs from t, to find the starting level.
1.1101 - level.set(_target,0);
1.1102 - bfs_queue.push(_target);
1.1103 -
1.1104 - while ( !bfs_queue.empty() ) {
1.1105 -
1.1106 - Node v=bfs_queue.front();
1.1107 - bfs_queue.pop();
1.1108 - int l=level[v]+1;
1.1109 -
1.1110 - for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
1.1111 - Node w=_g->source(e);
1.1112 - if ( level[w] == _node_num && w != _source ) {
1.1113 - bfs_queue.push(w);
1.1114 - Node z=level_list[l];
1.1115 - if ( z!=INVALID ) left.set(z,w);
1.1116 - right.set(w,z);
1.1117 - level_list[l]=w;
1.1118 - level.set(w, l);
1.1119 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1120 + Node v = _graph.target(e);
1.1121 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
1.1122 + reached.set(v, true);
1.1123 + _level->initAddItem(v);
1.1124 + nqueue.push_back(v);
1.1125 }
1.1126 }
1.1127 }
1.1128 -
1.1129 - //the starting flow
1.1130 - for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
1.1131 - Num c=(*_capacity)[e];
1.1132 - if ( !_surely.positive(c) ) continue;
1.1133 - Node w=_g->target(e);
1.1134 - if ( level[w] < _node_num ) {
1.1135 - if ( excess[w] == 0 && w!=_target ) { //putting into the stack
1.1136 - next.set(w,first[level[w]]);
1.1137 - first[level[w]]=w;
1.1138 - }
1.1139 - _flow->set(e, c);
1.1140 - excess.set(w, excess[w]+c);
1.1141 + queue.swap(nqueue);
1.1142 + }
1.1143 + _level->initFinish();
1.1144 +
1.1145 + for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
1.1146 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.1147 + if (_tolerance.positive(rem)) {
1.1148 + Node u = _graph.target(e);
1.1149 + if ((*_level)[u] == _level->maxLevel()) continue;
1.1150 + _flow->set(e, (*_capacity)[e]);
1.1151 + _excess->set(u, (*_excess)[u] + rem);
1.1152 + if (u != _target && !_level->active(u)) {
1.1153 + _level->activate(u);
1.1154 }
1.1155 }
1.1156 - break;
1.1157 + }
1.1158 + for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
1.1159 + Value rem = (*_flow)[e];
1.1160 + if (_tolerance.positive(rem)) {
1.1161 + Node v = _graph.source(e);
1.1162 + if ((*_level)[v] == _level->maxLevel()) continue;
1.1163 + _flow->set(e, 0);
1.1164 + _excess->set(v, (*_excess)[v] + rem);
1.1165 + if (v != _target && !_level->active(v)) {
1.1166 + _level->activate(v);
1.1167 + }
1.1168 + }
1.1169 + }
1.1170 + return true;
1.1171 + }
1.1172
1.1173 - case GEN_FLOW:
1.1174 - for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
1.1175 - {
1.1176 - Num exc=0;
1.1177 - for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
1.1178 - for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
1.1179 - if (!_surely.positive(exc)) {
1.1180 - exc = 0;
1.1181 - }
1.1182 - excess.set(_target,exc);
1.1183 + /// \brief Starts the first phase of the preflow algorithm.
1.1184 + ///
1.1185 + /// The preflow algorithm consists of two phases, this method runs
1.1186 + /// the first phase. After the first phase the maximum flow value
1.1187 + /// and a minimum value cut can already be computed, although a
1.1188 + /// maximum flow is not yet obtained. So after calling this method
1.1189 + /// \ref flowValue() returns the value of a maximum flow and \ref
1.1190 + /// minCut() returns a minimum cut.
1.1191 + /// \pre One of the \ref init() functions should be called.
1.1192 + void startFirstPhase() {
1.1193 + _phase = true;
1.1194 +
1.1195 + Node n = _level->highestActive();
1.1196 + int level = _level->highestActiveLevel();
1.1197 + while (n != INVALID) {
1.1198 + int num = _node_num;
1.1199 +
1.1200 + while (num > 0 && n != INVALID) {
1.1201 + Value excess = (*_excess)[n];
1.1202 + int new_level = _level->maxLevel();
1.1203 +
1.1204 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1205 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.1206 + if (!_tolerance.positive(rem)) continue;
1.1207 + Node v = _graph.target(e);
1.1208 + if ((*_level)[v] < level) {
1.1209 + if (!_level->active(v) && v != _target) {
1.1210 + _level->activate(v);
1.1211 + }
1.1212 + if (!_tolerance.less(rem, excess)) {
1.1213 + _flow->set(e, (*_flow)[e] + excess);
1.1214 + _excess->set(v, (*_excess)[v] + excess);
1.1215 + excess = 0;
1.1216 + goto no_more_push_1;
1.1217 + } else {
1.1218 + excess -= rem;
1.1219 + _excess->set(v, (*_excess)[v] + rem);
1.1220 + _flow->set(e, (*_capacity)[e]);
1.1221 + }
1.1222 + } else if (new_level > (*_level)[v]) {
1.1223 + new_level = (*_level)[v];
1.1224 + }
1.1225 + }
1.1226 +
1.1227 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1228 + Value rem = (*_flow)[e];
1.1229 + if (!_tolerance.positive(rem)) continue;
1.1230 + Node v = _graph.source(e);
1.1231 + if ((*_level)[v] < level) {
1.1232 + if (!_level->active(v) && v != _target) {
1.1233 + _level->activate(v);
1.1234 + }
1.1235 + if (!_tolerance.less(rem, excess)) {
1.1236 + _flow->set(e, (*_flow)[e] - excess);
1.1237 + _excess->set(v, (*_excess)[v] + excess);
1.1238 + excess = 0;
1.1239 + goto no_more_push_1;
1.1240 + } else {
1.1241 + excess -= rem;
1.1242 + _excess->set(v, (*_excess)[v] + rem);
1.1243 + _flow->set(e, 0);
1.1244 + }
1.1245 + } else if (new_level > (*_level)[v]) {
1.1246 + new_level = (*_level)[v];
1.1247 + }
1.1248 + }
1.1249 +
1.1250 + no_more_push_1:
1.1251 +
1.1252 + _excess->set(n, excess);
1.1253 +
1.1254 + if (excess != 0) {
1.1255 + if (new_level + 1 < _level->maxLevel()) {
1.1256 + _level->liftHighestActive(new_level + 1);
1.1257 + } else {
1.1258 + _level->liftHighestActiveToTop();
1.1259 + }
1.1260 + if (_level->emptyLevel(level)) {
1.1261 + _level->liftToTop(level);
1.1262 + }
1.1263 + } else {
1.1264 + _level->deactivate(n);
1.1265 + }
1.1266 +
1.1267 + n = _level->highestActive();
1.1268 + level = _level->highestActiveLevel();
1.1269 + --num;
1.1270 + }
1.1271 +
1.1272 + num = _node_num * 20;
1.1273 + while (num > 0 && n != INVALID) {
1.1274 + Value excess = (*_excess)[n];
1.1275 + int new_level = _level->maxLevel();
1.1276 +
1.1277 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1278 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.1279 + if (!_tolerance.positive(rem)) continue;
1.1280 + Node v = _graph.target(e);
1.1281 + if ((*_level)[v] < level) {
1.1282 + if (!_level->active(v) && v != _target) {
1.1283 + _level->activate(v);
1.1284 + }
1.1285 + if (!_tolerance.less(rem, excess)) {
1.1286 + _flow->set(e, (*_flow)[e] + excess);
1.1287 + _excess->set(v, (*_excess)[v] + excess);
1.1288 + excess = 0;
1.1289 + goto no_more_push_2;
1.1290 + } else {
1.1291 + excess -= rem;
1.1292 + _excess->set(v, (*_excess)[v] + rem);
1.1293 + _flow->set(e, (*_capacity)[e]);
1.1294 + }
1.1295 + } else if (new_level > (*_level)[v]) {
1.1296 + new_level = (*_level)[v];
1.1297 + }
1.1298 + }
1.1299 +
1.1300 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1301 + Value rem = (*_flow)[e];
1.1302 + if (!_tolerance.positive(rem)) continue;
1.1303 + Node v = _graph.source(e);
1.1304 + if ((*_level)[v] < level) {
1.1305 + if (!_level->active(v) && v != _target) {
1.1306 + _level->activate(v);
1.1307 + }
1.1308 + if (!_tolerance.less(rem, excess)) {
1.1309 + _flow->set(e, (*_flow)[e] - excess);
1.1310 + _excess->set(v, (*_excess)[v] + excess);
1.1311 + excess = 0;
1.1312 + goto no_more_push_2;
1.1313 + } else {
1.1314 + excess -= rem;
1.1315 + _excess->set(v, (*_excess)[v] + rem);
1.1316 + _flow->set(e, 0);
1.1317 + }
1.1318 + } else if (new_level > (*_level)[v]) {
1.1319 + new_level = (*_level)[v];
1.1320 + }
1.1321 + }
1.1322 +
1.1323 + no_more_push_2:
1.1324 +
1.1325 + _excess->set(n, excess);
1.1326 +
1.1327 + if (excess != 0) {
1.1328 + if (new_level + 1 < _level->maxLevel()) {
1.1329 + _level->liftActiveOn(level, new_level + 1);
1.1330 + } else {
1.1331 + _level->liftActiveToTop(level);
1.1332 + }
1.1333 + if (_level->emptyLevel(level)) {
1.1334 + _level->liftToTop(level);
1.1335 + }
1.1336 + } else {
1.1337 + _level->deactivate(n);
1.1338 + }
1.1339 +
1.1340 + while (level >= 0 && _level->activeFree(level)) {
1.1341 + --level;
1.1342 + }
1.1343 + if (level == -1) {
1.1344 + n = _level->highestActive();
1.1345 + level = _level->highestActiveLevel();
1.1346 + } else {
1.1347 + n = _level->activeOn(level);
1.1348 + }
1.1349 + --num;
1.1350 + }
1.1351 + }
1.1352 + }
1.1353 +
1.1354 + /// \brief Starts the second phase of the preflow algorithm.
1.1355 + ///
1.1356 + /// The preflow algorithm consists of two phases, this method runs
1.1357 + /// the second phase. After calling \ref init() and \ref
1.1358 + /// startFirstPhase() and then \ref startSecondPhase(), \ref
1.1359 + /// flowMap() return a maximum flow, \ref flowValue() returns the
1.1360 + /// value of a maximum flow, \ref minCut() returns a minimum cut
1.1361 + /// \pre The \ref init() and startFirstPhase() functions should be
1.1362 + /// called before.
1.1363 + void startSecondPhase() {
1.1364 + _phase = false;
1.1365 +
1.1366 + typename Graph::template NodeMap<bool> reached(_graph, false);
1.1367 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1368 + reached.set(n, (*_level)[n] < _level->maxLevel());
1.1369 + }
1.1370 +
1.1371 + _level->initStart();
1.1372 + _level->initAddItem(_source);
1.1373 +
1.1374 + std::vector<Node> queue;
1.1375 + queue.push_back(_source);
1.1376 + reached.set(_source, true);
1.1377 +
1.1378 + while (!queue.empty()) {
1.1379 + _level->initNewLevel();
1.1380 + std::vector<Node> nqueue;
1.1381 + for (int i = 0; i < int(queue.size()); ++i) {
1.1382 + Node n = queue[i];
1.1383 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1384 + Node v = _graph.target(e);
1.1385 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
1.1386 + reached.set(v, true);
1.1387 + _level->initAddItem(v);
1.1388 + nqueue.push_back(v);
1.1389 + }
1.1390 + }
1.1391 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1392 + Node u = _graph.source(e);
1.1393 + if (!reached[u] &&
1.1394 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
1.1395 + reached.set(u, true);
1.1396 + _level->initAddItem(u);
1.1397 + nqueue.push_back(u);
1.1398 + }
1.1399 + }
1.1400 + }
1.1401 + queue.swap(nqueue);
1.1402 + }
1.1403 + _level->initFinish();
1.1404 +
1.1405 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1406 + if ((*_excess)[n] > 0 && _target != n) {
1.1407 + _level->activate(n);
1.1408 + }
1.1409 + }
1.1410 +
1.1411 + Node n;
1.1412 + while ((n = _level->highestActive()) != INVALID) {
1.1413 + Value excess = (*_excess)[n];
1.1414 + int level = _level->highestActiveLevel();
1.1415 + int new_level = _level->maxLevel();
1.1416 +
1.1417 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1418 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.1419 + if (!_tolerance.positive(rem)) continue;
1.1420 + Node v = _graph.target(e);
1.1421 + if ((*_level)[v] < level) {
1.1422 + if (!_level->active(v) && v != _source) {
1.1423 + _level->activate(v);
1.1424 + }
1.1425 + if (!_tolerance.less(rem, excess)) {
1.1426 + _flow->set(e, (*_flow)[e] + excess);
1.1427 + _excess->set(v, (*_excess)[v] + excess);
1.1428 + excess = 0;
1.1429 + goto no_more_push;
1.1430 + } else {
1.1431 + excess -= rem;
1.1432 + _excess->set(v, (*_excess)[v] + rem);
1.1433 + _flow->set(e, (*_capacity)[e]);
1.1434 + }
1.1435 + } else if (new_level > (*_level)[v]) {
1.1436 + new_level = (*_level)[v];
1.1437 + }
1.1438 }
1.1439
1.1440 - //the starting flow
1.1441 - for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) {
1.1442 - Num rem=(*_capacity)[e]-(*_flow)[e];
1.1443 - if ( !_surely.positive(rem) ) continue;
1.1444 - Node w=_g->target(e);
1.1445 - if ( level[w] < _node_num ) {
1.1446 - if ( excess[w] == 0 && w!=_target ) { //putting into the stack
1.1447 - next.set(w,first[level[w]]);
1.1448 - first[level[w]]=w;
1.1449 - }
1.1450 - _flow->set(e, (*_capacity)[e]);
1.1451 - excess.set(w, excess[w]+rem);
1.1452 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.1453 + Value rem = (*_flow)[e];
1.1454 + if (!_tolerance.positive(rem)) continue;
1.1455 + Node v = _graph.source(e);
1.1456 + if ((*_level)[v] < level) {
1.1457 + if (!_level->active(v) && v != _source) {
1.1458 + _level->activate(v);
1.1459 + }
1.1460 + if (!_tolerance.less(rem, excess)) {
1.1461 + _flow->set(e, (*_flow)[e] - excess);
1.1462 + _excess->set(v, (*_excess)[v] + excess);
1.1463 + excess = 0;
1.1464 + goto no_more_push;
1.1465 + } else {
1.1466 + excess -= rem;
1.1467 + _excess->set(v, (*_excess)[v] + rem);
1.1468 + _flow->set(e, 0);
1.1469 + }
1.1470 + } else if (new_level > (*_level)[v]) {
1.1471 + new_level = (*_level)[v];
1.1472 }
1.1473 }
1.1474 -
1.1475 - for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
1.1476 - if ( !_surely.positive((*_flow)[e]) ) continue;
1.1477 - Node w=_g->source(e);
1.1478 - if ( level[w] < _node_num ) {
1.1479 - if ( excess[w] == 0 && w!=_target ) {
1.1480 - next.set(w,first[level[w]]);
1.1481 - first[level[w]]=w;
1.1482 - }
1.1483 - excess.set(w, excess[w]+(*_flow)[e]);
1.1484 - _flow->set(e, 0);
1.1485 +
1.1486 + no_more_push:
1.1487 +
1.1488 + _excess->set(n, excess);
1.1489 +
1.1490 + if (excess != 0) {
1.1491 + if (new_level + 1 < _level->maxLevel()) {
1.1492 + _level->liftHighestActive(new_level + 1);
1.1493 + } else {
1.1494 + // Calculation error
1.1495 + _level->liftHighestActiveToTop();
1.1496 }
1.1497 - }
1.1498 - break;
1.1499 -
1.1500 - case PRE_FLOW:
1.1501 - //the starting flow
1.1502 - for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
1.1503 - Num rem=(*_capacity)[e]-(*_flow)[e];
1.1504 - if ( !_surely.positive(rem) ) continue;
1.1505 - Node w=_g->target(e);
1.1506 - if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
1.1507 - }
1.1508 -
1.1509 - for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
1.1510 - if ( !_surely.positive((*_flow)[e]) ) continue;
1.1511 - Node w=_g->source(e);
1.1512 - if ( level[w] < _node_num ) _flow->set(e, 0);
1.1513 - }
1.1514 -
1.1515 - //computing the excess
1.1516 - for(NodeIt w(*_g); w!=INVALID; ++w) {
1.1517 - Num exc=0;
1.1518 - for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
1.1519 - for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
1.1520 - if (!_surely.positive(exc)) {
1.1521 - exc = 0;
1.1522 - }
1.1523 - excess.set(w,exc);
1.1524 -
1.1525 - //putting the active nodes into the stack
1.1526 - int lev=level[w];
1.1527 - if ( exc != 0 && lev < _node_num && Node(w) != _target ) {
1.1528 - next.set(w,first[lev]);
1.1529 - first[lev]=w;
1.1530 - }
1.1531 - }
1.1532 - break;
1.1533 - } //switch
1.1534 - } //preflowPreproc
1.1535 -
1.1536 -
1.1537 - void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
1.1538 - VecNode& level_list, NNMap& left,
1.1539 - NNMap& right, int& b, int& k, bool what_heur )
1.1540 - {
1.1541 -
1.1542 - int lev=level[w];
1.1543 -
1.1544 - Node right_n=right[w];
1.1545 - Node left_n=left[w];
1.1546 -
1.1547 - //unlacing starts
1.1548 - if ( right_n!=INVALID ) {
1.1549 - if ( left_n!=INVALID ) {
1.1550 - right.set(left_n, right_n);
1.1551 - left.set(right_n, left_n);
1.1552 + if (_level->emptyLevel(level)) {
1.1553 + // Calculation error
1.1554 + _level->liftToTop(level);
1.1555 + }
1.1556 } else {
1.1557 - level_list[lev]=right_n;
1.1558 - left.set(right_n, INVALID);
1.1559 - }
1.1560 - } else {
1.1561 - if ( left_n!=INVALID ) {
1.1562 - right.set(left_n, INVALID);
1.1563 - } else {
1.1564 - level_list[lev]=INVALID;
1.1565 - }
1.1566 - }
1.1567 - //unlacing ends
1.1568 -
1.1569 - if ( level_list[lev]==INVALID ) {
1.1570 -
1.1571 - //gapping starts
1.1572 - for (int i=lev; i!=k ; ) {
1.1573 - Node v=level_list[++i];
1.1574 - while ( v!=INVALID ) {
1.1575 - level.set(v,_node_num);
1.1576 - v=right[v];
1.1577 - }
1.1578 - level_list[i]=INVALID;
1.1579 - if ( !what_heur ) first[i]=INVALID;
1.1580 + _level->deactivate(n);
1.1581 }
1.1582
1.1583 - level.set(w,_node_num);
1.1584 - b=lev-1;
1.1585 - k=b;
1.1586 - //gapping ends
1.1587 + }
1.1588 + }
1.1589
1.1590 - } else {
1.1591 + /// \brief Runs the preflow algorithm.
1.1592 + ///
1.1593 + /// Runs the preflow algorithm.
1.1594 + /// \note pf.run() is just a shortcut of the following code.
1.1595 + /// \code
1.1596 + /// pf.init();
1.1597 + /// pf.startFirstPhase();
1.1598 + /// pf.startSecondPhase();
1.1599 + /// \endcode
1.1600 + void run() {
1.1601 + init();
1.1602 + startFirstPhase();
1.1603 + startSecondPhase();
1.1604 + }
1.1605
1.1606 - if ( newlevel == _node_num ) level.set(w,_node_num);
1.1607 - else {
1.1608 - level.set(w,++newlevel);
1.1609 - next.set(w,first[newlevel]);
1.1610 - first[newlevel]=w;
1.1611 - if ( what_heur ) b=newlevel;
1.1612 - if ( k < newlevel ) ++k; //now k=newlevel
1.1613 - Node z=level_list[newlevel];
1.1614 - if ( z!=INVALID ) left.set(z,w);
1.1615 - right.set(w,z);
1.1616 - left.set(w,INVALID);
1.1617 - level_list[newlevel]=w;
1.1618 - }
1.1619 + /// \brief Runs the preflow algorithm to compute the minimum cut.
1.1620 + ///
1.1621 + /// Runs the preflow algorithm to compute the minimum cut.
1.1622 + /// \note pf.runMinCut() is just a shortcut of the following code.
1.1623 + /// \code
1.1624 + /// pf.init();
1.1625 + /// pf.startFirstPhase();
1.1626 + /// \endcode
1.1627 + void runMinCut() {
1.1628 + init();
1.1629 + startFirstPhase();
1.1630 + }
1.1631 +
1.1632 + /// @}
1.1633 +
1.1634 + /// \name Query Functions
1.1635 + /// The result of the %Dijkstra algorithm can be obtained using these
1.1636 + /// functions.\n
1.1637 + /// Before the use of these functions,
1.1638 + /// either run() or start() must be called.
1.1639 +
1.1640 + ///@{
1.1641 +
1.1642 + /// \brief Returns the value of the maximum flow.
1.1643 + ///
1.1644 + /// Returns the value of the maximum flow by returning the excess
1.1645 + /// of the target node \c t. This value equals to the value of
1.1646 + /// the maximum flow already after the first phase.
1.1647 + Value flowValue() const {
1.1648 + return (*_excess)[_target];
1.1649 + }
1.1650 +
1.1651 + /// \brief Returns true when the node is on the source side of minimum cut.
1.1652 + ///
1.1653 + /// Returns true when the node is on the source side of minimum
1.1654 + /// cut. This method can be called both after running \ref
1.1655 + /// startFirstPhase() and \ref startSecondPhase().
1.1656 + bool minCut(const Node& node) const {
1.1657 + return ((*_level)[node] == _level->maxLevel()) == _phase;
1.1658 + }
1.1659 +
1.1660 + /// \brief Returns a minimum value cut.
1.1661 + ///
1.1662 + /// Sets the \c cutMap to the characteristic vector of a minimum value
1.1663 + /// cut. This method can be called both after running \ref
1.1664 + /// startFirstPhase() and \ref startSecondPhase(). The result after second
1.1665 + /// phase could be changed slightly if inexact computation is used.
1.1666 + /// \pre The \c cutMap should be a bool-valued node-map.
1.1667 + template <typename CutMap>
1.1668 + void minCutMap(CutMap& cutMap) const {
1.1669 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1670 + cutMap.set(n, minCut(n));
1.1671 }
1.1672 - } //relabel
1.1673 + }
1.1674
1.1675 - };
1.1676 + /// \brief Returns the flow on the edge.
1.1677 + ///
1.1678 + /// Sets the \c flowMap to the flow on the edges. This method can
1.1679 + /// be called after the second phase of algorithm.
1.1680 + Value flow(const Edge& edge) const {
1.1681 + return (*_flow)[edge];
1.1682 + }
1.1683 +
1.1684 + /// @}
1.1685 + };
1.1686 +}
1.1687
1.1688 - ///\ingroup max_flow
1.1689 - ///\brief Function type interface for Preflow algorithm.
1.1690 - ///
1.1691 - ///Function type interface for Preflow algorithm.
1.1692 - ///\sa Preflow
1.1693 - template<class GR, class CM, class FM>
1.1694 - Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
1.1695 - typename GR::Node source,
1.1696 - typename GR::Node target,
1.1697 - const CM &cap,
1.1698 - FM &flow
1.1699 - )
1.1700 - {
1.1701 - return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
1.1702 - }
1.1703 -
1.1704 -} //namespace lemon
1.1705 -
1.1706 -#endif //LEMON_PREFLOW_H
1.1707 +#endif