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