COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/circulation.h @ 2533:aea952a1af99

Last change on this file since 2533:aea952a1af99 was 2533:aea952a1af99, checked in by Peter Kovacs, 16 years ago

Bug fixes.

File size: 17.4 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2007
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
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.
12 *
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
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_CIRCULATION_H
20#define LEMON_CIRCULATION_H
21
22#include <lemon/graph_utils.h>
23#include <iostream>
24#include <queue>
25#include <lemon/tolerance.h>
26#include <lemon/elevator.h>
27
28///\ingroup max_flow
29///\file
30///\brief Push-prelabel algorithm for finding a feasible circulation.
31///
32namespace lemon {
33
34  /// \brief Default traits class of Circulation class.
35  ///
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 {
42
43    /// \brief The graph type the algorithm runs on.
44    typedef _Graph Graph;
45
46    /// \brief The type of the map that stores the circulation lower
47    /// bound.
48    ///
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;
52
53    /// \brief The type of the map that stores the circulation upper
54    /// bound.
55    ///
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;
59
60    /// \brief The type of the map that stores the upper bound of
61    /// node excess.
62    ///
63    /// The type of the map that stores the lower bound of node
64    /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
65    /// concept.
66    typedef _DeltaMap DeltaMap;
67
68    /// \brief The type of the length of the edges.
69    typedef typename DeltaMap::Value Value;
70
71    /// \brief The map type that stores the flow values.
72    ///
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;
76
77    /// \brief Instantiates a FlowMap.
78    ///
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);
83    }
84
85    /// \brief The eleavator type used by Circulation algorithm.
86    ///
87    /// The elevator type used by Circulation algorithm.
88    ///
89    /// \sa Elevator
90    /// \sa LinkedElevator
91    typedef Elevator<Graph, typename Graph::Node> Elevator;
92   
93    /// \brief Instantiates an Elevator.
94    ///
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);
100    }
101
102    /// \brief The tolerance used by the algorithm
103    ///
104    /// The tolerance used by the algorithm to handle inexact computation.
105    typedef Tolerance<Value> Tolerance;
106
107  };
108 
109  ///Push-relabel algorithm for the Network Circulation Problem.
110 
111  ///\ingroup max_flow
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]
117  ///
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> >
125  class Circulation {
126
127    typedef _Traits Traits;
128    typedef typename Traits::Graph Graph;
129    GRAPH_TYPEDEFS(typename Graph);
130
131    typedef typename Traits::Value Value;
132
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;
139
140    typedef typename Graph::template NodeMap<Value> ExcessMap;
141
142    const Graph &_g;
143    int _node_num;
144
145    const LCapMap *_lo;
146    const UCapMap *_up;
147    const DeltaMap *_delta;
148
149    FlowMap *_flow;
150    bool _local_flow;
151
152    Elevator* _level;
153    bool _local_level;
154
155    ExcessMap* _excess;
156
157    Tolerance _tol;
158    int _el;
159
160  public:
161
162    typedef Circulation Create;
163
164    ///\name Named template parameters
165
166    ///@{
167
168    template <typename _FlowMap>
169    struct DefFlowMapTraits : public Traits {
170      typedef _FlowMap FlowMap;
171      static FlowMap *createFlowMap(const Graph&) {
172        throw UninitializedParameter();
173      }
174    };
175
176    /// \brief \ref named-templ-param "Named parameter" for setting
177    /// FlowMap type
178    ///
179    /// \ref named-templ-param "Named parameter" for setting FlowMap
180    /// type
181    template <typename _FlowMap>
182    struct DefFlowMap
183      : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
184                           DefFlowMapTraits<_FlowMap> > {
185      typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
186                          DefFlowMapTraits<_FlowMap> > Create;
187    };
188
189    template <typename _Elevator>
190    struct DefElevatorTraits : public Traits {
191      typedef _Elevator Elevator;
192      static Elevator *createElevator(const Graph&, int) {
193        throw UninitializedParameter();
194      }
195    };
196
197    /// \brief \ref named-templ-param "Named parameter" for setting
198    /// Elevator type
199    ///
200    /// \ref named-templ-param "Named parameter" for setting Elevator
201    /// type
202    template <typename _Elevator>
203    struct DefElevator
204      : public Circulation<Graph, LCapMap, UCapMap, DeltaMap,
205                           DefElevatorTraits<_Elevator> > {
206      typedef Circulation<Graph, LCapMap, UCapMap, DeltaMap,
207                          DefElevatorTraits<_Elevator> > Create;
208    };
209
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);
215      }
216    };
217
218    /// \brief \ref named-templ-param "Named parameter" for setting
219    /// Elevator type
220    ///
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;
230    };   
231
232    /// @}
233
234  protected:
235
236    Circulation() {}
237 
238  public:
239
240    /// The constructor of the class.
241
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() {}
252
253    /// Destrcutor.
254    ~Circulation() {
255      destroyStructures();
256    }
257
258  private:
259
260    void createStructures() {
261      _node_num = _el = countNodes(_g);
262
263      if (!_flow) {
264        _flow = Traits::createFlowMap(_g);
265        _local_flow = true;
266      }
267      if (!_level) {
268        _level = Traits::createElevator(_g, _node_num);
269        _local_level = true;
270      }
271      if (!_excess) {
272        _excess = new ExcessMap(_g);
273      }
274    }
275
276    void destroyStructures() {
277      if (_local_flow) {
278        delete _flow;
279      }
280      if (_local_level) {
281        delete _level;
282      }
283      if (_excess) {
284        delete _excess;
285      }
286    }
287
288  public:
289
290    /// Sets the lower bound capacity map.
291
292    /// Sets the lower bound capacity map.
293    /// \return \c (*this)
294    Circulation& lowerCapMap(const LCapMap& map) {
295      _lo = &map;
296      return *this;
297    }
298
299    /// Sets the upper bound capacity map.
300
301    /// Sets the upper bound capacity map.
302    /// \return \c (*this)
303    Circulation& upperCapMap(const LCapMap& map) {
304      _up = &map;
305      return *this;
306    }
307
308    /// Sets the lower bound map on excess.
309
310    /// Sets the lower bound map on excess.
311    /// \return \c (*this)
312    Circulation& deltaMap(const DeltaMap& map) {
313      _delta = &map;
314      return *this;
315    }
316
317    /// Sets the flow map.
318
319    /// Sets the flow map.
320    /// \return \c (*this)
321    Circulation& flowMap(FlowMap& map) {
322      if (_local_flow) {
323        delete _flow;
324        _local_flow = false;
325      }
326      _flow = &map;
327      return *this;
328    }
329
330    /// Returns the flow map.
331
332    /// \return The flow map.
333    ///
334    const FlowMap& flowMap() {
335      return *_flow;
336    }
337
338    /// Sets the elevator.
339
340    /// Sets the elevator.
341    /// \return \c (*this)
342    Circulation& elevator(Elevator& elevator) {
343      if (_local_level) {
344        delete _level;
345        _local_level = false;
346      }
347      _level = &elevator;
348      return *this;
349    }
350
351    /// Returns the elevator.
352
353    /// \return The elevator.
354    ///
355    const Elevator& elevator() {
356      return *_level;
357    }
358   
359    /// Sets the tolerance used by algorithm.
360
361    /// Sets the tolerance used by algorithm.
362    ///
363    Circulation& tolerance(const Tolerance& tolerance) const {
364      _tol = tolerance;
365      return *this;
366    }
367
368    /// Returns the tolerance used by algorithm.
369
370    /// Returns the tolerance used by algorithm.
371    ///
372    const Tolerance& tolerance() const {
373      return tolerance;
374    }
375   
376    /// \name Execution control
377    /// The simplest way to execute the algorithm is to use one of the
378    /// member functions called \c run().
379    /// \n
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()
382    /// function.
383   
384    ///@{
385 
386    /// Initializes the internal data structures.
387
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.
392    void init()
393    {
394      createStructures();
395
396      for(NodeIt n(_g);n!=INVALID;++n) {
397        _excess->set(n, (*_delta)[n]);
398      }
399     
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]);
404      }
405
406      typename Graph::template NodeMap<bool> reached(_g, false);
407
408
409      // global relabeling tested, but in general case it provides
410      // worse performance for random graphs
411      _level->initStart();
412      for(NodeIt n(_g);n!=INVALID;++n)
413        _level->initAddItem(n);
414      _level->initFinish();
415      for(NodeIt n(_g);n!=INVALID;++n)
416        if(_tol.positive((*_excess)[n]))
417          _level->activate(n);
418    }
419
420    /// Initializes the internal data structures.
421
422    /// Initializes the internal data structures. This functions uses
423    /// greedy approach to construct the initial solution.
424    void greedyInit()
425    {
426      createStructures();
427     
428      for(NodeIt n(_g);n!=INVALID;++n) {
429        _excess->set(n, (*_delta)[n]);
430      }
431     
432      for (EdgeIt e(_g);e!=INVALID;++e) {
433        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
434          _flow->set(e, (*_up)[e]);
435          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
436          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
437        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
438          _flow->set(e, (*_lo)[e]);
439          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
440          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
441        } else {
442          Value fc = -(*_excess)[_g.target(e)];
443          _flow->set(e, fc);
444          _excess->set(_g.target(e), 0);
445          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
446        }
447      }
448     
449      _level->initStart();
450      for(NodeIt n(_g);n!=INVALID;++n)
451        _level->initAddItem(n);
452      _level->initFinish();
453      for(NodeIt n(_g);n!=INVALID;++n)
454        if(_tol.positive((*_excess)[n]))
455          _level->activate(n);
456    }
457
458    ///Starts the algorithm
459
460    ///This function starts the algorithm.
461    ///\return This function returns true if it found a feasible circulation.
462    ///
463    ///\sa barrier()
464    bool start()
465    {
466     
467      Node act;
468      Node bact=INVALID;
469      Node last_activated=INVALID;
470      while((act=_level->highestActive())!=INVALID) {
471        int actlevel=(*_level)[act];
472        int mlevel=_node_num;
473        Value exc=(*_excess)[act];
474       
475        for(OutEdgeIt e(_g,act);e!=INVALID; ++e) {
476          Node v = _g.target(e);
477          Value fc=(*_up)[e]-(*_flow)[e];
478          if(!_tol.positive(fc)) continue;
479          if((*_level)[v]<actlevel) {
480            if(!_tol.less(fc, exc)) {
481              _flow->set(e, (*_flow)[e] + exc);
482              _excess->set(v, (*_excess)[v] + exc);
483              if(!_level->active(v) && _tol.positive((*_excess)[v]))
484                _level->activate(v);
485              _excess->set(act,0);
486              _level->deactivate(act);
487              goto next_l;
488            }
489            else {
490              _flow->set(e, (*_up)[e]);
491              _excess->set(v, (*_excess)[v] + exc);
492              if(!_level->active(v) && _tol.positive((*_excess)[v]))
493                _level->activate(v);
494              exc-=fc;
495            }
496          }
497          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
498        }
499        for(InEdgeIt e(_g,act);e!=INVALID; ++e) {
500          Node v = _g.source(e);
501          Value fc=(*_flow)[e]-(*_lo)[e];
502          if(!_tol.positive(fc)) continue;
503          if((*_level)[v]<actlevel) {
504            if(!_tol.less(fc, exc)) {
505              _flow->set(e, (*_flow)[e] - exc);
506              _excess->set(v, (*_excess)[v] + exc);
507              if(!_level->active(v) && _tol.positive((*_excess)[v]))
508                _level->activate(v);
509              _excess->set(act,0);
510              _level->deactivate(act);
511              goto next_l;
512            }
513            else {
514              _flow->set(e, (*_lo)[e]);
515              _excess->set(v, (*_excess)[v] + fc);
516              if(!_level->active(v) && _tol.positive((*_excess)[v]))
517                _level->activate(v);
518              exc-=fc;
519            }
520          }
521          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
522        }
523   
524        _excess->set(act, exc);
525        if(!_tol.positive(exc)) _level->deactivate(act);
526        else if(mlevel==_node_num) {
527          _level->liftHighestActiveToTop();
528          _el = _node_num;
529          return false;
530        }
531        else {
532          _level->liftHighestActive(mlevel+1);
533          if(_level->onLevel(actlevel)==0) {
534            _el = actlevel;
535            return false;
536          }
537        }
538      next_l:
539        ;
540      }
541      return true;
542    }
543
544    /// Runs the circulation algorithm. 
545
546    /// Runs the circulation algorithm.
547    /// \note fc.run() is just a shortcut of the following code.
548    /// \code
549    ///   fc.greedyInit();
550    ///   return fc.start();
551    /// \endcode
552    bool run() {
553      greedyInit();
554      return start();
555    }
556
557    /// @}
558
559    /// \name Query Functions
560    /// The result of the %Circulation algorithm can be obtained using
561    /// these functions.
562    /// \n
563    /// Before the use of these functions,
564    /// either run() or start() must be called.
565   
566    ///@{
567   
568    ///Returns a barrier
569   
570    ///Barrier is a set \e B of nodes for which
571    /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
572    ///holds. The existence of a set with this property prooves that a feasible
573    ///flow cannot exists.
574    ///\sa checkBarrier()
575    ///\sa run()
576    template<class GT>
577    void barrierMap(GT &bar)
578    {
579      for(NodeIt n(_g);n!=INVALID;++n)
580        bar.set(n, (*_level)[n] >= _el);
581    } 
582
583    ///Returns true if the node is in the barrier
584
585    ///Returns true if the node is in the barrier
586    ///\sa barrierMap()
587    bool barrier(const Node& node)
588    {
589      return (*_level)[node] >= _el;
590    } 
591
592    /// \brief Returns the flow on the edge.
593    ///
594    /// Sets the \c flowMap to the flow on the edges. This method can
595    /// be called after the second phase of algorithm.
596    Value flow(const Edge& edge) const {
597      return (*_flow)[edge];
598    }
599
600    /// @}
601
602    /// \name Checker Functions
603    /// The feasibility  of the results can be checked using
604    /// these functions.
605    /// \n
606    /// Before the use of these functions,
607    /// either run() or start() must be called.
608   
609    ///@{
610
611    ///Check if the  \c flow is a feasible circulation
612    bool checkFlow() {
613      for(EdgeIt e(_g);e!=INVALID;++e)
614        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
615      for(NodeIt n(_g);n!=INVALID;++n)
616        {
617          Value dif=-(*_delta)[n];
618          for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
619          for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
620          if(_tol.negative(dif)) return false;
621        }
622      return true;
623    }
624
625    ///Check whether or not the last execution provides a barrier
626
627    ///Check whether or not the last execution provides a barrier
628    ///\sa barrier()
629    bool checkBarrier()
630    {
631      Value delta=0;
632      for(NodeIt n(_g);n!=INVALID;++n)
633        if(barrier(n))
634          delta-=(*_delta)[n];
635      for(EdgeIt e(_g);e!=INVALID;++e)
636        {
637          Node s=_g.source(e);
638          Node t=_g.target(e);
639          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
640          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
641        }
642      return _tol.negative(delta);
643    }
644
645    /// @}
646
647  };
648 
649}
650
651#endif
Note: See TracBrowser for help on using the repository browser.