Revert to long long int since currently I don't know a better solution.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
19 #ifndef LEMON_CIRCULATION_H
20 #define LEMON_CIRCULATION_H
22 #include <lemon/graph_utils.h>
25 #include <lemon/tolerance.h>
26 #include <lemon/elevator.h>
30 ///\brief Push-prelabel algorithm for finding a feasible circulation.
34 /// \brief Default traits class of Circulation class.
36 /// Default traits class of Circulation class.
37 /// \param _Graph Graph type.
38 /// \param _CapacityMap Type of capacity map.
39 template <typename _Graph, typename _LCapMap,
40 typename _UCapMap, typename _DeltaMap>
41 struct CirculationDefaultTraits {
43 /// \brief The graph type the algorithm runs on.
46 /// \brief The type of the map that stores the circulation lower
49 /// The type of the map that stores the circulation lower bound.
50 /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
51 typedef _LCapMap LCapMap;
53 /// \brief The type of the map that stores the circulation upper
56 /// The type of the map that stores the circulation upper bound.
57 /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
58 typedef _UCapMap UCapMap;
60 /// \brief The type of the map that stores the upper bound of
63 /// The type of the map that stores the lower bound of node
64 /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
66 typedef _DeltaMap DeltaMap;
68 /// \brief The type of the length of the edges.
69 typedef typename DeltaMap::Value Value;
71 /// \brief The map type that stores the flow values.
73 /// The map type that stores the flow values.
74 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
75 typedef typename Graph::template EdgeMap<Value> FlowMap;
77 /// \brief Instantiates a FlowMap.
79 /// This function instantiates a \ref FlowMap.
80 /// \param graph The graph, to which we would like to define the flow map.
81 static FlowMap* createFlowMap(const Graph& graph) {
82 return new FlowMap(graph);
85 /// \brief The eleavator type used by Circulation algorithm.
87 /// The elevator type used by Circulation algorithm.
90 /// \sa LinkedElevator
91 typedef Elevator<Graph, typename Graph::Node> Elevator;
93 /// \brief Instantiates an Elevator.
95 /// This function instantiates a \ref Elevator.
96 /// \param graph The graph, to which we would like to define the elevator.
97 /// \param max_level The maximum level of the elevator.
98 static Elevator* createElevator(const Graph& graph, int max_level) {
99 return new Elevator(graph, max_level);
102 /// \brief The tolerance used by the algorithm
104 /// The tolerance used by the algorithm to handle inexact computation.
105 typedef Tolerance<Value> Tolerance;
109 ///Push-relabel algorithm for the Network Circulation Problem.
112 ///This class implements a push-relabel algorithm
113 ///for the Network Circulation Problem.
114 ///The exact formulation of this problem is the following.
115 /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq -delta(v)\quad \forall v\in V \f]
116 /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
118 template<class _Graph,
119 class _LCapMap=typename _Graph::template EdgeMap<int>,
120 class _UCapMap=_LCapMap,
121 class _DeltaMap=typename _Graph::template NodeMap<
122 typename _UCapMap::Value>,
123 class _Traits=CirculationDefaultTraits<_Graph, _LCapMap,
124 _UCapMap, _DeltaMap> >
127 typedef _Traits Traits;
128 typedef typename Traits::Graph Graph;
129 GRAPH_TYPEDEFS(typename Graph);
131 typedef typename Traits::Value Value;
133 typedef typename Traits::LCapMap LCapMap;
134 typedef typename Traits::UCapMap UCapMap;
135 typedef typename Traits::DeltaMap DeltaMap;
136 typedef typename Traits::FlowMap FlowMap;
137 typedef typename Traits::Elevator Elevator;
138 typedef typename Traits::Tolerance Tolerance;
140 typedef typename Graph::template NodeMap<Value> ExcessMap;
147 const DeltaMap *_delta;
162 typedef Circulation Create;
164 ///\name Named template parameters
168 template <typename _FlowMap>
169 struct DefFlowMapTraits : public Traits {
170 typedef _FlowMap FlowMap;
171 static FlowMap *createFlowMap(const Graph&) {
172 throw UninitializedParameter();
176 /// \brief \ref named-templ-param "Named parameter" for setting
179 /// \ref named-templ-param "Named parameter" for setting FlowMap
181 template <typename _FlowMap>
183 : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
184 DefFlowMapTraits<_FlowMap> > {
185 typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
186 DefFlowMapTraits<_FlowMap> > Create;
189 template <typename _Elevator>
190 struct DefElevatorTraits : public Traits {
191 typedef _Elevator Elevator;
192 static Elevator *createElevator(const Graph&, int) {
193 throw UninitializedParameter();
197 /// \brief \ref named-templ-param "Named parameter" for setting
200 /// \ref named-templ-param "Named parameter" for setting Elevator
202 template <typename _Elevator>
204 : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
205 DefElevatorTraits<_Elevator> > {
206 typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
207 DefElevatorTraits<_Elevator> > Create;
210 template <typename _Elevator>
211 struct DefStandardElevatorTraits : public Traits {
212 typedef _Elevator Elevator;
213 static Elevator *createElevator(const Graph& graph, int max_level) {
214 return new Elevator(graph, max_level);
218 /// \brief \ref named-templ-param "Named parameter" for setting
221 /// \ref named-templ-param "Named parameter" for setting Elevator
222 /// type. The Elevator should be standard constructor interface, ie.
223 /// the graph and the maximum level should be passed to it.
224 template <typename _Elevator>
225 struct DefStandardElevator
226 : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
227 DefStandardElevatorTraits<_Elevator> > {
228 typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
229 DefStandardElevatorTraits<_Elevator> > Create;
240 /// The constructor of the class.
242 /// The constructor of the class.
243 /// \param g The directed graph the algorithm runs on.
244 /// \param lo The lower bound capacity of the edges.
245 /// \param up The upper bound capacity of the edges.
246 /// \param delta The lower bound on node excess.
247 Circulation(const Graph &g,const LCapMap &lo,
248 const UCapMap &up,const DeltaMap &delta)
249 : _g(g), _node_num(),
250 _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
251 _level(0), _local_level(false), _excess(0), _el() {}
260 void createStructures() {
261 _node_num = _el = countNodes(_g);
264 _flow = Traits::createFlowMap(_g);
268 _level = Traits::createElevator(_g, _node_num);
272 _excess = new ExcessMap(_g);
276 void destroyStructures() {
290 /// Sets the lower bound capacity map.
292 /// Sets the lower bound capacity map.
293 /// \return \c (*this)
294 Circulation& lowerCapMap(const LCapMap& map) {
299 /// Sets the upper bound capacity map.
301 /// Sets the upper bound capacity map.
302 /// \return \c (*this)
303 Circulation& upperCapMap(const LCapMap& map) {
308 /// Sets the lower bound map on excess.
310 /// Sets the lower bound map on excess.
311 /// \return \c (*this)
312 Circulation& deltaMap(const DeltaMap& map) {
317 /// Sets the flow map.
319 /// Sets the flow map.
320 /// \return \c (*this)
321 Circulation& flowMap(FlowMap& map) {
330 /// Returns the flow map.
332 /// \return The flow map.
334 const FlowMap& flowMap() {
338 /// Sets the elevator.
340 /// Sets the elevator.
341 /// \return \c (*this)
342 Circulation& elevator(Elevator& elevator) {
345 _local_level = false;
351 /// Returns the elevator.
353 /// \return The elevator.
355 const Elevator& elevator() {
359 /// Sets the tolerance used by algorithm.
361 /// Sets the tolerance used by algorithm.
363 Circulation& tolerance(const Tolerance& tolerance) const {
368 /// Returns the tolerance used by algorithm.
370 /// Returns the tolerance used by algorithm.
372 const Tolerance& tolerance() const {
376 /// \name Execution control
377 /// The simplest way to execute the algorithm is to use one of the
378 /// member functions called \c run().
380 /// If you need more control on initial solution or execution then
381 /// you have to call one \ref init() function and then the start()
386 /// Initializes the internal data structures.
388 /// Initializes the internal data structures. This function sets
389 /// all flow values to the lower bound.
390 /// \return This function returns false if the initialization
391 /// process found a barrier.
396 for(NodeIt n(_g);n!=INVALID;++n) {
397 _excess->set(n, (*_delta)[n]);
400 for (EdgeIt e(_g);e!=INVALID;++e) {
401 _flow->set(e, (*_lo)[e]);
402 _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
403 _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
406 // global relabeling tested, but in general case it provides
407 // worse performance for random graphs
409 for(NodeIt n(_g);n!=INVALID;++n)
410 _level->initAddItem(n);
411 _level->initFinish();
412 for(NodeIt n(_g);n!=INVALID;++n)
413 if(_tol.positive((*_excess)[n]))
417 /// Initializes the internal data structures.
419 /// Initializes the internal data structures. This functions uses
420 /// greedy approach to construct the initial solution.
425 for(NodeIt n(_g);n!=INVALID;++n) {
426 _excess->set(n, (*_delta)[n]);
429 for (EdgeIt e(_g);e!=INVALID;++e) {
430 if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
431 _flow->set(e, (*_up)[e]);
432 _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
433 _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
434 } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
435 _flow->set(e, (*_lo)[e]);
436 _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
437 _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
439 Value fc = -(*_excess)[_g.target(e)];
441 _excess->set(_g.target(e), 0);
442 _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
447 for(NodeIt n(_g);n!=INVALID;++n)
448 _level->initAddItem(n);
449 _level->initFinish();
450 for(NodeIt n(_g);n!=INVALID;++n)
451 if(_tol.positive((*_excess)[n]))
455 ///Starts the algorithm
457 ///This function starts the algorithm.
458 ///\return This function returns true if it found a feasible circulation.
466 Node last_activated=INVALID;
467 while((act=_level->highestActive())!=INVALID) {
468 int actlevel=(*_level)[act];
469 int mlevel=_node_num;
470 Value exc=(*_excess)[act];
472 for(OutEdgeIt e(_g,act);e!=INVALID; ++e) {
473 Node v = _g.target(e);
474 Value fc=(*_up)[e]-(*_flow)[e];
475 if(!_tol.positive(fc)) continue;
476 if((*_level)[v]<actlevel) {
477 if(!_tol.less(fc, exc)) {
478 _flow->set(e, (*_flow)[e] + exc);
479 _excess->set(v, (*_excess)[v] + exc);
480 if(!_level->active(v) && _tol.positive((*_excess)[v]))
483 _level->deactivate(act);
487 _flow->set(e, (*_up)[e]);
488 _excess->set(v, (*_excess)[v] + fc);
489 if(!_level->active(v) && _tol.positive((*_excess)[v]))
494 else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
496 for(InEdgeIt e(_g,act);e!=INVALID; ++e) {
497 Node v = _g.source(e);
498 Value fc=(*_flow)[e]-(*_lo)[e];
499 if(!_tol.positive(fc)) continue;
500 if((*_level)[v]<actlevel) {
501 if(!_tol.less(fc, exc)) {
502 _flow->set(e, (*_flow)[e] - exc);
503 _excess->set(v, (*_excess)[v] + exc);
504 if(!_level->active(v) && _tol.positive((*_excess)[v]))
507 _level->deactivate(act);
511 _flow->set(e, (*_lo)[e]);
512 _excess->set(v, (*_excess)[v] + fc);
513 if(!_level->active(v) && _tol.positive((*_excess)[v]))
518 else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
521 _excess->set(act, exc);
522 if(!_tol.positive(exc)) _level->deactivate(act);
523 else if(mlevel==_node_num) {
524 _level->liftHighestActiveToTop();
529 _level->liftHighestActive(mlevel+1);
530 if(_level->onLevel(actlevel)==0) {
541 /// Runs the circulation algorithm.
543 /// Runs the circulation algorithm.
544 /// \note fc.run() is just a shortcut of the following code.
547 /// return fc.start();
556 /// \name Query Functions
557 /// The result of the %Circulation algorithm can be obtained using
560 /// Before the use of these functions,
561 /// either run() or start() must be called.
567 ///Barrier is a set \e B of nodes for which
568 /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
569 ///holds. The existence of a set with this property prooves that a feasible
570 ///flow cannot exists.
571 ///\sa checkBarrier()
574 void barrierMap(GT &bar)
576 for(NodeIt n(_g);n!=INVALID;++n)
577 bar.set(n, (*_level)[n] >= _el);
580 ///Returns true if the node is in the barrier
582 ///Returns true if the node is in the barrier
584 bool barrier(const Node& node)
586 return (*_level)[node] >= _el;
589 /// \brief Returns the flow on the edge.
591 /// Sets the \c flowMap to the flow on the edges. This method can
592 /// be called after the second phase of algorithm.
593 Value flow(const Edge& edge) const {
594 return (*_flow)[edge];
599 /// \name Checker Functions
600 /// The feasibility of the results can be checked using
603 /// Before the use of these functions,
604 /// either run() or start() must be called.
608 ///Check if the \c flow is a feasible circulation
610 for(EdgeIt e(_g);e!=INVALID;++e)
611 if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
612 for(NodeIt n(_g);n!=INVALID;++n)
614 Value dif=-(*_delta)[n];
615 for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
616 for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
617 if(_tol.negative(dif)) return false;
622 ///Check whether or not the last execution provides a barrier
624 ///Check whether or not the last execution provides a barrier
629 for(NodeIt n(_g);n!=INVALID;++n)
632 for(EdgeIt e(_g);e!=INVALID;++e)
636 if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
637 else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
639 return _tol.negative(delta);