COIN-OR::LEMON - Graph Library

source: lemon-main/lemon/circulation.h @ 399:5a7dbeaed70e

Last change on this file since 399:5a7dbeaed70e was 399:5a7dbeaed70e, checked in by Alpar Juttner <alpar@…>, 15 years ago

Port Circulation from svn -r3516 (#175)
Namely,

  • port the files
  • apply the migrate script
  • apply the unify script
  • fix the compilation
  • strip the demo input file
  • break long lines
File size: 18.6 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
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 <iostream>
23#include <queue>
24#include <lemon/tolerance.h>
25#include <lemon/elevator.h>
26
27///\ingroup max_flow
28///\file
29///\brief Push-prelabel algorithm for finding a feasible circulation.
30///
31namespace lemon {
32
33  /// \brief Default traits class of Circulation class.
34  ///
35  /// Default traits class of Circulation class.
36  /// \param _Graph Digraph type.
37  /// \param _CapacityMap Type of capacity map.
38  template <typename _Graph, typename _LCapMap,
39            typename _UCapMap, typename _DeltaMap>
40  struct CirculationDefaultTraits {
41
42    /// \brief The digraph type the algorithm runs on.
43    typedef _Graph Digraph;
44
45    /// \brief The type of the map that stores the circulation lower
46    /// bound.
47    ///
48    /// The type of the map that stores the circulation lower bound.
49    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
50    typedef _LCapMap LCapMap;
51
52    /// \brief The type of the map that stores the circulation upper
53    /// bound.
54    ///
55    /// The type of the map that stores the circulation upper bound.
56    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
57    typedef _UCapMap UCapMap;
58
59    /// \brief The type of the map that stores the upper bound of
60    /// node excess.
61    ///
62    /// The type of the map that stores the lower bound of node
63    /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
64    /// concept.
65    typedef _DeltaMap DeltaMap;
66
67    /// \brief The type of the length of the arcs.
68    typedef typename DeltaMap::Value Value;
69
70    /// \brief The map type that stores the flow values.
71    ///
72    /// The map type that stores the flow values.
73    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
74    typedef typename Digraph::template ArcMap<Value> FlowMap;
75
76    /// \brief Instantiates a FlowMap.
77    ///
78    /// This function instantiates a \ref FlowMap.
79    /// \param digraph The digraph, to which we would like to define
80    /// the flow map.
81    static FlowMap* createFlowMap(const Digraph& digraph) {
82      return new FlowMap(digraph);
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 lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
92
93    /// \brief Instantiates an Elevator.
94    ///
95    /// This function instantiates a \ref Elevator.
96    /// \param digraph The digraph, to which we would like to define
97    /// the elevator.
98    /// \param max_level The maximum level of the elevator.
99    static Elevator* createElevator(const Digraph& digraph, int max_level) {
100      return new Elevator(digraph, max_level);
101    }
102
103    /// \brief The tolerance used by the algorithm
104    ///
105    /// The tolerance used by the algorithm to handle inexact computation.
106    typedef lemon::Tolerance<Value> Tolerance;
107
108  };
109
110  ///Push-relabel algorithm for the Network Circulation Problem.
111
112  /**
113     \ingroup max_flow
114     This class implements a push-relabel algorithm
115     or the Network Circulation Problem.
116     The exact formulation of this problem is the following.
117     \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq
118     -delta(v)\quad \forall v\in V \f]
119     \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
120  */
121  template<class _Graph,
122           class _LCapMap=typename _Graph::template ArcMap<int>,
123           class _UCapMap=_LCapMap,
124           class _DeltaMap=typename _Graph::template NodeMap<
125             typename _UCapMap::Value>,
126           class _Traits=CirculationDefaultTraits<_Graph, _LCapMap,
127                                                  _UCapMap, _DeltaMap> >
128  class Circulation {
129
130    typedef _Traits Traits;
131    typedef typename Traits::Digraph Digraph;
132    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
133
134    typedef typename Traits::Value Value;
135
136    typedef typename Traits::LCapMap LCapMap;
137    typedef typename Traits::UCapMap UCapMap;
138    typedef typename Traits::DeltaMap DeltaMap;
139    typedef typename Traits::FlowMap FlowMap;
140    typedef typename Traits::Elevator Elevator;
141    typedef typename Traits::Tolerance Tolerance;
142
143    typedef typename Digraph::template NodeMap<Value> ExcessMap;
144
145    const Digraph &_g;
146    int _node_num;
147
148    const LCapMap *_lo;
149    const UCapMap *_up;
150    const DeltaMap *_delta;
151
152    FlowMap *_flow;
153    bool _local_flow;
154
155    Elevator* _level;
156    bool _local_level;
157
158    ExcessMap* _excess;
159
160    Tolerance _tol;
161    int _el;
162
163  public:
164
165    typedef Circulation Create;
166
167    ///\name Named template parameters
168
169    ///@{
170
171    template <typename _FlowMap>
172    struct DefFlowMapTraits : public Traits {
173      typedef _FlowMap FlowMap;
174      static FlowMap *createFlowMap(const Digraph&) {
175        LEMON_ASSERT(false, "FlowMap is not initialized");
176        return 0; // ignore warnings
177      }
178    };
179
180    /// \brief \ref named-templ-param "Named parameter" for setting
181    /// FlowMap type
182    ///
183    /// \ref named-templ-param "Named parameter" for setting FlowMap
184    /// type
185    template <typename _FlowMap>
186    struct DefFlowMap
187      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
188                           DefFlowMapTraits<_FlowMap> > {
189      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
190                          DefFlowMapTraits<_FlowMap> > Create;
191    };
192
193    template <typename _Elevator>
194    struct DefElevatorTraits : public Traits {
195      typedef _Elevator Elevator;
196      static Elevator *createElevator(const Digraph&, int) {
197        LEMON_ASSERT(false, "Elevator is not initialized");
198        return 0; // ignore warnings
199      }
200    };
201
202    /// \brief \ref named-templ-param "Named parameter" for setting
203    /// Elevator type
204    ///
205    /// \ref named-templ-param "Named parameter" for setting Elevator
206    /// type
207    template <typename _Elevator>
208    struct DefElevator
209      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
210                           DefElevatorTraits<_Elevator> > {
211      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
212                          DefElevatorTraits<_Elevator> > Create;
213    };
214
215    template <typename _Elevator>
216    struct DefStandardElevatorTraits : public Traits {
217      typedef _Elevator Elevator;
218      static Elevator *createElevator(const Digraph& digraph, int max_level) {
219        return new Elevator(digraph, max_level);
220      }
221    };
222
223    /// \brief \ref named-templ-param "Named parameter" for setting
224    /// Elevator type
225    ///
226    /// \ref named-templ-param "Named parameter" for setting Elevator
227    /// type. The Elevator should be standard constructor interface, ie.
228    /// the digraph and the maximum level should be passed to it.
229    template <typename _Elevator>
230    struct DefStandardElevator
231      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
232                       DefStandardElevatorTraits<_Elevator> > {
233      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
234                      DefStandardElevatorTraits<_Elevator> > Create;
235    };
236
237    /// @}
238
239  protected:
240
241    Circulation() {}
242
243  public:
244
245    /// The constructor of the class.
246
247    /// The constructor of the class.
248    /// \param g The digraph the algorithm runs on.
249    /// \param lo The lower bound capacity of the arcs.
250    /// \param up The upper bound capacity of the arcs.
251    /// \param delta The lower bound on node excess.
252    Circulation(const Digraph &g,const LCapMap &lo,
253                const UCapMap &up,const DeltaMap &delta)
254      : _g(g), _node_num(),
255        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
256        _level(0), _local_level(false), _excess(0), _el() {}
257
258    /// Destrcutor.
259    ~Circulation() {
260      destroyStructures();
261    }
262
263  private:
264
265    void createStructures() {
266      _node_num = _el = countNodes(_g);
267
268      if (!_flow) {
269        _flow = Traits::createFlowMap(_g);
270        _local_flow = true;
271      }
272      if (!_level) {
273        _level = Traits::createElevator(_g, _node_num);
274        _local_level = true;
275      }
276      if (!_excess) {
277        _excess = new ExcessMap(_g);
278      }
279    }
280
281    void destroyStructures() {
282      if (_local_flow) {
283        delete _flow;
284      }
285      if (_local_level) {
286        delete _level;
287      }
288      if (_excess) {
289        delete _excess;
290      }
291    }
292
293  public:
294
295    /// Sets the lower bound capacity map.
296
297    /// Sets the lower bound capacity map.
298    /// \return \c (*this)
299    Circulation& lowerCapMap(const LCapMap& map) {
300      _lo = &map;
301      return *this;
302    }
303
304    /// Sets the upper bound capacity map.
305
306    /// Sets the upper bound capacity map.
307    /// \return \c (*this)
308    Circulation& upperCapMap(const LCapMap& map) {
309      _up = &map;
310      return *this;
311    }
312
313    /// Sets the lower bound map on excess.
314
315    /// Sets the lower bound map on excess.
316    /// \return \c (*this)
317    Circulation& deltaMap(const DeltaMap& map) {
318      _delta = &map;
319      return *this;
320    }
321
322    /// Sets the flow map.
323
324    /// Sets the flow map.
325    /// \return \c (*this)
326    Circulation& flowMap(FlowMap& map) {
327      if (_local_flow) {
328        delete _flow;
329        _local_flow = false;
330      }
331      _flow = &map;
332      return *this;
333    }
334
335    /// Returns the flow map.
336
337    /// \return The flow map.
338    ///
339    const FlowMap& flowMap() {
340      return *_flow;
341    }
342
343    /// Sets the elevator.
344
345    /// Sets the elevator.
346    /// \return \c (*this)
347    Circulation& elevator(Elevator& elevator) {
348      if (_local_level) {
349        delete _level;
350        _local_level = false;
351      }
352      _level = &elevator;
353      return *this;
354    }
355
356    /// Returns the elevator.
357
358    /// \return The elevator.
359    ///
360    const Elevator& elevator() {
361      return *_level;
362    }
363
364    /// Sets the tolerance used by algorithm.
365
366    /// Sets the tolerance used by algorithm.
367    ///
368    Circulation& tolerance(const Tolerance& tolerance) const {
369      _tol = tolerance;
370      return *this;
371    }
372
373    /// Returns the tolerance used by algorithm.
374
375    /// Returns the tolerance used by algorithm.
376    ///
377    const Tolerance& tolerance() const {
378      return tolerance;
379    }
380
381    /// \name Execution control
382    /// The simplest way to execute the algorithm is to use one of the
383    /// member functions called \c run().
384    /// \n
385    /// If you need more control on initial solution or execution then
386    /// you have to call one \ref init() function and then the start()
387    /// function.
388
389    ///@{
390
391    /// Initializes the internal data structures.
392
393    /// Initializes the internal data structures. This function sets
394    /// all flow values to the lower bound.
395    /// \return This function returns false if the initialization
396    /// process found a barrier.
397    void init()
398    {
399      createStructures();
400
401      for(NodeIt n(_g);n!=INVALID;++n) {
402        _excess->set(n, (*_delta)[n]);
403      }
404
405      for (ArcIt e(_g);e!=INVALID;++e) {
406        _flow->set(e, (*_lo)[e]);
407        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
408        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
409      }
410
411      // global relabeling tested, but in general case it provides
412      // worse performance for random digraphs
413      _level->initStart();
414      for(NodeIt n(_g);n!=INVALID;++n)
415        _level->initAddItem(n);
416      _level->initFinish();
417      for(NodeIt n(_g);n!=INVALID;++n)
418        if(_tol.positive((*_excess)[n]))
419          _level->activate(n);
420    }
421
422    /// Initializes the internal data structures.
423
424    /// Initializes the internal data structures. This functions uses
425    /// greedy approach to construct the initial solution.
426    void greedyInit()
427    {
428      createStructures();
429
430      for(NodeIt n(_g);n!=INVALID;++n) {
431        _excess->set(n, (*_delta)[n]);
432      }
433
434      for (ArcIt e(_g);e!=INVALID;++e) {
435        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
436          _flow->set(e, (*_up)[e]);
437          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
438          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
439        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
440          _flow->set(e, (*_lo)[e]);
441          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
442          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
443        } else {
444          Value fc = -(*_excess)[_g.target(e)];
445          _flow->set(e, fc);
446          _excess->set(_g.target(e), 0);
447          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
448        }
449      }
450
451      _level->initStart();
452      for(NodeIt n(_g);n!=INVALID;++n)
453        _level->initAddItem(n);
454      _level->initFinish();
455      for(NodeIt n(_g);n!=INVALID;++n)
456        if(_tol.positive((*_excess)[n]))
457          _level->activate(n);
458    }
459
460    ///Starts the algorithm
461
462    ///This function starts the algorithm.
463    ///\return This function returns true if it found a feasible circulation.
464    ///
465    ///\sa barrier()
466    bool start()
467    {
468
469      Node act;
470      Node bact=INVALID;
471      Node last_activated=INVALID;
472      while((act=_level->highestActive())!=INVALID) {
473        int actlevel=(*_level)[act];
474        int mlevel=_node_num;
475        Value exc=(*_excess)[act];
476
477        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
478          Node v = _g.target(e);
479          Value fc=(*_up)[e]-(*_flow)[e];
480          if(!_tol.positive(fc)) continue;
481          if((*_level)[v]<actlevel) {
482            if(!_tol.less(fc, exc)) {
483              _flow->set(e, (*_flow)[e] + exc);
484              _excess->set(v, (*_excess)[v] + exc);
485              if(!_level->active(v) && _tol.positive((*_excess)[v]))
486                _level->activate(v);
487              _excess->set(act,0);
488              _level->deactivate(act);
489              goto next_l;
490            }
491            else {
492              _flow->set(e, (*_up)[e]);
493              _excess->set(v, (*_excess)[v] + fc);
494              if(!_level->active(v) && _tol.positive((*_excess)[v]))
495                _level->activate(v);
496              exc-=fc;
497            }
498          }
499          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
500        }
501        for(InArcIt e(_g,act);e!=INVALID; ++e) {
502          Node v = _g.source(e);
503          Value fc=(*_flow)[e]-(*_lo)[e];
504          if(!_tol.positive(fc)) continue;
505          if((*_level)[v]<actlevel) {
506            if(!_tol.less(fc, exc)) {
507              _flow->set(e, (*_flow)[e] - exc);
508              _excess->set(v, (*_excess)[v] + exc);
509              if(!_level->active(v) && _tol.positive((*_excess)[v]))
510                _level->activate(v);
511              _excess->set(act,0);
512              _level->deactivate(act);
513              goto next_l;
514            }
515            else {
516              _flow->set(e, (*_lo)[e]);
517              _excess->set(v, (*_excess)[v] + fc);
518              if(!_level->active(v) && _tol.positive((*_excess)[v]))
519                _level->activate(v);
520              exc-=fc;
521            }
522          }
523          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
524        }
525
526        _excess->set(act, exc);
527        if(!_tol.positive(exc)) _level->deactivate(act);
528        else if(mlevel==_node_num) {
529          _level->liftHighestActiveToTop();
530          _el = _node_num;
531          return false;
532        }
533        else {
534          _level->liftHighestActive(mlevel+1);
535          if(_level->onLevel(actlevel)==0) {
536            _el = actlevel;
537            return false;
538          }
539        }
540      next_l:
541        ;
542      }
543      return true;
544    }
545
546    /// Runs the circulation algorithm.
547
548    /// Runs the circulation algorithm.
549    /// \note fc.run() is just a shortcut of the following code.
550    /// \code
551    ///   fc.greedyInit();
552    ///   return fc.start();
553    /// \endcode
554    bool run() {
555      greedyInit();
556      return start();
557    }
558
559    /// @}
560
561    /// \name Query Functions
562    /// The result of the %Circulation algorithm can be obtained using
563    /// these functions.
564    /// \n
565    /// Before the use of these functions,
566    /// either run() or start() must be called.
567
568    ///@{
569
570    /**
571       \brief Returns a barrier
572       
573       Barrier is a set \e B of nodes for which
574       \f[ \sum_{v\in B}-delta(v)<
575       \sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
576       holds. The existence of a set with this property prooves that a feasible
577       flow cannot exists.
578       \sa checkBarrier()
579       \sa run()
580    */
581    template<class GT>
582    void barrierMap(GT &bar)
583    {
584      for(NodeIt n(_g);n!=INVALID;++n)
585        bar.set(n, (*_level)[n] >= _el);
586    }
587
588    ///Returns true if the node is in the barrier
589
590    ///Returns true if the node is in the barrier
591    ///\sa barrierMap()
592    bool barrier(const Node& node)
593    {
594      return (*_level)[node] >= _el;
595    }
596
597    /// \brief Returns the flow on the arc.
598    ///
599    /// Sets the \c flowMap to the flow on the arcs. This method can
600    /// be called after the second phase of algorithm.
601    Value flow(const Arc& arc) const {
602      return (*_flow)[arc];
603    }
604
605    /// @}
606
607    /// \name Checker Functions
608    /// The feasibility  of the results can be checked using
609    /// these functions.
610    /// \n
611    /// Before the use of these functions,
612    /// either run() or start() must be called.
613
614    ///@{
615
616    ///Check if the  \c flow is a feasible circulation
617    bool checkFlow() {
618      for(ArcIt e(_g);e!=INVALID;++e)
619        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
620      for(NodeIt n(_g);n!=INVALID;++n)
621        {
622          Value dif=-(*_delta)[n];
623          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
624          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
625          if(_tol.negative(dif)) return false;
626        }
627      return true;
628    }
629
630    ///Check whether or not the last execution provides a barrier
631
632    ///Check whether or not the last execution provides a barrier
633    ///\sa barrier()
634    bool checkBarrier()
635    {
636      Value delta=0;
637      for(NodeIt n(_g);n!=INVALID;++n)
638        if(barrier(n))
639          delta-=(*_delta)[n];
640      for(ArcIt e(_g);e!=INVALID;++e)
641        {
642          Node s=_g.source(e);
643          Node t=_g.target(e);
644          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
645          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
646        }
647      return _tol.negative(delta);
648    }
649
650    /// @}
651
652  };
653
654}
655
656#endif
Note: See TracBrowser for help on using the repository browser.