Last change on this file since 2375:e30a0fdad0d7 was 2375:e30a0fdad0d7, checked in by Alpar Juttner, 17 years ago

A preflow based general network circulation algorithm and a simple demo

File size: 7.8 KB
Line
1/* -*- C++ -*-
2 * lemon/preflow_matching.h - Part of LEMON, a generic C++ optimization library
3 *
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5 * (Egervary Research Group on Combinatorial Optimization, EGRES).
6 *
7 * Permission to use, modify and distribute this software is granted
8 * provided that this copyright notice appears in all copies. For
9 * precise terms see the accompanying LICENSE file.
10 *
11 * This software is provided "AS IS" with no warranty of any kind,
12 * express or implied, and with no claim as to its suitability for any
13 * purpose.
14 *
15 */
16
17#ifndef LEMON_CIRCULATION_H
18#define LEMON_CIRCULATION_H
19
20#include <lemon/graph_utils.h>
21#include <iostream>
22#include <queue>
23#include <lemon/tolerance.h>
24#include <lemon/elevator.h>
25
26///\ingroup flowalgs
27///\file
28///\brief Push-prelabel algorithm for finding a feasible circulation.
29///
30namespace lemon {
31
32  ///Preflow algorithm for the Network Circulation Problem.
33
34  ///\ingroup flowalgs
35  ///This class implements a preflow algorithm
36  ///for the Network Circulation Problem.
37  ///The exact formulation of this problem is the following.
38  /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq delta(v)\quad \forall v\in V \f]
39  /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
40  ///
41  template<class Graph,
42           class Value,
43           class FlowMap=typename Graph::template EdgeMap<Value>,
44           class LCapMap=typename Graph::template EdgeMap<Value>,
45           class UCapMap=LCapMap,
46           class DeltaMap=typename Graph::template NodeMap<Value>
47            >
48  class Circulation {
49    typedef typename Graph::Node Node;
50    typedef typename Graph::NodeIt NodeIt;
51    typedef typename Graph::Edge Edge;
52    typedef typename Graph::EdgeIt EdgeIt;
53    typedef typename Graph::InEdgeIt InEdgeIt;
54    typedef typename Graph::OutEdgeIt OutEdgeIt;
55
56
57    const Graph &_g;
58    int _node_num;
59
60    const LCapMap &_lo;
61    const UCapMap &_up;
62    const DeltaMap &_delta;
63    FlowMap &_x;
64    Tolerance<Value> _tol;
65    Elevator<Graph,typename Graph::Node> _levels;
66    typename Graph::template NodeMap<Value> _excess;
67
68  public:
69    ///\e
70    Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up,
71                const DeltaMap &delta,FlowMap &x) :
72      _g(g),
73      _node_num(countNodes(g)),
74      _lo(lo),_up(up),_delta(delta),_x(x),
75      _levels(g,_node_num),
76      _excess(g)
77    {
78    }
79
80  private:
81
83    {
84      if(_tol.positive(_excess[n]+=v))
85        {
86          if(!_levels.active(n)) _levels.activate(n);
87        }
88      else if(_levels.active(n)) _levels.deactivate(n);
89    }
90
91    void init()
92    {
93
94      _x=_lo;
95
96      for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=-_delta[n];
97
98      for(EdgeIt e(_g);e!=INVALID;++e)
99        {
100          _excess[_g.target(e)]+=_x[e];
101          _excess[_g.source(e)]-=_x[e];
102        }
103
104      _levels.initStart();
105      for(NodeIt n(_g);n!=INVALID;++n)
107      _levels.initFinish();
108      for(NodeIt n(_g);n!=INVALID;++n)
109        if(_tol.positive(_excess[n]))
110          _levels.activate(n);
111    }
112
113  public:
114    ///Check if \c x is a feasible circulation
115    template<class FT>
116    bool checkX(FT &x) {
117      for(EdgeIt e(_g);e!=INVALID;++e)
118        if(x[e]<_lo[e]||x[e]>_up[e]) return false;
119      for(NodeIt n(_g);n!=INVALID;++n)
120        {
121          Value dif=_delta[n];
122          for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=x[e];
123          for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=x[e];
124          if(_tol.negative(dif)) return false;
125        }
126      return true;
127    };
128    ///Check if the default \c x is a feasible circulation
129    bool checkX() { return checkX(_x); }
130
131    ///Chech if \c bar is a real barrier
132
133    ///Chech if \c bar is a real barrier
134    ///\sa barrier()
135    template<class GT>
136    bool checkBarrier(GT &bar)
137    {
138      Value delta=0;
139      for(NodeIt n(_g);n!=INVALID;++n)
140        if(bar[n])
141          delta+=_delta[n];
142      for(EdgeIt e(_g);e!=INVALID;++e)
143        {
144          Node s=_g.source(e);
145          Node t=_g.target(e);
146          if(bar[s]&&!bar[t]) delta+=_up[e];
147          else if(bar[t]&&!bar[s]) delta-=_lo[e];
148        }
149      return _tol.negative(delta);
150    }
151    ///Chech whether or not the last execution provides a barrier
152
153    ///Chech whether or not the last execution provides a barrier
154    ///\sa barrier()
155    bool checkBarrier()
156    {
157      typename Graph:: template NodeMap<bool> bar(_g);
158      barrier(bar);
159      return checkBarrier(bar);
160    }
161    ///Run the algorithm
162
163    ///This function runs the algorithm.
164    ///\return This function returns -1 if it found a feasible circulation.
165    /// nonnegative values (including 0) mean that no feasible solution is
166    /// found. In this case the return value means an "empty level".
167    ///
168    ///\sa barrier()
169    int run()
170    {
171      init();
172
173#ifdef LEMON_CIRCULATION_DEBUG
174      for(NodeIt n(_g);n!=INVALID;++n)
175        std::cerr<< _levels[n] << ' ';
176      std::cerr << std::endl;
177#endif
178      Node act;
179      Node bact=INVALID;
180      Node last_activated=INVALID;
181      while((act=_levels.highestActive())!=INVALID) {
182        int actlevel=_levels[act];
183        int tlevel;
184        int mlevel=_node_num;
185        Value exc=_excess[act];
186        Value fc;
187
188#ifdef LEMON_CIRCULATION_DEBUG
189        for(NodeIt n(_g);n!=INVALID;++n)
190          std::cerr<< _levels[n] << ' ';
191        std::cerr << std::endl;
192        std::cerr << "Process node " << _g.id(act)
193                  << " on level " << actlevel
194                  << " with excess " << exc
195                  << std::endl;
196#endif
197        for(OutEdgeIt e(_g,act);e!=INVALID; ++e)
198          if((fc=_up[e]-_x[e])>0)
199            if((tlevel=_levels[_g.target(e)])<actlevel)
200              if(fc<=exc) {
201                _x[e]=_up[e];
203                exc-=fc;
204#ifdef LEMON_CIRCULATION_DEBUG
205                std::cerr << "  Push " << fc
206                          << " toward " << _g.id(_g.target(e)) << std::endl;
207#endif
208              }
209              else {
210                _x[e]+=exc;
212                //exc=0;
213                _excess[act]=0;
214                _levels.deactivate(act);
215#ifdef LEMON_CIRCULATION_DEBUG
216                std::cerr << "  Push " << exc
217                          << " toward " << _g.id(_g.target(e)) << std::endl;
218                std::cerr << "  Deactivate\n";
219#endif
220                goto next_l;
221              }
222            else if(tlevel<mlevel) mlevel=tlevel;
223
224        for(InEdgeIt e(_g,act);e!=INVALID; ++e)
225          if((fc=_x[e]-_lo[e])>0)
226            if((tlevel=_levels[_g.source(e)])<actlevel)
227              if(fc<=exc) {
228                _x[e]=_lo[e];
230                exc-=fc;
231#ifdef LEMON_CIRCULATION_DEBUG
232                std::cerr << "  Push " << fc
233                          << " toward " << _g.id(_g.source(e)) << std::endl;
234#endif
235              }
236              else {
237                _x[e]-=exc;
239                //exc=0;
240                _excess[act]=0;
241                _levels.deactivate(act);
242#ifdef LEMON_CIRCULATION_DEBUG
243                std::cerr << "  Push " << exc
244                          << " toward " << _g.id(_g.source(e)) << std::endl;
245                std::cerr << "  Deactivate\n";
246#endif
247                goto next_l;
248              }
249            else if(tlevel<mlevel) mlevel=tlevel;
250
251        _excess[act]=exc;
252        if(!_tol.positive(exc)) _levels.deactivate(act);
253        else if(mlevel==_node_num) {
254          _levels.liftHighestActiveTo(_node_num);
255#ifdef LEMON_CIRCULATION_DEBUG
256          std::cerr << "  Lift to level " << _node_num << std::endl;
257#endif
258          return _levels.onLevel(_node_num-1)==0?_node_num-1:actlevel;
259        }
260        else {
261          _levels.liftHighestActiveTo(mlevel+1);
262#ifdef LEMON_CIRCULATION_DEBUG
263          std::cerr << "  Lift to level " << mlevel+1 << std::endl;
264#endif
265          if(_levels.onLevel(actlevel)==0)
266            return actlevel;
267        }
268      next_l:
269        ;
270      }
271#ifdef LEMON_CIRCULATION_DEBUG
272      std::cerr << "Feasible flow found.\n";
273#endif
274      return -1;
275    }
276
277    ///Return a barrier
278
279    ///Barrier is a set \e B of nodes for which
280    /// \f[ \sum_{v\in B}delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
281    ///holds. The existence of a set with this property prooves that a feasible
282    ///flow cannot exists.
283    ///\pre The run() must have been executed, and its return value was -1.
284    ///\sa checkBarrier()
285    ///\sa run()
286    template<class GT>
287    void barrier(GT &bar,int empty_level=-1)
288    {
289      if(empty_level==-1)
290        for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ;
291      for(NodeIt n(_g);n!=INVALID;++n)
292        bar[n] = _levels[n]>empty_level;
293    }
294
295  };
296
297}
298
299#endif
Note: See TracBrowser for help on using the repository browser.