COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/circulation.h @ 2526:b7727edd44f2

Last change on this file since 2526:b7727edd44f2 was 2526:b7727edd44f2, checked in by Balazs Dezso, 12 years ago

Redesign Circulation interface according to new flow interface
New greedy approach initialization

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    /// The constructor of the class.
235
236    /// The constructor of the class.
237    /// \param g The directed graph the algorithm runs on.
238    /// \param lo The lower bound capacity of the edges.
239    /// \param up The upper bound capacity of the edges.
240    /// \param delta The lower bound on node excess.
241    Circulation(const Graph &g,const LCapMap &lo,
242                const UCapMap &up,const DeltaMap &delta)
243      : _g(g), _node_num(),
244        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
245        _level(0), _local_level(false), _excess(0), _el() {}
246
247    /// Destrcutor.
248    ~Circulation() {
249      destroyStructures();
250    }
251
252  private:
253
254    void createStructures() {
255      _node_num = _el = countNodes(_g);
256
257      if (!_flow) {
258        _flow = Traits::createFlowMap(_g);
259        _local_flow = true;
260      }
261      if (!_level) {
262        _level = Traits::createElevator(_g, _node_num);
263        _local_level = true;
264      }
265      if (!_excess) {
266        _excess = new ExcessMap(_g);
267      }
268    }
269
270    void destroyStructures() {
271      if (_local_flow) {
272        delete _flow;
273      }
274      if (_local_level) {
275        delete _level;
276      }
277      if (_excess) {
278        delete _excess;
279      }
280    }
281
282  public:
283
284    /// Sets the lower bound capacity map.
285
286    /// Sets the lower bound capacity map.
287    /// \return \c (*this)
288    Circulation& lowerCapMap(const LCapMap& map) {
289      _lo = &map;
290      return *this;
291    }
292
293    /// Sets the upper bound capacity map.
294
295    /// Sets the upper bound capacity map.
296    /// \return \c (*this)
297    Circulation& upperCapMap(const LCapMap& map) {
298      _up = &map;
299      return *this;
300    }
301
302    /// Sets the lower bound map on excess.
303
304    /// Sets the lower bound map on excess.
305    /// \return \c (*this)
306    Circulation& deltaMap(const DeltaMap& map) {
307      _delta = &map;
308      return *this;
309    }
310
311    /// Sets the flow map.
312
313    /// Sets the flow map.
314    /// \return \c (*this)
315    Circulation& flowMap(FlowMap& map) {
316      if (_local_flow) {
317        delete _flow;
318        _local_flow = false;
319      }
320      _flow = &map;
321      return *this;
322    }
323
324    /// Returns the flow map.
325
326    /// \return The flow map.
327    ///
328    const FlowMap& flowMap() {
329      return *_flow;
330    }
331
332    /// Sets the elevator.
333
334    /// Sets the elevator.
335    /// \return \c (*this)
336    Circulation& elevator(Elevator& elevator) {
337      if (_local_level) {
338        delete _level;
339        _local_level = false;
340      }
341      _level = &elevator;
342      return *this;
343    }
344
345    /// Returns the elevator.
346
347    /// \return The elevator.
348    ///
349    const Elevator& elevator() {
350      return *_level;
351    }
352   
353    /// Sets the tolerance used by algorithm.
354
355    /// Sets the tolerance used by algorithm.
356    ///
357    Circulation& tolerance(const Tolerance& tolerance) const {
358      _tol = tolerance;
359      return *this;
360    }
361
362    /// Returns the tolerance used by algorithm.
363
364    /// Returns the tolerance used by algorithm.
365    ///
366    const Tolerance& tolerance() const {
367      return tolerance;
368    }
369   
370    /// \name Execution control The simplest way to execute the
371    /// algorithm is to use one of the member functions called \c
372    /// run(). 
373    /// \n
374    /// If you need more control on initial solution or execution then
375    /// you have to call one \ref init() function and then the start()
376    /// function.
377   
378    ///@{
379 
380    /// Initializes the internal data structures.
381
382    /// Initializes the internal data structures. This function sets
383    /// all flow values to the lower bound. 
384    /// \return This function returns false if the initialization
385    /// process found a barrier.
386    void init()
387    {
388      createStructures();
389
390      for(NodeIt n(_g);n!=INVALID;++n) {
391        _excess->set(n, (*_delta)[n]);
392      }
393     
394      for (EdgeIt e(_g);e!=INVALID;++e) {
395        _flow->set(e, (*_lo)[e]);
396        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
397        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
398      }
399
400      typename Graph::template NodeMap<bool> reached(_g, false);
401
402
403      // global relabeling tested, but in general case it provides
404      // worse performance for random graphs
405      _level->initStart();
406      for(NodeIt n(_g);n!=INVALID;++n)
407        _level->initAddItem(n);
408      _level->initFinish();
409      for(NodeIt n(_g);n!=INVALID;++n)
410        if(_tol.positive((*_excess)[n]))
411          _level->activate(n);
412    }
413
414    /// Initializes the internal data structures.
415
416    /// Initializes the internal data structures. This functions uses
417    /// greedy approach to construct the initial solution.
418    void greedyInit()
419    {
420      createStructures();
421     
422      for(NodeIt n(_g);n!=INVALID;++n) {
423        _excess->set(n, (*_delta)[n]);
424      }
425     
426      for (EdgeIt e(_g);e!=INVALID;++e) {
427        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
428          _flow->set(e, (*_up)[e]);
429          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
430          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
431        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
432          _flow->set(e, (*_lo)[e]);
433          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
434          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
435        } else {
436          Value fc = -(*_excess)[_g.target(e)];
437          _flow->set(e, fc);
438          _excess->set(_g.target(e), 0);
439          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
440        }
441      }
442     
443      _level->initStart();
444      for(NodeIt n(_g);n!=INVALID;++n)
445        _level->initAddItem(n);
446      _level->initFinish();
447      for(NodeIt n(_g);n!=INVALID;++n)
448        if(_tol.positive((*_excess)[n]))
449          _level->activate(n);
450    }
451
452    ///Starts the algorithm
453
454    ///This function starts the algorithm.
455    ///\return This function returns true if it found a feasible circulation.
456    ///
457    ///\sa barrier()
458    bool start()
459    {
460     
461      Node act;
462      Node bact=INVALID;
463      Node last_activated=INVALID;
464      while((act=_level->highestActive())!=INVALID) {
465        int actlevel=(*_level)[act];
466        int mlevel=_node_num;
467        Value exc=(*_excess)[act];
468       
469        for(OutEdgeIt e(_g,act);e!=INVALID; ++e) {
470          Node v = _g.target(e);
471          Value fc=(*_up)[e]-(*_flow)[e];
472          if(!_tol.positive(fc)) continue;
473          if((*_level)[v]<actlevel) {
474            if(!_tol.less(fc, exc)) {
475              _flow->set(e, (*_flow)[e] + exc);
476              _excess->set(v, (*_excess)[v] + exc);
477              if(!_level->active(v) && _tol.positive((*_excess)[v]))
478                _level->activate(v);
479              _excess->set(act,0);
480              _level->deactivate(act);
481              goto next_l;
482            }
483            else {
484              _flow->set(e, (*_up)[e]);
485              _excess->set(v, (*_excess)[v] + exc);
486              if(!_level->active(v) && _tol.positive((*_excess)[v]))
487                _level->activate(v);
488              exc-=fc;
489            }
490          }
491          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
492        }
493        for(InEdgeIt e(_g,act);e!=INVALID; ++e) {
494          Node v = _g.source(e);
495          Value fc=(*_flow)[e]-(*_lo)[e];
496          if(!_tol.positive(fc)) continue;
497          if((*_level)[v]<actlevel) {
498            if(!_tol.less(fc, exc)) {
499              _flow->set(e, (*_flow)[e] - exc);
500              _excess->set(v, (*_excess)[v] + exc);
501              if(!_level->active(v) && _tol.positive((*_excess)[v]))
502                _level->activate(v);
503              _excess->set(act,0);
504              _level->deactivate(act);
505              goto next_l;
506            }
507            else {
508              _flow->set(e, (*_lo)[e]);
509              _excess->set(v, (*_excess)[v] + fc);
510              if(!_level->active(v) && _tol.positive((*_excess)[v]))
511                _level->activate(v);
512              exc-=fc;
513            }
514          }
515          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
516        }
517   
518        _excess->set(act, exc);
519        if(!_tol.positive(exc)) _level->deactivate(act);
520        else if(mlevel==_node_num) {
521          _level->liftHighestActiveToTop();
522          _el = _node_num;
523          return false;
524        }
525        else {
526          _level->liftHighestActive(mlevel+1);
527          if(_level->onLevel(actlevel)==0) {
528            _el = actlevel;
529            return false;
530          }
531        }
532      next_l:
533        ;
534      }
535      return true;
536    }
537
538    /// Runs the circulation algorithm. 
539
540    /// Runs the circulation algorithm.
541    /// \note fc.run() is just a shortcut of the following code.
542    /// \code
543    ///   fc.greedyInit();
544    ///   return fc.start();
545    /// \endcode
546    bool run() {
547      greedyInit();
548      return start();
549    }
550
551    /// @}
552
553    /// \name Query Functions
554    /// The result of the %Circulation algorithm can be obtained using
555    /// these functions.
556    /// \n
557    /// Before the use of these functions,
558    /// either run() or start() must be called.
559   
560    ///@{
561   
562    ///Returns a barrier
563   
564    ///Barrier is a set \e B of nodes for which
565    /// \f[ \sum_{v\in B}-delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
566    ///holds. The existence of a set with this property prooves that a feasible
567    ///flow cannot exists.
568    ///\sa checkBarrier()
569    ///\sa run()
570    template<class GT>
571    void barrierMap(GT &bar)
572    {
573      for(NodeIt n(_g);n!=INVALID;++n)
574        bar.set(n, (*_level)[n] >= _el);
575    } 
576
577    ///Returns true if the node is in the barrier
578
579    ///Returns true if the node is in the barrier
580    ///\sa barrierMap()
581    bool barrier(const Node& node)
582    {
583      return (*_level)[node] >= _el;
584    } 
585
586    /// \brief Returns the flow on the edge.
587    ///
588    /// Sets the \c flowMap to the flow on the edges. This method can
589    /// be called after the second phase of algorithm.
590    Value flow(const Edge& edge) const {
591      return (*_flow)[edge];
592    }
593
594    /// @}
595
596    /// \name Checker Functions
597    /// The feasibility  of the results can be checked using
598    /// these functions.
599    /// \n
600    /// Before the use of these functions,
601    /// either run() or start() must be called.
602   
603    ///@{
604
605    ///Check if the  \c flow is a feasible circulation
606    bool checkFlow() {
607      for(EdgeIt e(_g);e!=INVALID;++e)
608        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
609      for(NodeIt n(_g);n!=INVALID;++n)
610        {
611          Value dif=-(*_delta)[n];
612          for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
613          for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
614          if(_tol.negative(dif)) return false;
615        }
616      return true;
617    }
618
619    ///Check whether or not the last execution provides a barrier
620
621    ///Check whether or not the last execution provides a barrier
622    ///\sa barrier()
623    bool checkBarrier()
624    {
625      Value delta=0;
626      for(NodeIt n(_g);n!=INVALID;++n)
627        if(barrier(n))
628          delta-=(*_delta)[n];
629      for(EdgeIt e(_g);e!=INVALID;++e)
630        {
631          Node s=_g.source(e);
632          Node t=_g.target(e);
633          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
634          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
635        }
636      return _tol.negative(delta);
637    }
638
639    /// @}
640
641  };
642 
643}
644
645#endif
Note: See TracBrowser for help on using the repository browser.