- Gamma distributon random variable.
- Test file for random.h
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2007
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
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.
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
19 #ifndef LEMON_CIRCULATION_H
20 #define LEMON_CIRCULATION_H
22 #include <lemon/graph_utils.h>
25 #include <lemon/tolerance.h>
26 #include <lemon/elevator.h>
30 ///\brief Push-prelabel algorithm for finding a feasible circulation.
34 ///Preflow algorithm for the Network Circulation Problem.
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]
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>
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;
64 const DeltaMap &_delta;
66 Tolerance<Value> _tol;
67 Elevator<Graph,typename Graph::Node> _levels;
68 typename Graph::template NodeMap<Value> _excess;
72 Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up,
73 const DeltaMap &delta,FlowMap &x) :
75 _node_num(countNodes(g)),
76 _lo(lo),_up(up),_delta(delta),_x(x),
84 void addExcess(Node n,Value v)
86 if(_tol.positive(_excess[n]+=v))
88 if(!_levels.active(n)) _levels.activate(n);
90 else if(_levels.active(n)) _levels.deactivate(n);
98 for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=_delta[n];
100 for(EdgeIt e(_g);e!=INVALID;++e)
102 _excess[_g.target(e)]+=_x[e];
103 _excess[_g.source(e)]-=_x[e];
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]))
116 ///Check if \c x is a feasible circulation
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)
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;
130 ///Check if the default \c x is a feasible circulation
131 bool checkX() { return checkX(_x); }
133 ///Check if \c bar is a real barrier
135 ///Check if \c bar is a real barrier
138 bool checkBarrier(GT &bar)
141 for(NodeIt n(_g);n!=INVALID;++n)
144 for(EdgeIt e(_g);e!=INVALID;++e)
148 if(bar[s]&&!bar[t]) delta+=_up[e];
149 else if(bar[t]&&!bar[s]) delta-=_lo[e];
151 return _tol.negative(delta);
153 ///Check whether or not the last execution provides a barrier
155 ///Check whether or not the last execution provides a barrier
159 typename Graph:: template NodeMap<bool> bar(_g);
161 return checkBarrier(bar);
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".
175 #ifdef LEMON_CIRCULATION_DEBUG
176 for(NodeIt n(_g);n!=INVALID;++n)
177 std::cerr<< _levels[n] << ' ';
178 std::cerr << std::endl;
182 Node last_activated=INVALID;
183 while((act=_levels.highestActive())!=INVALID) {
184 int actlevel=_levels[act];
186 int mlevel=_node_num;
187 Value exc=_excess[act];
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
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)
204 addExcess(_g.target(e),fc);
206 #ifdef LEMON_CIRCULATION_DEBUG
207 std::cerr << " Push " << fc
208 << " toward " << _g.id(_g.target(e)) << std::endl;
213 addExcess(_g.target(e),exc);
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";
224 else if(tlevel<mlevel) mlevel=tlevel;
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)
231 addExcess(_g.source(e),fc);
233 #ifdef LEMON_CIRCULATION_DEBUG
234 std::cerr << " Push " << fc
235 << " toward " << _g.id(_g.source(e)) << std::endl;
240 addExcess(_g.source(e),exc);
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";
251 else if(tlevel<mlevel) mlevel=tlevel;
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;
260 return _levels.onLevel(_node_num-1)==0?_node_num-1:actlevel;
263 _levels.liftHighestActiveTo(mlevel+1);
264 #ifdef LEMON_CIRCULATION_DEBUG
265 std::cerr << " Lift to level " << mlevel+1 << std::endl;
267 if(_levels.onLevel(actlevel)==0)
273 #ifdef LEMON_CIRCULATION_DEBUG
274 std::cerr << "Feasible flow found.\n";
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()
289 void barrier(GT &bar,int 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;