lemon/circulation.h
author Alpar Juttner <alpar@cs.elte.hu>
Mon, 01 Dec 2008 14:11:31 +0000
changeset 400 fa341dd6ab23
child 401 26fd85a3087e
permissions -rw-r--r--
Transform circulation demo to test
     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 ///
    31 namespace 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