1.1 --- a/demo/circulation_demo.cc Wed Nov 28 17:40:41 2007 +0000
1.2 +++ b/demo/circulation_demo.cc Wed Nov 28 17:51:02 2007 +0000
1.3 @@ -52,7 +52,6 @@
1.4 Graph g;
1.5 EdgeMap lo(g);
1.6 EdgeMap up(g);
1.7 - EdgeMap x(g);
1.8 NodeMap delta(g);
1.9 NodeMap nid(g);
1.10 EdgeMap eid(g);
1.11 @@ -69,16 +68,16 @@
1.12 readNodeMap("coordinates_y", cy).
1.13 run();
1.14
1.15 - Circulation<Graph,int> gen(g,lo,up,delta,x);
1.16 - int ret=gen.run();
1.17 - if(ret==-1)
1.18 + Circulation<Graph> gen(g,lo,up,delta);
1.19 + bool ret=gen.run();
1.20 + if(ret)
1.21 {
1.22 std::cout << "\nA feasible flow has been found.\n";
1.23 - if(!gen.checkX(x)) std::cerr << "Oops!!!\n";
1.24 + if(!gen.checkFlow()) std::cerr << "Oops!!!\n";
1.25 GraphWriter<Graph>("circulation-output.lgf", g).
1.26 writeEdgeMap("lo_cap", lo).
1.27 writeEdgeMap("up_cap", up).
1.28 - writeEdgeMap("flow", x).
1.29 + writeEdgeMap("flow", gen.flowMap()).
1.30 writeNodeMap("delta", delta).
1.31 writeEdgeMap("label", eid).
1.32 writeNodeMap("label", nid).
1.33 @@ -89,8 +88,8 @@
1.34 else {
1.35 std::cout << "\nThere is no such a flow\n";
1.36 Graph::NodeMap<int> bar(g);
1.37 - gen.barrier(bar,ret);
1.38 - if(!gen.checkBarrier(bar)) std::cerr << "Dual Oops!!!\n";
1.39 + gen.barrierMap(bar);
1.40 + if(!gen.checkBarrier()) std::cerr << "Dual Oops!!!\n";
1.41
1.42 GraphWriter<Graph>("circulation-output.lgf", g).
1.43 writeEdgeMap("lo_cap", lo).
2.1 --- a/lemon/circulation.h Wed Nov 28 17:40:41 2007 +0000
2.2 +++ b/lemon/circulation.h Wed Nov 28 17:51:02 2007 +0000
2.3 @@ -31,268 +31,612 @@
2.4 ///
2.5 namespace lemon {
2.6
2.7 - ///Preflow algorithm for the Network Circulation Problem.
2.8 + /// \brief Default traits class of Circulation class.
2.9 + ///
2.10 + /// Default traits class of Circulation class.
2.11 + /// \param _Graph Graph type.
2.12 + /// \param _CapacityMap Type of capacity map.
2.13 + template <typename _Graph, typename _LCapMap,
2.14 + typename _UCapMap, typename _DeltaMap>
2.15 + struct CirculationDefaultTraits {
2.16 +
2.17 + /// \brief The graph type the algorithm runs on.
2.18 + typedef _Graph Graph;
2.19 +
2.20 + /// \brief The type of the map that stores the circulation lower
2.21 + /// bound.
2.22 + ///
2.23 + /// The type of the map that stores the circulation lower bound.
2.24 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
2.25 + typedef _LCapMap LCapMap;
2.26 +
2.27 + /// \brief The type of the map that stores the circulation upper
2.28 + /// bound.
2.29 + ///
2.30 + /// The type of the map that stores the circulation upper bound.
2.31 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
2.32 + typedef _UCapMap UCapMap;
2.33 +
2.34 + /// \brief The type of the map that stores the upper bound of
2.35 + /// node excess.
2.36 + ///
2.37 + /// The type of the map that stores the lower bound of node
2.38 + /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
2.39 + /// concept.
2.40 + typedef _DeltaMap DeltaMap;
2.41 +
2.42 + /// \brief The type of the length of the edges.
2.43 + typedef typename DeltaMap::Value Value;
2.44 +
2.45 + /// \brief The map type that stores the flow values.
2.46 + ///
2.47 + /// The map type that stores the flow values.
2.48 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
2.49 + typedef typename Graph::template EdgeMap<Value> FlowMap;
2.50 +
2.51 + /// \brief Instantiates a FlowMap.
2.52 + ///
2.53 + /// This function instantiates a \ref FlowMap.
2.54 + /// \param graph The graph, to which we would like to define the flow map.
2.55 + static FlowMap* createFlowMap(const Graph& graph) {
2.56 + return new FlowMap(graph);
2.57 + }
2.58 +
2.59 + /// \brief The eleavator type used by Circulation algorithm.
2.60 + ///
2.61 + /// The elevator type used by Circulation algorithm.
2.62 + ///
2.63 + /// \sa Elevator
2.64 + /// \sa LinkedElevator
2.65 + typedef Elevator<Graph, typename Graph::Node> Elevator;
2.66 +
2.67 + /// \brief Instantiates an Elevator.
2.68 + ///
2.69 + /// This function instantiates a \ref Elevator.
2.70 + /// \param graph The graph, to which we would like to define the elevator.
2.71 + /// \param max_level The maximum level of the elevator.
2.72 + static Elevator* createElevator(const Graph& graph, int max_level) {
2.73 + return new Elevator(graph, max_level);
2.74 + }
2.75 +
2.76 + /// \brief The tolerance used by the algorithm
2.77 + ///
2.78 + /// The tolerance used by the algorithm to handle inexact computation.
2.79 + typedef Tolerance<Value> Tolerance;
2.80 +
2.81 + };
2.82 +
2.83 + ///Push-relabel algorithm for the Network Circulation Problem.
2.84
2.85 ///\ingroup max_flow
2.86 - ///This class implements a preflow algorithm
2.87 + ///This class implements a push-relabel algorithm
2.88 ///for the Network Circulation Problem.
2.89 ///The exact formulation of this problem is the following.
2.90 /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq -delta(v)\quad \forall v\in V \f]
2.91 /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
2.92 ///
2.93 - template<class Graph,
2.94 - class Value,
2.95 - class FlowMap=typename Graph::template EdgeMap<Value>,
2.96 - class LCapMap=typename Graph::template EdgeMap<Value>,
2.97 - class UCapMap=LCapMap,
2.98 - class DeltaMap=typename Graph::template NodeMap<Value>
2.99 - >
2.100 + template<class _Graph,
2.101 + class _LCapMap=typename _Graph::template EdgeMap<int>,
2.102 + class _UCapMap=_LCapMap,
2.103 + class _DeltaMap=typename _Graph::template NodeMap<
2.104 + typename _UCapMap::Value>,
2.105 + class _Traits=CirculationDefaultTraits<_Graph, _LCapMap,
2.106 + _UCapMap, _DeltaMap> >
2.107 class Circulation {
2.108 - typedef typename Graph::Node Node;
2.109 - typedef typename Graph::NodeIt NodeIt;
2.110 - typedef typename Graph::Edge Edge;
2.111 - typedef typename Graph::EdgeIt EdgeIt;
2.112 - typedef typename Graph::InEdgeIt InEdgeIt;
2.113 - typedef typename Graph::OutEdgeIt OutEdgeIt;
2.114 -
2.115 +
2.116 + typedef _Traits Traits;
2.117 + typedef typename Traits::Graph Graph;
2.118 + GRAPH_TYPEDEFS(typename Graph);
2.119 +
2.120 + typedef typename Traits::Value Value;
2.121 +
2.122 + typedef typename Traits::LCapMap LCapMap;
2.123 + typedef typename Traits::UCapMap UCapMap;
2.124 + typedef typename Traits::DeltaMap DeltaMap;
2.125 + typedef typename Traits::FlowMap FlowMap;
2.126 + typedef typename Traits::Elevator Elevator;
2.127 + typedef typename Traits::Tolerance Tolerance;
2.128 +
2.129 + typedef typename Graph::template NodeMap<Value> ExcessMap;
2.130
2.131 const Graph &_g;
2.132 int _node_num;
2.133
2.134 - const LCapMap &_lo;
2.135 - const UCapMap &_up;
2.136 - const DeltaMap &_delta;
2.137 - FlowMap &_x;
2.138 - Tolerance<Value> _tol;
2.139 - Elevator<Graph,typename Graph::Node> _levels;
2.140 - typename Graph::template NodeMap<Value> _excess;
2.141 + const LCapMap *_lo;
2.142 + const UCapMap *_up;
2.143 + const DeltaMap *_delta;
2.144 +
2.145 + FlowMap *_flow;
2.146 + bool _local_flow;
2.147 +
2.148 + Elevator* _level;
2.149 + bool _local_level;
2.150 +
2.151 + ExcessMap* _excess;
2.152 +
2.153 + Tolerance _tol;
2.154 + int _el;
2.155
2.156 public:
2.157 - ///\e
2.158 - Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up,
2.159 - const DeltaMap &delta,FlowMap &x) :
2.160 - _g(g),
2.161 - _node_num(countNodes(g)),
2.162 - _lo(lo),_up(up),_delta(delta),_x(x),
2.163 - _levels(g,_node_num),
2.164 - _excess(g)
2.165 - {
2.166 +
2.167 + typedef Circulation Create;
2.168 +
2.169 + ///\name Named template parameters
2.170 +
2.171 + ///@{
2.172 +
2.173 + template <typename _FlowMap>
2.174 + struct DefFlowMapTraits : public Traits {
2.175 + typedef _FlowMap FlowMap;
2.176 + static FlowMap *createFlowMap(const Graph&) {
2.177 + throw UninitializedParameter();
2.178 + }
2.179 + };
2.180 +
2.181 + /// \brief \ref named-templ-param "Named parameter" for setting
2.182 + /// FlowMap type
2.183 + ///
2.184 + /// \ref named-templ-param "Named parameter" for setting FlowMap
2.185 + /// type
2.186 + template <typename _FlowMap>
2.187 + struct DefFlowMap
2.188 + : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
2.189 + DefFlowMapTraits<_FlowMap> > {
2.190 + typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
2.191 + DefFlowMapTraits<_FlowMap> > Create;
2.192 + };
2.193 +
2.194 + template <typename _Elevator>
2.195 + struct DefElevatorTraits : public Traits {
2.196 + typedef _Elevator Elevator;
2.197 + static Elevator *createElevator(const Graph&, int) {
2.198 + throw UninitializedParameter();
2.199 + }
2.200 + };
2.201 +
2.202 + /// \brief \ref named-templ-param "Named parameter" for setting
2.203 + /// Elevator type
2.204 + ///
2.205 + /// \ref named-templ-param "Named parameter" for setting Elevator
2.206 + /// type
2.207 + template <typename _Elevator>
2.208 + struct DefElevator
2.209 + : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
2.210 + DefElevatorTraits<_Elevator> > {
2.211 + typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
2.212 + DefElevatorTraits<_Elevator> > Create;
2.213 + };
2.214 +
2.215 + template <typename _Elevator>
2.216 + struct DefStandardElevatorTraits : public Traits {
2.217 + typedef _Elevator Elevator;
2.218 + static Elevator *createElevator(const Graph& graph, int max_level) {
2.219 + return new Elevator(graph, max_level);
2.220 + }
2.221 + };
2.222 +
2.223 + /// \brief \ref named-templ-param "Named parameter" for setting
2.224 + /// Elevator type
2.225 + ///
2.226 + /// \ref named-templ-param "Named parameter" for setting Elevator
2.227 + /// type. The Elevator should be standard constructor interface, ie.
2.228 + /// the graph and the maximum level should be passed to it.
2.229 + template <typename _Elevator>
2.230 + struct DefStandardElevator
2.231 + : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
2.232 + DefStandardElevatorTraits<_Elevator> > {
2.233 + typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
2.234 + DefStandardElevatorTraits<_Elevator> > Create;
2.235 + };
2.236 +
2.237 + /// @}
2.238 +
2.239 + /// The constructor of the class.
2.240 +
2.241 + /// The constructor of the class.
2.242 + /// \param g The directed graph the algorithm runs on.
2.243 + /// \param lo The lower bound capacity of the edges.
2.244 + /// \param up The upper bound capacity of the edges.
2.245 + /// \param delta The lower bound on node excess.
2.246 + Circulation(const Graph &g,const LCapMap &lo,
2.247 + const UCapMap &up,const DeltaMap &delta)
2.248 + : _g(g), _node_num(),
2.249 + _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
2.250 + _level(0), _local_level(false), _excess(0), _el() {}
2.251 +
2.252 + /// Destrcutor.
2.253 + ~Circulation() {
2.254 + destroyStructures();
2.255 }
2.256 -
2.257 +
2.258 private:
2.259
2.260 - void addExcess(Node n,Value v)
2.261 - {
2.262 - if(_tol.positive(_excess[n]+=v))
2.263 - {
2.264 - if(!_levels.active(n)) _levels.activate(n);
2.265 - }
2.266 - else if(_levels.active(n)) _levels.deactivate(n);
2.267 + void createStructures() {
2.268 + _node_num = _el = countNodes(_g);
2.269 +
2.270 + if (!_flow) {
2.271 + _flow = Traits::createFlowMap(_g);
2.272 + _local_flow = true;
2.273 + }
2.274 + if (!_level) {
2.275 + _level = Traits::createElevator(_g, _node_num);
2.276 + _local_level = true;
2.277 + }
2.278 + if (!_excess) {
2.279 + _excess = new ExcessMap(_g);
2.280 + }
2.281 }
2.282 -
2.283 - void init()
2.284 - {
2.285 -
2.286 - _x=_lo;
2.287
2.288 - for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=_delta[n];
2.289 -
2.290 - for(EdgeIt e(_g);e!=INVALID;++e)
2.291 - {
2.292 - _excess[_g.target(e)]+=_x[e];
2.293 - _excess[_g.source(e)]-=_x[e];
2.294 - }
2.295 -
2.296 - _levels.initStart();
2.297 - for(NodeIt n(_g);n!=INVALID;++n)
2.298 - _levels.initAddItem(n);
2.299 - _levels.initFinish();
2.300 - for(NodeIt n(_g);n!=INVALID;++n)
2.301 - if(_tol.positive(_excess[n]))
2.302 - _levels.activate(n);
2.303 + void destroyStructures() {
2.304 + if (_local_flow) {
2.305 + delete _flow;
2.306 + }
2.307 + if (_local_level) {
2.308 + delete _level;
2.309 + }
2.310 + if (_excess) {
2.311 + delete _excess;
2.312 + }
2.313 }
2.314
2.315 public:
2.316 - ///Check if \c x is a feasible circulation
2.317 - template<class FT>
2.318 - bool checkX(FT &x) {
2.319 +
2.320 + /// Sets the lower bound capacity map.
2.321 +
2.322 + /// Sets the lower bound capacity map.
2.323 + /// \return \c (*this)
2.324 + Circulation& lowerCapMap(const LCapMap& map) {
2.325 + _lo = ↦
2.326 + return *this;
2.327 + }
2.328 +
2.329 + /// Sets the upper bound capacity map.
2.330 +
2.331 + /// Sets the upper bound capacity map.
2.332 + /// \return \c (*this)
2.333 + Circulation& upperCapMap(const LCapMap& map) {
2.334 + _up = ↦
2.335 + return *this;
2.336 + }
2.337 +
2.338 + /// Sets the lower bound map on excess.
2.339 +
2.340 + /// Sets the lower bound map on excess.
2.341 + /// \return \c (*this)
2.342 + Circulation& deltaMap(const DeltaMap& map) {
2.343 + _delta = ↦
2.344 + return *this;
2.345 + }
2.346 +
2.347 + /// Sets the flow map.
2.348 +
2.349 + /// Sets the flow map.
2.350 + /// \return \c (*this)
2.351 + Circulation& flowMap(FlowMap& map) {
2.352 + if (_local_flow) {
2.353 + delete _flow;
2.354 + _local_flow = false;
2.355 + }
2.356 + _flow = ↦
2.357 + return *this;
2.358 + }
2.359 +
2.360 + /// Returns the flow map.
2.361 +
2.362 + /// \return The flow map.
2.363 + ///
2.364 + const FlowMap& flowMap() {
2.365 + return *_flow;
2.366 + }
2.367 +
2.368 + /// Sets the elevator.
2.369 +
2.370 + /// Sets the elevator.
2.371 + /// \return \c (*this)
2.372 + Circulation& elevator(Elevator& elevator) {
2.373 + if (_local_level) {
2.374 + delete _level;
2.375 + _local_level = false;
2.376 + }
2.377 + _level = &elevator;
2.378 + return *this;
2.379 + }
2.380 +
2.381 + /// Returns the elevator.
2.382 +
2.383 + /// \return The elevator.
2.384 + ///
2.385 + const Elevator& elevator() {
2.386 + return *_level;
2.387 + }
2.388 +
2.389 + /// Sets the tolerance used by algorithm.
2.390 +
2.391 + /// Sets the tolerance used by algorithm.
2.392 + ///
2.393 + Circulation& tolerance(const Tolerance& tolerance) const {
2.394 + _tol = tolerance;
2.395 + return *this;
2.396 + }
2.397 +
2.398 + /// Returns the tolerance used by algorithm.
2.399 +
2.400 + /// Returns the tolerance used by algorithm.
2.401 + ///
2.402 + const Tolerance& tolerance() const {
2.403 + return tolerance;
2.404 + }
2.405 +
2.406 + /// \name Execution control The simplest way to execute the
2.407 + /// algorithm is to use one of the member functions called \c
2.408 + /// run().
2.409 + /// \n
2.410 + /// If you need more control on initial solution or execution then
2.411 + /// you have to call one \ref init() function and then the start()
2.412 + /// function.
2.413 +
2.414 + ///@{
2.415 +
2.416 + /// Initializes the internal data structures.
2.417 +
2.418 + /// Initializes the internal data structures. This function sets
2.419 + /// all flow values to the lower bound.
2.420 + /// \return This function returns false if the initialization
2.421 + /// process found a barrier.
2.422 + void init()
2.423 + {
2.424 + createStructures();
2.425 +
2.426 + for(NodeIt n(_g);n!=INVALID;++n) {
2.427 + _excess->set(n, (*_delta)[n]);
2.428 + }
2.429 +
2.430 + for (EdgeIt e(_g);e!=INVALID;++e) {
2.431 + _flow->set(e, (*_lo)[e]);
2.432 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
2.433 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
2.434 + }
2.435 +
2.436 + typename Graph::template NodeMap<bool> reached(_g, false);
2.437 +
2.438 +
2.439 + // global relabeling tested, but in general case it provides
2.440 + // worse performance for random graphs
2.441 + _level->initStart();
2.442 + for(NodeIt n(_g);n!=INVALID;++n)
2.443 + _level->initAddItem(n);
2.444 + _level->initFinish();
2.445 + for(NodeIt n(_g);n!=INVALID;++n)
2.446 + if(_tol.positive((*_excess)[n]))
2.447 + _level->activate(n);
2.448 + }
2.449 +
2.450 + /// Initializes the internal data structures.
2.451 +
2.452 + /// Initializes the internal data structures. This functions uses
2.453 + /// greedy approach to construct the initial solution.
2.454 + void greedyInit()
2.455 + {
2.456 + createStructures();
2.457 +
2.458 + for(NodeIt n(_g);n!=INVALID;++n) {
2.459 + _excess->set(n, (*_delta)[n]);
2.460 + }
2.461 +
2.462 + for (EdgeIt e(_g);e!=INVALID;++e) {
2.463 + if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
2.464 + _flow->set(e, (*_up)[e]);
2.465 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
2.466 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
2.467 + } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
2.468 + _flow->set(e, (*_lo)[e]);
2.469 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
2.470 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
2.471 + } else {
2.472 + Value fc = -(*_excess)[_g.target(e)];
2.473 + _flow->set(e, fc);
2.474 + _excess->set(_g.target(e), 0);
2.475 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
2.476 + }
2.477 + }
2.478 +
2.479 + _level->initStart();
2.480 + for(NodeIt n(_g);n!=INVALID;++n)
2.481 + _level->initAddItem(n);
2.482 + _level->initFinish();
2.483 + for(NodeIt n(_g);n!=INVALID;++n)
2.484 + if(_tol.positive((*_excess)[n]))
2.485 + _level->activate(n);
2.486 + }
2.487 +
2.488 + ///Starts the algorithm
2.489 +
2.490 + ///This function starts the algorithm.
2.491 + ///\return This function returns true if it found a feasible circulation.
2.492 + ///
2.493 + ///\sa barrier()
2.494 + bool start()
2.495 + {
2.496 +
2.497 + Node act;
2.498 + Node bact=INVALID;
2.499 + Node last_activated=INVALID;
2.500 + while((act=_level->highestActive())!=INVALID) {
2.501 + int actlevel=(*_level)[act];
2.502 + int mlevel=_node_num;
2.503 + Value exc=(*_excess)[act];
2.504 +
2.505 + for(OutEdgeIt e(_g,act);e!=INVALID; ++e) {
2.506 + Node v = _g.target(e);
2.507 + Value fc=(*_up)[e]-(*_flow)[e];
2.508 + if(!_tol.positive(fc)) continue;
2.509 + if((*_level)[v]<actlevel) {
2.510 + if(!_tol.less(fc, exc)) {
2.511 + _flow->set(e, (*_flow)[e] + exc);
2.512 + _excess->set(v, (*_excess)[v] + exc);
2.513 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
2.514 + _level->activate(v);
2.515 + _excess->set(act,0);
2.516 + _level->deactivate(act);
2.517 + goto next_l;
2.518 + }
2.519 + else {
2.520 + _flow->set(e, (*_up)[e]);
2.521 + _excess->set(v, (*_excess)[v] + exc);
2.522 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
2.523 + _level->activate(v);
2.524 + exc-=fc;
2.525 + }
2.526 + }
2.527 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
2.528 + }
2.529 + for(InEdgeIt e(_g,act);e!=INVALID; ++e) {
2.530 + Node v = _g.source(e);
2.531 + Value fc=(*_flow)[e]-(*_lo)[e];
2.532 + if(!_tol.positive(fc)) continue;
2.533 + if((*_level)[v]<actlevel) {
2.534 + if(!_tol.less(fc, exc)) {
2.535 + _flow->set(e, (*_flow)[e] - exc);
2.536 + _excess->set(v, (*_excess)[v] + exc);
2.537 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
2.538 + _level->activate(v);
2.539 + _excess->set(act,0);
2.540 + _level->deactivate(act);
2.541 + goto next_l;
2.542 + }
2.543 + else {
2.544 + _flow->set(e, (*_lo)[e]);
2.545 + _excess->set(v, (*_excess)[v] + fc);
2.546 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
2.547 + _level->activate(v);
2.548 + exc-=fc;
2.549 + }
2.550 + }
2.551 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
2.552 + }
2.553 +
2.554 + _excess->set(act, exc);
2.555 + if(!_tol.positive(exc)) _level->deactivate(act);
2.556 + else if(mlevel==_node_num) {
2.557 + _level->liftHighestActiveToTop();
2.558 + _el = _node_num;
2.559 + return false;
2.560 + }
2.561 + else {
2.562 + _level->liftHighestActive(mlevel+1);
2.563 + if(_level->onLevel(actlevel)==0) {
2.564 + _el = actlevel;
2.565 + return false;
2.566 + }
2.567 + }
2.568 + next_l:
2.569 + ;
2.570 + }
2.571 + return true;
2.572 + }
2.573 +
2.574 + /// Runs the circulation algorithm.
2.575 +
2.576 + /// Runs the circulation algorithm.
2.577 + /// \note fc.run() is just a shortcut of the following code.
2.578 + /// \code
2.579 + /// fc.greedyInit();
2.580 + /// return fc.start();
2.581 + /// \endcode
2.582 + bool run() {
2.583 + greedyInit();
2.584 + return start();
2.585 + }
2.586 +
2.587 + /// @}
2.588 +
2.589 + /// \name Query Functions
2.590 + /// The result of the %Circulation algorithm can be obtained using
2.591 + /// these functions.
2.592 + /// \n
2.593 + /// Before the use of these functions,
2.594 + /// either run() or start() must be called.
2.595 +
2.596 + ///@{
2.597 +
2.598 + ///Returns a barrier
2.599 +
2.600 + ///Barrier is a set \e B of nodes for which
2.601 + /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
2.602 + ///holds. The existence of a set with this property prooves that a feasible
2.603 + ///flow cannot exists.
2.604 + ///\sa checkBarrier()
2.605 + ///\sa run()
2.606 + template<class GT>
2.607 + void barrierMap(GT &bar)
2.608 + {
2.609 + for(NodeIt n(_g);n!=INVALID;++n)
2.610 + bar.set(n, (*_level)[n] >= _el);
2.611 + }
2.612 +
2.613 + ///Returns true if the node is in the barrier
2.614 +
2.615 + ///Returns true if the node is in the barrier
2.616 + ///\sa barrierMap()
2.617 + bool barrier(const Node& node)
2.618 + {
2.619 + return (*_level)[node] >= _el;
2.620 + }
2.621 +
2.622 + /// \brief Returns the flow on the edge.
2.623 + ///
2.624 + /// Sets the \c flowMap to the flow on the edges. This method can
2.625 + /// be called after the second phase of algorithm.
2.626 + Value flow(const Edge& edge) const {
2.627 + return (*_flow)[edge];
2.628 + }
2.629 +
2.630 + /// @}
2.631 +
2.632 + /// \name Checker Functions
2.633 + /// The feasibility of the results can be checked using
2.634 + /// these functions.
2.635 + /// \n
2.636 + /// Before the use of these functions,
2.637 + /// either run() or start() must be called.
2.638 +
2.639 + ///@{
2.640 +
2.641 + ///Check if the \c flow is a feasible circulation
2.642 + bool checkFlow() {
2.643 for(EdgeIt e(_g);e!=INVALID;++e)
2.644 - if(x[e]<_lo[e]||x[e]>_up[e]) return false;
2.645 + if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
2.646 for(NodeIt n(_g);n!=INVALID;++n)
2.647 {
2.648 - Value dif=-_delta[n];
2.649 - for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=x[e];
2.650 - for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=x[e];
2.651 + Value dif=-(*_delta)[n];
2.652 + for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
2.653 + for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
2.654 if(_tol.negative(dif)) return false;
2.655 }
2.656 return true;
2.657 - };
2.658 - ///Check if the default \c x is a feasible circulation
2.659 - bool checkX() { return checkX(_x); }
2.660 + }
2.661
2.662 - ///Check if \c bar is a real barrier
2.663 -
2.664 - ///Check if \c bar is a real barrier
2.665 - ///\sa barrier()
2.666 - template<class GT>
2.667 - bool checkBarrier(GT &bar)
2.668 - {
2.669 - Value delta=0;
2.670 - for(NodeIt n(_g);n!=INVALID;++n)
2.671 - if(bar[n])
2.672 - delta-=_delta[n];
2.673 - for(EdgeIt e(_g);e!=INVALID;++e)
2.674 - {
2.675 - Node s=_g.source(e);
2.676 - Node t=_g.target(e);
2.677 - if(bar[s]&&!bar[t]) delta+=_up[e];
2.678 - else if(bar[t]&&!bar[s]) delta-=_lo[e];
2.679 - }
2.680 - return _tol.negative(delta);
2.681 - }
2.682 ///Check whether or not the last execution provides a barrier
2.683
2.684 ///Check whether or not the last execution provides a barrier
2.685 ///\sa barrier()
2.686 bool checkBarrier()
2.687 {
2.688 - typename Graph:: template NodeMap<bool> bar(_g);
2.689 - barrier(bar);
2.690 - return checkBarrier(bar);
2.691 + Value delta=0;
2.692 + for(NodeIt n(_g);n!=INVALID;++n)
2.693 + if(barrier(n))
2.694 + delta-=(*_delta)[n];
2.695 + for(EdgeIt e(_g);e!=INVALID;++e)
2.696 + {
2.697 + Node s=_g.source(e);
2.698 + Node t=_g.target(e);
2.699 + if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
2.700 + else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
2.701 + }
2.702 + return _tol.negative(delta);
2.703 }
2.704 - ///Run the algorithm
2.705
2.706 - ///This function runs the algorithm.
2.707 - ///\return This function returns -1 if it found a feasible circulation.
2.708 - /// nonnegative values (including 0) mean that no feasible solution is
2.709 - /// found. In this case the return value means an "empty level".
2.710 - ///
2.711 - ///\sa barrier()
2.712 - int run()
2.713 - {
2.714 - init();
2.715 -
2.716 -#ifdef LEMON_CIRCULATION_DEBUG
2.717 - for(NodeIt n(_g);n!=INVALID;++n)
2.718 - std::cerr<< _levels[n] << ' ';
2.719 - std::cerr << std::endl;
2.720 -#endif
2.721 - Node act;
2.722 - Node bact=INVALID;
2.723 - Node last_activated=INVALID;
2.724 - while((act=_levels.highestActive())!=INVALID) {
2.725 - int actlevel=_levels[act];
2.726 - int tlevel;
2.727 - int mlevel=_node_num;
2.728 - Value exc=_excess[act];
2.729 - Value fc;
2.730 -
2.731 -#ifdef LEMON_CIRCULATION_DEBUG
2.732 - for(NodeIt n(_g);n!=INVALID;++n)
2.733 - std::cerr<< _levels[n] << ' ';
2.734 - std::cerr << std::endl;
2.735 - std::cerr << "Process node " << _g.id(act)
2.736 - << " on level " << actlevel
2.737 - << " with excess " << exc
2.738 - << std::endl;
2.739 -#endif
2.740 - for(OutEdgeIt e(_g,act);e!=INVALID; ++e)
2.741 - if((fc=_up[e]-_x[e])>0)
2.742 - if((tlevel=_levels[_g.target(e)])<actlevel)
2.743 - if(fc<=exc) {
2.744 - _x[e]=_up[e];
2.745 - addExcess(_g.target(e),fc);
2.746 - exc-=fc;
2.747 -#ifdef LEMON_CIRCULATION_DEBUG
2.748 - std::cerr << " Push " << fc
2.749 - << " toward " << _g.id(_g.target(e)) << std::endl;
2.750 -#endif
2.751 - }
2.752 - else {
2.753 - _x[e]+=exc;
2.754 - addExcess(_g.target(e),exc);
2.755 - //exc=0;
2.756 - _excess[act]=0;
2.757 - _levels.deactivate(act);
2.758 -#ifdef LEMON_CIRCULATION_DEBUG
2.759 - std::cerr << " Push " << exc
2.760 - << " toward " << _g.id(_g.target(e)) << std::endl;
2.761 - std::cerr << " Deactivate\n";
2.762 -#endif
2.763 - goto next_l;
2.764 - }
2.765 - else if(tlevel<mlevel) mlevel=tlevel;
2.766 -
2.767 - for(InEdgeIt e(_g,act);e!=INVALID; ++e)
2.768 - if((fc=_x[e]-_lo[e])>0)
2.769 - if((tlevel=_levels[_g.source(e)])<actlevel)
2.770 - if(fc<=exc) {
2.771 - _x[e]=_lo[e];
2.772 - addExcess(_g.source(e),fc);
2.773 - exc-=fc;
2.774 -#ifdef LEMON_CIRCULATION_DEBUG
2.775 - std::cerr << " Push " << fc
2.776 - << " toward " << _g.id(_g.source(e)) << std::endl;
2.777 -#endif
2.778 - }
2.779 - else {
2.780 - _x[e]-=exc;
2.781 - addExcess(_g.source(e),exc);
2.782 - //exc=0;
2.783 - _excess[act]=0;
2.784 - _levels.deactivate(act);
2.785 -#ifdef LEMON_CIRCULATION_DEBUG
2.786 - std::cerr << " Push " << exc
2.787 - << " toward " << _g.id(_g.source(e)) << std::endl;
2.788 - std::cerr << " Deactivate\n";
2.789 -#endif
2.790 - goto next_l;
2.791 - }
2.792 - else if(tlevel<mlevel) mlevel=tlevel;
2.793 -
2.794 - _excess[act]=exc;
2.795 - if(!_tol.positive(exc)) _levels.deactivate(act);
2.796 - else if(mlevel==_node_num) {
2.797 - _levels.liftHighestActiveToTop();
2.798 -#ifdef LEMON_CIRCULATION_DEBUG
2.799 - std::cerr << " Lift to level " << _node_num << std::endl;
2.800 -#endif
2.801 - return _levels.onLevel(_node_num-1)==0?_node_num-1:actlevel;
2.802 - }
2.803 - else {
2.804 - _levels.liftHighestActive(mlevel+1);
2.805 -#ifdef LEMON_CIRCULATION_DEBUG
2.806 - std::cerr << " Lift to level " << mlevel+1 << std::endl;
2.807 -#endif
2.808 - if(_levels.onLevel(actlevel)==0)
2.809 - return actlevel;
2.810 - }
2.811 - next_l:
2.812 - ;
2.813 - }
2.814 -#ifdef LEMON_CIRCULATION_DEBUG
2.815 - std::cerr << "Feasible flow found.\n";
2.816 -#endif
2.817 - return -1;
2.818 - }
2.819 -
2.820 - ///Return a barrier
2.821 -
2.822 - ///Barrier is a set \e B of nodes for which
2.823 - /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
2.824 - ///holds. The existence of a set with this property prooves that a feasible
2.825 - ///flow cannot exists.
2.826 - ///\pre The run() must have been executed, and its return value was -1.
2.827 - ///\sa checkBarrier()
2.828 - ///\sa run()
2.829 - template<class GT>
2.830 - void barrier(GT &bar,int empty_level=-1)
2.831 - {
2.832 - if(empty_level==-1)
2.833 - for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ;
2.834 - for(NodeIt n(_g);n!=INVALID;++n)
2.835 - bar[n] = _levels[n]>empty_level;
2.836 - }
2.837 + /// @}
2.838
2.839 };
2.840
3.1 --- a/lemon/cycle_canceling.h Wed Nov 28 17:40:41 2007 +0000
3.2 +++ b/lemon/cycle_canceling.h Wed Nov 28 17:51:02 2007 +0000
3.3 @@ -323,11 +323,12 @@
3.4 if (sum != 0) return false;
3.5
3.6 // Finding a feasible flow
3.7 - Circulation< Graph, Capacity, FlowMap, ConstMap<Edge, Capacity>,
3.8 - CapacityRefMap, SupplyMap >
3.9 + Circulation< Graph, Capacity, ConstMap<Edge, Capacity>,
3.10 + CapacityRefMap, SupplyMap >::DefFlowMap<FlowMap>::Create
3.11 circulation( graph, constMap<Edge>((Capacity)0),
3.12 - capacity, supply, flow );
3.13 - return circulation.run() == -1;
3.14 + capacity, supply);
3.15 + circulation.flowMap(flowMap);
3.16 + return circulation.run();
3.17 }
3.18
3.19 #ifdef LIMITED_CYCLE_CANCELING