COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/circulation.h @ 2588:4d3bc1d04c1d

Last change on this file since 2588:4d3bc1d04c1d was 2553:bfced05fa852, checked in by Alpar Juttner, 16 years ago

Happy New Year to LEMON (+ better update-copyright-header script)

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-2008
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      // global relabeling tested, but in general case it provides
407      // worse performance for random graphs
408      _level->initStart();
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]))
414          _level->activate(n);
415    }
416
417    /// Initializes the internal data structures.
418
419    /// Initializes the internal data structures. This functions uses
420    /// greedy approach to construct the initial solution.
421    void greedyInit()
422    {
423      createStructures();
424     
425      for(NodeIt n(_g);n!=INVALID;++n) {
426        _excess->set(n, (*_delta)[n]);
427      }
428     
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]);
438        } else {
439          Value fc = -(*_excess)[_g.target(e)];
440          _flow->set(e, fc);
441          _excess->set(_g.target(e), 0);
442          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
443        }
444      }
445     
446      _level->initStart();
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]))
452          _level->activate(n);
453    }
454
455    ///Starts the algorithm
456
457    ///This function starts the algorithm.
458    ///\return This function returns true if it found a feasible circulation.
459    ///
460    ///\sa barrier()
461    bool start()
462    {
463     
464      Node act;
465      Node bact=INVALID;
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];
471
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]))
481                _level->activate(v);
482              _excess->set(act,0);
483              _level->deactivate(act);
484              goto next_l;
485            }
486            else {
487              _flow->set(e, (*_up)[e]);
488              _excess->set(v, (*_excess)[v] + fc);
489              if(!_level->active(v) && _tol.positive((*_excess)[v]))
490                _level->activate(v);
491              exc-=fc;
492            }
493          }
494          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
495        }
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]))
505                _level->activate(v);
506              _excess->set(act,0);
507              _level->deactivate(act);
508              goto next_l;
509            }
510            else {
511              _flow->set(e, (*_lo)[e]);
512              _excess->set(v, (*_excess)[v] + fc);
513              if(!_level->active(v) && _tol.positive((*_excess)[v]))
514                _level->activate(v);
515              exc-=fc;
516            }
517          }
518          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
519        }
520   
521        _excess->set(act, exc);
522        if(!_tol.positive(exc)) _level->deactivate(act);
523        else if(mlevel==_node_num) {
524          _level->liftHighestActiveToTop();
525          _el = _node_num;
526          return false;
527        }
528        else {
529          _level->liftHighestActive(mlevel+1);
530          if(_level->onLevel(actlevel)==0) {
531            _el = actlevel;
532            return false;
533          }
534        }
535      next_l:
536        ;
537      }
538      return true;
539    }
540
541    /// Runs the circulation algorithm. 
542
543    /// Runs the circulation algorithm.
544    /// \note fc.run() is just a shortcut of the following code.
545    /// \code
546    ///   fc.greedyInit();
547    ///   return fc.start();
548    /// \endcode
549    bool run() {
550      greedyInit();
551      return start();
552    }
553
554    /// @}
555
556    /// \name Query Functions
557    /// The result of the %Circulation algorithm can be obtained using
558    /// these functions.
559    /// \n
560    /// Before the use of these functions,
561    /// either run() or start() must be called.
562   
563    ///@{
564   
565    ///Returns a barrier
566   
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()
572    ///\sa run()
573    template<class GT>
574    void barrierMap(GT &bar)
575    {
576      for(NodeIt n(_g);n!=INVALID;++n)
577        bar.set(n, (*_level)[n] >= _el);
578    } 
579
580    ///Returns true if the node is in the barrier
581
582    ///Returns true if the node is in the barrier
583    ///\sa barrierMap()
584    bool barrier(const Node& node)
585    {
586      return (*_level)[node] >= _el;
587    } 
588
589    /// \brief Returns the flow on the edge.
590    ///
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];
595    }
596
597    /// @}
598
599    /// \name Checker Functions
600    /// The feasibility  of the results can be checked using
601    /// these functions.
602    /// \n
603    /// Before the use of these functions,
604    /// either run() or start() must be called.
605   
606    ///@{
607
608    ///Check if the  \c flow is a feasible circulation
609    bool checkFlow() {
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)
613        {
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;
618        }
619      return true;
620    }
621
622    ///Check whether or not the last execution provides a barrier
623
624    ///Check whether or not the last execution provides a barrier
625    ///\sa barrier()
626    bool checkBarrier()
627    {
628      Value delta=0;
629      for(NodeIt n(_g);n!=INVALID;++n)
630        if(barrier(n))
631          delta-=(*_delta)[n];
632      for(EdgeIt e(_g);e!=INVALID;++e)
633        {
634          Node s=_g.source(e);
635          Node t=_g.target(e);
636          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
637          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
638        }
639      return _tol.negative(delta);
640    }
641
642    /// @}
643
644  };
645 
646}
647
648#endif
Note: See TracBrowser for help on using the repository browser.