deba@2440
|
1 |
/* -*- C++ -*-
|
deba@2440
|
2 |
*
|
deba@2440
|
3 |
* This file is a part of LEMON, a generic C++ optimization library
|
deba@2440
|
4 |
*
|
alpar@2553
|
5 |
* Copyright (C) 2003-2008
|
deba@2440
|
6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
|
deba@2440
|
7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES).
|
deba@2440
|
8 |
*
|
deba@2440
|
9 |
* Permission to use, modify and distribute this software is granted
|
deba@2440
|
10 |
* provided that this copyright notice appears in all copies. For
|
deba@2440
|
11 |
* precise terms see the accompanying LICENSE file.
|
deba@2440
|
12 |
*
|
deba@2440
|
13 |
* This software is provided "AS IS" with no warranty of any kind,
|
deba@2440
|
14 |
* express or implied, and with no claim as to its suitability for any
|
deba@2440
|
15 |
* purpose.
|
deba@2440
|
16 |
*
|
deba@2440
|
17 |
*/
|
deba@2440
|
18 |
|
deba@2440
|
19 |
#ifndef LEMON_CAPACITY_SCALING_H
|
deba@2440
|
20 |
#define LEMON_CAPACITY_SCALING_H
|
deba@2440
|
21 |
|
deba@2440
|
22 |
/// \ingroup min_cost_flow
|
deba@2440
|
23 |
///
|
deba@2440
|
24 |
/// \file
|
kpeter@2574
|
25 |
/// \brief Capacity scaling algorithm for finding a minimum cost flow.
|
kpeter@2574
|
26 |
|
kpeter@2574
|
27 |
#include <vector>
|
deba@2440
|
28 |
|
kpeter@2535
|
29 |
#include <lemon/graph_adaptor.h>
|
kpeter@2535
|
30 |
#include <lemon/bin_heap.h>
|
deba@2457
|
31 |
|
deba@2440
|
32 |
namespace lemon {
|
deba@2440
|
33 |
|
deba@2440
|
34 |
/// \addtogroup min_cost_flow
|
deba@2440
|
35 |
/// @{
|
deba@2440
|
36 |
|
kpeter@2574
|
37 |
/// \brief Implementation of the capacity scaling algorithm for
|
kpeter@2574
|
38 |
/// finding a minimum cost flow.
|
deba@2440
|
39 |
///
|
kpeter@2535
|
40 |
/// \ref CapacityScaling implements the capacity scaling version
|
kpeter@2535
|
41 |
/// of the successive shortest path algorithm for finding a minimum
|
kpeter@2535
|
42 |
/// cost flow.
|
deba@2440
|
43 |
///
|
kpeter@2574
|
44 |
/// \tparam Graph The directed graph type the algorithm runs on.
|
kpeter@2574
|
45 |
/// \tparam LowerMap The type of the lower bound map.
|
kpeter@2574
|
46 |
/// \tparam CapacityMap The type of the capacity (upper bound) map.
|
kpeter@2574
|
47 |
/// \tparam CostMap The type of the cost (length) map.
|
kpeter@2574
|
48 |
/// \tparam SupplyMap The type of the supply map.
|
deba@2440
|
49 |
///
|
deba@2440
|
50 |
/// \warning
|
kpeter@2574
|
51 |
/// - Edge capacities and costs should be \e non-negative \e integers.
|
kpeter@2574
|
52 |
/// - Supply values should be \e signed \e integers.
|
kpeter@2574
|
53 |
/// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
|
kpeter@2574
|
54 |
/// - \c CapacityMap::Value and \c SupplyMap::Value must be
|
kpeter@2574
|
55 |
/// convertible to each other.
|
kpeter@2574
|
56 |
/// - All value types must be convertible to \c CostMap::Value, which
|
kpeter@2574
|
57 |
/// must be signed type.
|
deba@2440
|
58 |
///
|
deba@2440
|
59 |
/// \author Peter Kovacs
|
deba@2440
|
60 |
|
kpeter@2533
|
61 |
template < typename Graph,
|
kpeter@2535
|
62 |
typename LowerMap = typename Graph::template EdgeMap<int>,
|
kpeter@2574
|
63 |
typename CapacityMap = typename Graph::template EdgeMap<int>,
|
kpeter@2535
|
64 |
typename CostMap = typename Graph::template EdgeMap<int>,
|
kpeter@2574
|
65 |
typename SupplyMap = typename Graph::template NodeMap<int> >
|
deba@2440
|
66 |
class CapacityScaling
|
deba@2440
|
67 |
{
|
kpeter@2556
|
68 |
GRAPH_TYPEDEFS(typename Graph);
|
deba@2440
|
69 |
|
deba@2440
|
70 |
typedef typename CapacityMap::Value Capacity;
|
deba@2440
|
71 |
typedef typename CostMap::Value Cost;
|
deba@2440
|
72 |
typedef typename SupplyMap::Value Supply;
|
kpeter@2556
|
73 |
typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
|
kpeter@2556
|
74 |
typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
|
kpeter@2535
|
75 |
typedef typename Graph::template NodeMap<Edge> PredMap;
|
deba@2440
|
76 |
|
deba@2440
|
77 |
public:
|
deba@2440
|
78 |
|
kpeter@2556
|
79 |
/// The type of the flow map.
|
kpeter@2556
|
80 |
typedef typename Graph::template EdgeMap<Capacity> FlowMap;
|
kpeter@2556
|
81 |
/// The type of the potential map.
|
deba@2440
|
82 |
typedef typename Graph::template NodeMap<Cost> PotentialMap;
|
deba@2440
|
83 |
|
kpeter@2574
|
84 |
private:
|
deba@2440
|
85 |
|
kpeter@2535
|
86 |
/// \brief Special implementation of the \ref Dijkstra algorithm
|
kpeter@2574
|
87 |
/// for finding shortest paths in the residual network.
|
kpeter@2574
|
88 |
///
|
kpeter@2574
|
89 |
/// \ref ResidualDijkstra is a special implementation of the
|
kpeter@2574
|
90 |
/// \ref Dijkstra algorithm for finding shortest paths in the
|
kpeter@2574
|
91 |
/// residual network of the graph with respect to the reduced edge
|
kpeter@2574
|
92 |
/// costs and modifying the node potentials according to the
|
kpeter@2574
|
93 |
/// distance of the nodes.
|
kpeter@2535
|
94 |
class ResidualDijkstra
|
deba@2440
|
95 |
{
|
kpeter@2535
|
96 |
typedef typename Graph::template NodeMap<Cost> CostNodeMap;
|
kpeter@2535
|
97 |
typedef typename Graph::template NodeMap<Edge> PredMap;
|
deba@2440
|
98 |
|
kpeter@2535
|
99 |
typedef typename Graph::template NodeMap<int> HeapCrossRef;
|
kpeter@2535
|
100 |
typedef BinHeap<Cost, HeapCrossRef> Heap;
|
kpeter@2535
|
101 |
|
kpeter@2574
|
102 |
private:
|
kpeter@2535
|
103 |
|
kpeter@2574
|
104 |
// The directed graph the algorithm runs on
|
kpeter@2574
|
105 |
const Graph &_graph;
|
kpeter@2535
|
106 |
|
kpeter@2574
|
107 |
// The main maps
|
kpeter@2574
|
108 |
const FlowMap &_flow;
|
kpeter@2574
|
109 |
const CapacityEdgeMap &_res_cap;
|
kpeter@2574
|
110 |
const CostMap &_cost;
|
kpeter@2574
|
111 |
const SupplyNodeMap &_excess;
|
kpeter@2574
|
112 |
PotentialMap &_potential;
|
kpeter@2535
|
113 |
|
kpeter@2574
|
114 |
// The distance map
|
kpeter@2574
|
115 |
CostNodeMap _dist;
|
kpeter@2574
|
116 |
// The pred edge map
|
kpeter@2574
|
117 |
PredMap &_pred;
|
kpeter@2574
|
118 |
// The processed (i.e. permanently labeled) nodes
|
kpeter@2574
|
119 |
std::vector<Node> _proc_nodes;
|
deba@2440
|
120 |
|
deba@2440
|
121 |
public:
|
deba@2440
|
122 |
|
kpeter@2556
|
123 |
/// The constructor of the class.
|
kpeter@2574
|
124 |
ResidualDijkstra( const Graph &graph,
|
kpeter@2574
|
125 |
const FlowMap &flow,
|
kpeter@2574
|
126 |
const CapacityEdgeMap &res_cap,
|
kpeter@2574
|
127 |
const CostMap &cost,
|
kpeter@2574
|
128 |
const SupplyMap &excess,
|
kpeter@2574
|
129 |
PotentialMap &potential,
|
kpeter@2574
|
130 |
PredMap &pred ) :
|
kpeter@2574
|
131 |
_graph(graph), _flow(flow), _res_cap(res_cap), _cost(cost),
|
kpeter@2574
|
132 |
_excess(excess), _potential(potential), _dist(graph),
|
kpeter@2574
|
133 |
_pred(pred)
|
kpeter@2535
|
134 |
{}
|
deba@2440
|
135 |
|
kpeter@2556
|
136 |
/// Runs the algorithm from the given source node.
|
kpeter@2535
|
137 |
Node run(Node s, Capacity delta) {
|
kpeter@2574
|
138 |
HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
|
kpeter@2535
|
139 |
Heap heap(heap_cross_ref);
|
kpeter@2535
|
140 |
heap.push(s, 0);
|
kpeter@2574
|
141 |
_pred[s] = INVALID;
|
kpeter@2574
|
142 |
_proc_nodes.clear();
|
kpeter@2535
|
143 |
|
kpeter@2535
|
144 |
// Processing nodes
|
kpeter@2574
|
145 |
while (!heap.empty() && _excess[heap.top()] > -delta) {
|
kpeter@2535
|
146 |
Node u = heap.top(), v;
|
kpeter@2574
|
147 |
Cost d = heap.prio() + _potential[u], nd;
|
kpeter@2574
|
148 |
_dist[u] = heap.prio();
|
kpeter@2535
|
149 |
heap.pop();
|
kpeter@2574
|
150 |
_proc_nodes.push_back(u);
|
kpeter@2535
|
151 |
|
kpeter@2535
|
152 |
// Traversing outgoing edges
|
kpeter@2574
|
153 |
for (OutEdgeIt e(_graph, u); e != INVALID; ++e) {
|
kpeter@2574
|
154 |
if (_res_cap[e] >= delta) {
|
kpeter@2574
|
155 |
v = _graph.target(e);
|
kpeter@2535
|
156 |
switch(heap.state(v)) {
|
kpeter@2535
|
157 |
case Heap::PRE_HEAP:
|
kpeter@2574
|
158 |
heap.push(v, d + _cost[e] - _potential[v]);
|
kpeter@2574
|
159 |
_pred[v] = e;
|
kpeter@2535
|
160 |
break;
|
kpeter@2535
|
161 |
case Heap::IN_HEAP:
|
kpeter@2574
|
162 |
nd = d + _cost[e] - _potential[v];
|
kpeter@2535
|
163 |
if (nd < heap[v]) {
|
kpeter@2535
|
164 |
heap.decrease(v, nd);
|
kpeter@2574
|
165 |
_pred[v] = e;
|
kpeter@2535
|
166 |
}
|
kpeter@2535
|
167 |
break;
|
kpeter@2535
|
168 |
case Heap::POST_HEAP:
|
kpeter@2535
|
169 |
break;
|
kpeter@2535
|
170 |
}
|
kpeter@2535
|
171 |
}
|
kpeter@2535
|
172 |
}
|
kpeter@2535
|
173 |
|
kpeter@2535
|
174 |
// Traversing incoming edges
|
kpeter@2574
|
175 |
for (InEdgeIt e(_graph, u); e != INVALID; ++e) {
|
kpeter@2574
|
176 |
if (_flow[e] >= delta) {
|
kpeter@2574
|
177 |
v = _graph.source(e);
|
kpeter@2535
|
178 |
switch(heap.state(v)) {
|
kpeter@2535
|
179 |
case Heap::PRE_HEAP:
|
kpeter@2574
|
180 |
heap.push(v, d - _cost[e] - _potential[v]);
|
kpeter@2574
|
181 |
_pred[v] = e;
|
kpeter@2535
|
182 |
break;
|
kpeter@2535
|
183 |
case Heap::IN_HEAP:
|
kpeter@2574
|
184 |
nd = d - _cost[e] - _potential[v];
|
kpeter@2535
|
185 |
if (nd < heap[v]) {
|
kpeter@2535
|
186 |
heap.decrease(v, nd);
|
kpeter@2574
|
187 |
_pred[v] = e;
|
kpeter@2535
|
188 |
}
|
kpeter@2535
|
189 |
break;
|
kpeter@2535
|
190 |
case Heap::POST_HEAP:
|
kpeter@2535
|
191 |
break;
|
kpeter@2535
|
192 |
}
|
kpeter@2535
|
193 |
}
|
kpeter@2535
|
194 |
}
|
kpeter@2535
|
195 |
}
|
kpeter@2535
|
196 |
if (heap.empty()) return INVALID;
|
kpeter@2535
|
197 |
|
kpeter@2535
|
198 |
// Updating potentials of processed nodes
|
kpeter@2535
|
199 |
Node t = heap.top();
|
kpeter@2574
|
200 |
Cost t_dist = heap.prio();
|
kpeter@2574
|
201 |
for (int i = 0; i < int(_proc_nodes.size()); ++i)
|
kpeter@2574
|
202 |
_potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
|
kpeter@2535
|
203 |
|
kpeter@2535
|
204 |
return t;
|
deba@2440
|
205 |
}
|
deba@2440
|
206 |
|
kpeter@2535
|
207 |
}; //class ResidualDijkstra
|
deba@2440
|
208 |
|
kpeter@2574
|
209 |
private:
|
deba@2440
|
210 |
|
kpeter@2574
|
211 |
// The directed graph the algorithm runs on
|
kpeter@2574
|
212 |
const Graph &_graph;
|
kpeter@2574
|
213 |
// The original lower bound map
|
kpeter@2574
|
214 |
const LowerMap *_lower;
|
kpeter@2574
|
215 |
// The modified capacity map
|
kpeter@2574
|
216 |
CapacityEdgeMap _capacity;
|
kpeter@2574
|
217 |
// The original cost map
|
kpeter@2574
|
218 |
const CostMap &_cost;
|
kpeter@2574
|
219 |
// The modified supply map
|
kpeter@2574
|
220 |
SupplyNodeMap _supply;
|
kpeter@2574
|
221 |
bool _valid_supply;
|
deba@2440
|
222 |
|
kpeter@2574
|
223 |
// Edge map of the current flow
|
kpeter@2574
|
224 |
FlowMap _flow;
|
kpeter@2574
|
225 |
// Node map of the current potentials
|
kpeter@2574
|
226 |
PotentialMap _potential;
|
deba@2440
|
227 |
|
kpeter@2574
|
228 |
// The residual capacity map
|
kpeter@2574
|
229 |
CapacityEdgeMap _res_cap;
|
kpeter@2574
|
230 |
// The excess map
|
kpeter@2574
|
231 |
SupplyNodeMap _excess;
|
kpeter@2574
|
232 |
// The excess nodes (i.e. nodes with positive excess)
|
kpeter@2574
|
233 |
std::vector<Node> _excess_nodes;
|
kpeter@2574
|
234 |
// The deficit nodes (i.e. nodes with negative excess)
|
kpeter@2574
|
235 |
std::vector<Node> _deficit_nodes;
|
deba@2440
|
236 |
|
kpeter@2574
|
237 |
// The delta parameter used for capacity scaling
|
kpeter@2574
|
238 |
Capacity _delta;
|
kpeter@2574
|
239 |
// The maximum number of phases
|
kpeter@2574
|
240 |
int _phase_num;
|
deba@2440
|
241 |
|
kpeter@2574
|
242 |
// The pred edge map
|
kpeter@2574
|
243 |
PredMap _pred;
|
kpeter@2574
|
244 |
// Implementation of the Dijkstra algorithm for finding augmenting
|
kpeter@2574
|
245 |
// shortest paths in the residual network
|
kpeter@2574
|
246 |
ResidualDijkstra _dijkstra;
|
deba@2440
|
247 |
|
deba@2440
|
248 |
public :
|
deba@2440
|
249 |
|
deba@2440
|
250 |
/// \brief General constructor of the class (with lower bounds).
|
deba@2440
|
251 |
///
|
deba@2440
|
252 |
/// General constructor of the class (with lower bounds).
|
deba@2440
|
253 |
///
|
kpeter@2574
|
254 |
/// \param graph The directed graph the algorithm runs on.
|
kpeter@2574
|
255 |
/// \param lower The lower bounds of the edges.
|
kpeter@2574
|
256 |
/// \param capacity The capacities (upper bounds) of the edges.
|
kpeter@2574
|
257 |
/// \param cost The cost (length) values of the edges.
|
kpeter@2574
|
258 |
/// \param supply The supply values of the nodes (signed).
|
kpeter@2574
|
259 |
CapacityScaling( const Graph &graph,
|
kpeter@2574
|
260 |
const LowerMap &lower,
|
kpeter@2574
|
261 |
const CapacityMap &capacity,
|
kpeter@2574
|
262 |
const CostMap &cost,
|
kpeter@2574
|
263 |
const SupplyMap &supply ) :
|
kpeter@2574
|
264 |
_graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
|
kpeter@2574
|
265 |
_supply(graph), _flow(graph, 0), _potential(graph, 0),
|
kpeter@2574
|
266 |
_res_cap(graph), _excess(graph), _pred(graph),
|
kpeter@2574
|
267 |
_dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
|
deba@2440
|
268 |
{
|
kpeter@2556
|
269 |
// Removing non-zero lower bounds
|
kpeter@2574
|
270 |
_capacity = subMap(capacity, lower);
|
kpeter@2574
|
271 |
_res_cap = _capacity;
|
deba@2440
|
272 |
Supply sum = 0;
|
kpeter@2574
|
273 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@2574
|
274 |
Supply s = supply[n];
|
kpeter@2574
|
275 |
for (InEdgeIt e(_graph, n); e != INVALID; ++e)
|
kpeter@2574
|
276 |
s += lower[e];
|
kpeter@2574
|
277 |
for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
|
kpeter@2574
|
278 |
s -= lower[e];
|
kpeter@2574
|
279 |
_supply[n] = s;
|
kpeter@2535
|
280 |
sum += s;
|
deba@2440
|
281 |
}
|
kpeter@2574
|
282 |
_valid_supply = sum == 0;
|
deba@2440
|
283 |
}
|
deba@2440
|
284 |
|
deba@2440
|
285 |
/// \brief General constructor of the class (without lower bounds).
|
deba@2440
|
286 |
///
|
deba@2440
|
287 |
/// General constructor of the class (without lower bounds).
|
deba@2440
|
288 |
///
|
kpeter@2574
|
289 |
/// \param graph The directed graph the algorithm runs on.
|
kpeter@2574
|
290 |
/// \param capacity The capacities (upper bounds) of the edges.
|
kpeter@2574
|
291 |
/// \param cost The cost (length) values of the edges.
|
kpeter@2574
|
292 |
/// \param supply The supply values of the nodes (signed).
|
kpeter@2574
|
293 |
CapacityScaling( const Graph &graph,
|
kpeter@2574
|
294 |
const CapacityMap &capacity,
|
kpeter@2574
|
295 |
const CostMap &cost,
|
kpeter@2574
|
296 |
const SupplyMap &supply ) :
|
kpeter@2574
|
297 |
_graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
|
kpeter@2574
|
298 |
_supply(supply), _flow(graph, 0), _potential(graph, 0),
|
kpeter@2574
|
299 |
_res_cap(capacity), _excess(graph), _pred(graph),
|
kpeter@2574
|
300 |
_dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
|
deba@2440
|
301 |
{
|
deba@2440
|
302 |
// Checking the sum of supply values
|
deba@2440
|
303 |
Supply sum = 0;
|
kpeter@2574
|
304 |
for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
|
kpeter@2574
|
305 |
_valid_supply = sum == 0;
|
deba@2440
|
306 |
}
|
deba@2440
|
307 |
|
deba@2440
|
308 |
/// \brief Simple constructor of the class (with lower bounds).
|
deba@2440
|
309 |
///
|
deba@2440
|
310 |
/// Simple constructor of the class (with lower bounds).
|
deba@2440
|
311 |
///
|
kpeter@2574
|
312 |
/// \param graph The directed graph the algorithm runs on.
|
kpeter@2574
|
313 |
/// \param lower The lower bounds of the edges.
|
kpeter@2574
|
314 |
/// \param capacity The capacities (upper bounds) of the edges.
|
kpeter@2574
|
315 |
/// \param cost The cost (length) values of the edges.
|
kpeter@2574
|
316 |
/// \param s The source node.
|
kpeter@2574
|
317 |
/// \param t The target node.
|
kpeter@2574
|
318 |
/// \param flow_value The required amount of flow from node \c s
|
kpeter@2574
|
319 |
/// to node \c t (i.e. the supply of \c s and the demand of \c t).
|
kpeter@2574
|
320 |
CapacityScaling( const Graph &graph,
|
kpeter@2574
|
321 |
const LowerMap &lower,
|
kpeter@2574
|
322 |
const CapacityMap &capacity,
|
kpeter@2574
|
323 |
const CostMap &cost,
|
kpeter@2574
|
324 |
Node s, Node t,
|
kpeter@2574
|
325 |
Supply flow_value ) :
|
kpeter@2574
|
326 |
_graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
|
kpeter@2574
|
327 |
_supply(graph), _flow(graph, 0), _potential(graph, 0),
|
kpeter@2574
|
328 |
_res_cap(graph), _excess(graph), _pred(graph),
|
kpeter@2574
|
329 |
_dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
|
deba@2440
|
330 |
{
|
kpeter@2556
|
331 |
// Removing non-zero lower bounds
|
kpeter@2574
|
332 |
_capacity = subMap(capacity, lower);
|
kpeter@2574
|
333 |
_res_cap = _capacity;
|
kpeter@2574
|
334 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@2574
|
335 |
Supply sum = 0;
|
kpeter@2574
|
336 |
if (n == s) sum = flow_value;
|
kpeter@2574
|
337 |
if (n == t) sum = -flow_value;
|
kpeter@2574
|
338 |
for (InEdgeIt e(_graph, n); e != INVALID; ++e)
|
kpeter@2574
|
339 |
sum += lower[e];
|
kpeter@2574
|
340 |
for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
|
kpeter@2574
|
341 |
sum -= lower[e];
|
kpeter@2574
|
342 |
_supply[n] = sum;
|
deba@2440
|
343 |
}
|
kpeter@2574
|
344 |
_valid_supply = true;
|
deba@2440
|
345 |
}
|
deba@2440
|
346 |
|
deba@2440
|
347 |
/// \brief Simple constructor of the class (without lower bounds).
|
deba@2440
|
348 |
///
|
deba@2440
|
349 |
/// Simple constructor of the class (without lower bounds).
|
deba@2440
|
350 |
///
|
kpeter@2574
|
351 |
/// \param graph The directed graph the algorithm runs on.
|
kpeter@2574
|
352 |
/// \param capacity The capacities (upper bounds) of the edges.
|
kpeter@2574
|
353 |
/// \param cost The cost (length) values of the edges.
|
kpeter@2574
|
354 |
/// \param s The source node.
|
kpeter@2574
|
355 |
/// \param t The target node.
|
kpeter@2574
|
356 |
/// \param flow_value The required amount of flow from node \c s
|
kpeter@2574
|
357 |
/// to node \c t (i.e. the supply of \c s and the demand of \c t).
|
kpeter@2574
|
358 |
CapacityScaling( const Graph &graph,
|
kpeter@2574
|
359 |
const CapacityMap &capacity,
|
kpeter@2574
|
360 |
const CostMap &cost,
|
kpeter@2574
|
361 |
Node s, Node t,
|
kpeter@2574
|
362 |
Supply flow_value ) :
|
kpeter@2574
|
363 |
_graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
|
kpeter@2574
|
364 |
_supply(graph, 0), _flow(graph, 0), _potential(graph, 0),
|
kpeter@2574
|
365 |
_res_cap(capacity), _excess(graph), _pred(graph),
|
kpeter@2574
|
366 |
_dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
|
deba@2440
|
367 |
{
|
kpeter@2574
|
368 |
_supply[s] = flow_value;
|
kpeter@2574
|
369 |
_supply[t] = -flow_value;
|
kpeter@2574
|
370 |
_valid_supply = true;
|
deba@2440
|
371 |
}
|
deba@2440
|
372 |
|
kpeter@2556
|
373 |
/// \brief Runs the algorithm.
|
kpeter@2556
|
374 |
///
|
kpeter@2556
|
375 |
/// Runs the algorithm.
|
kpeter@2556
|
376 |
///
|
kpeter@2574
|
377 |
/// \param scaling Enable or disable capacity scaling.
|
kpeter@2556
|
378 |
/// If the maximum edge capacity and/or the amount of total supply
|
kpeter@2574
|
379 |
/// is rather small, the algorithm could be slightly faster without
|
kpeter@2556
|
380 |
/// scaling.
|
kpeter@2556
|
381 |
///
|
kpeter@2556
|
382 |
/// \return \c true if a feasible flow can be found.
|
kpeter@2574
|
383 |
bool run(bool scaling = true) {
|
kpeter@2574
|
384 |
return init(scaling) && start();
|
kpeter@2556
|
385 |
}
|
kpeter@2556
|
386 |
|
kpeter@2574
|
387 |
/// \brief Returns a const reference to the edge map storing the
|
kpeter@2574
|
388 |
/// found flow.
|
deba@2440
|
389 |
///
|
kpeter@2574
|
390 |
/// Returns a const reference to the edge map storing the found flow.
|
deba@2440
|
391 |
///
|
deba@2440
|
392 |
/// \pre \ref run() must be called before using this function.
|
deba@2440
|
393 |
const FlowMap& flowMap() const {
|
kpeter@2574
|
394 |
return _flow;
|
deba@2440
|
395 |
}
|
deba@2440
|
396 |
|
kpeter@2574
|
397 |
/// \brief Returns a const reference to the node map storing the
|
kpeter@2574
|
398 |
/// found potentials (the dual solution).
|
deba@2440
|
399 |
///
|
kpeter@2574
|
400 |
/// Returns a const reference to the node map storing the found
|
kpeter@2574
|
401 |
/// potentials (the dual solution).
|
deba@2440
|
402 |
///
|
deba@2440
|
403 |
/// \pre \ref run() must be called before using this function.
|
deba@2440
|
404 |
const PotentialMap& potentialMap() const {
|
kpeter@2574
|
405 |
return _potential;
|
deba@2440
|
406 |
}
|
deba@2440
|
407 |
|
deba@2440
|
408 |
/// \brief Returns the total cost of the found flow.
|
deba@2440
|
409 |
///
|
deba@2440
|
410 |
/// Returns the total cost of the found flow. The complexity of the
|
deba@2440
|
411 |
/// function is \f$ O(e) \f$.
|
deba@2440
|
412 |
///
|
deba@2440
|
413 |
/// \pre \ref run() must be called before using this function.
|
deba@2440
|
414 |
Cost totalCost() const {
|
deba@2440
|
415 |
Cost c = 0;
|
kpeter@2574
|
416 |
for (EdgeIt e(_graph); e != INVALID; ++e)
|
kpeter@2574
|
417 |
c += _flow[e] * _cost[e];
|
deba@2440
|
418 |
return c;
|
deba@2440
|
419 |
}
|
deba@2440
|
420 |
|
kpeter@2574
|
421 |
private:
|
deba@2440
|
422 |
|
kpeter@2556
|
423 |
/// Initializes the algorithm.
|
kpeter@2574
|
424 |
bool init(bool scaling) {
|
kpeter@2574
|
425 |
if (!_valid_supply) return false;
|
kpeter@2574
|
426 |
_excess = _supply;
|
deba@2440
|
427 |
|
deba@2440
|
428 |
// Initilaizing delta value
|
kpeter@2574
|
429 |
if (scaling) {
|
kpeter@2535
|
430 |
// With scaling
|
kpeter@2535
|
431 |
Supply max_sup = 0, max_dem = 0;
|
kpeter@2574
|
432 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@2574
|
433 |
if ( _supply[n] > max_sup) max_sup = _supply[n];
|
kpeter@2574
|
434 |
if (-_supply[n] > max_dem) max_dem = -_supply[n];
|
kpeter@2535
|
435 |
}
|
kpeter@2535
|
436 |
if (max_dem < max_sup) max_sup = max_dem;
|
kpeter@2574
|
437 |
_phase_num = 0;
|
kpeter@2574
|
438 |
for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
|
kpeter@2574
|
439 |
++_phase_num;
|
kpeter@2535
|
440 |
} else {
|
kpeter@2535
|
441 |
// Without scaling
|
kpeter@2574
|
442 |
_delta = 1;
|
deba@2440
|
443 |
}
|
deba@2440
|
444 |
return true;
|
deba@2440
|
445 |
}
|
deba@2440
|
446 |
|
kpeter@2535
|
447 |
bool start() {
|
kpeter@2574
|
448 |
if (_delta > 1)
|
kpeter@2535
|
449 |
return startWithScaling();
|
kpeter@2535
|
450 |
else
|
kpeter@2535
|
451 |
return startWithoutScaling();
|
kpeter@2535
|
452 |
}
|
kpeter@2535
|
453 |
|
kpeter@2574
|
454 |
/// Executes the capacity scaling algorithm.
|
kpeter@2535
|
455 |
bool startWithScaling() {
|
kpeter@2535
|
456 |
// Processing capacity scaling phases
|
kpeter@2535
|
457 |
Node s, t;
|
kpeter@2535
|
458 |
int phase_cnt = 0;
|
kpeter@2535
|
459 |
int factor = 4;
|
kpeter@2535
|
460 |
while (true) {
|
kpeter@2535
|
461 |
// Saturating all edges not satisfying the optimality condition
|
kpeter@2574
|
462 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
kpeter@2574
|
463 |
Node u = _graph.source(e), v = _graph.target(e);
|
kpeter@2574
|
464 |
Cost c = _cost[e] + _potential[u] - _potential[v];
|
kpeter@2574
|
465 |
if (c < 0 && _res_cap[e] >= _delta) {
|
kpeter@2574
|
466 |
_excess[u] -= _res_cap[e];
|
kpeter@2574
|
467 |
_excess[v] += _res_cap[e];
|
kpeter@2574
|
468 |
_flow[e] = _capacity[e];
|
kpeter@2574
|
469 |
_res_cap[e] = 0;
|
kpeter@2535
|
470 |
}
|
kpeter@2574
|
471 |
else if (c > 0 && _flow[e] >= _delta) {
|
kpeter@2574
|
472 |
_excess[u] += _flow[e];
|
kpeter@2574
|
473 |
_excess[v] -= _flow[e];
|
kpeter@2574
|
474 |
_flow[e] = 0;
|
kpeter@2574
|
475 |
_res_cap[e] = _capacity[e];
|
kpeter@2535
|
476 |
}
|
kpeter@2535
|
477 |
}
|
kpeter@2535
|
478 |
|
kpeter@2535
|
479 |
// Finding excess nodes and deficit nodes
|
kpeter@2574
|
480 |
_excess_nodes.clear();
|
kpeter@2574
|
481 |
_deficit_nodes.clear();
|
kpeter@2574
|
482 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@2574
|
483 |
if (_excess[n] >= _delta) _excess_nodes.push_back(n);
|
kpeter@2574
|
484 |
if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
|
kpeter@2535
|
485 |
}
|
kpeter@2556
|
486 |
int next_node = 0;
|
kpeter@2535
|
487 |
|
kpeter@2535
|
488 |
// Finding augmenting shortest paths
|
kpeter@2574
|
489 |
while (next_node < int(_excess_nodes.size())) {
|
kpeter@2535
|
490 |
// Checking deficit nodes
|
kpeter@2574
|
491 |
if (_delta > 1) {
|
kpeter@2535
|
492 |
bool delta_deficit = false;
|
kpeter@2574
|
493 |
for (int i = 0; i < int(_deficit_nodes.size()); ++i) {
|
kpeter@2574
|
494 |
if (_excess[_deficit_nodes[i]] <= -_delta) {
|
kpeter@2535
|
495 |
delta_deficit = true;
|
kpeter@2535
|
496 |
break;
|
kpeter@2535
|
497 |
}
|
kpeter@2535
|
498 |
}
|
kpeter@2535
|
499 |
if (!delta_deficit) break;
|
kpeter@2535
|
500 |
}
|
kpeter@2535
|
501 |
|
kpeter@2535
|
502 |
// Running Dijkstra
|
kpeter@2574
|
503 |
s = _excess_nodes[next_node];
|
kpeter@2574
|
504 |
if ((t = _dijkstra.run(s, _delta)) == INVALID) {
|
kpeter@2574
|
505 |
if (_delta > 1) {
|
kpeter@2535
|
506 |
++next_node;
|
kpeter@2535
|
507 |
continue;
|
kpeter@2535
|
508 |
}
|
kpeter@2535
|
509 |
return false;
|
kpeter@2535
|
510 |
}
|
kpeter@2535
|
511 |
|
kpeter@2535
|
512 |
// Augmenting along a shortest path from s to t.
|
kpeter@2574
|
513 |
Capacity d = _excess[s] < -_excess[t] ? _excess[s] : -_excess[t];
|
kpeter@2535
|
514 |
Node u = t;
|
kpeter@2535
|
515 |
Edge e;
|
kpeter@2574
|
516 |
if (d > _delta) {
|
kpeter@2574
|
517 |
while ((e = _pred[u]) != INVALID) {
|
kpeter@2535
|
518 |
Capacity rc;
|
kpeter@2574
|
519 |
if (u == _graph.target(e)) {
|
kpeter@2574
|
520 |
rc = _res_cap[e];
|
kpeter@2574
|
521 |
u = _graph.source(e);
|
kpeter@2535
|
522 |
} else {
|
kpeter@2574
|
523 |
rc = _flow[e];
|
kpeter@2574
|
524 |
u = _graph.target(e);
|
kpeter@2535
|
525 |
}
|
kpeter@2535
|
526 |
if (rc < d) d = rc;
|
kpeter@2535
|
527 |
}
|
kpeter@2535
|
528 |
}
|
kpeter@2535
|
529 |
u = t;
|
kpeter@2574
|
530 |
while ((e = _pred[u]) != INVALID) {
|
kpeter@2574
|
531 |
if (u == _graph.target(e)) {
|
kpeter@2574
|
532 |
_flow[e] += d;
|
kpeter@2574
|
533 |
_res_cap[e] -= d;
|
kpeter@2574
|
534 |
u = _graph.source(e);
|
kpeter@2535
|
535 |
} else {
|
kpeter@2574
|
536 |
_flow[e] -= d;
|
kpeter@2574
|
537 |
_res_cap[e] += d;
|
kpeter@2574
|
538 |
u = _graph.target(e);
|
kpeter@2535
|
539 |
}
|
kpeter@2535
|
540 |
}
|
kpeter@2574
|
541 |
_excess[s] -= d;
|
kpeter@2574
|
542 |
_excess[t] += d;
|
kpeter@2535
|
543 |
|
kpeter@2574
|
544 |
if (_excess[s] < _delta) ++next_node;
|
kpeter@2535
|
545 |
}
|
kpeter@2535
|
546 |
|
kpeter@2574
|
547 |
if (_delta == 1) break;
|
kpeter@2574
|
548 |
if (++phase_cnt > _phase_num / 4) factor = 2;
|
kpeter@2574
|
549 |
_delta = _delta <= factor ? 1 : _delta / factor;
|
kpeter@2535
|
550 |
}
|
kpeter@2535
|
551 |
|
kpeter@2556
|
552 |
// Handling non-zero lower bounds
|
kpeter@2574
|
553 |
if (_lower) {
|
kpeter@2574
|
554 |
for (EdgeIt e(_graph); e != INVALID; ++e)
|
kpeter@2574
|
555 |
_flow[e] += (*_lower)[e];
|
kpeter@2535
|
556 |
}
|
kpeter@2535
|
557 |
return true;
|
kpeter@2535
|
558 |
}
|
kpeter@2535
|
559 |
|
kpeter@2574
|
560 |
/// Executes the successive shortest path algorithm.
|
kpeter@2535
|
561 |
bool startWithoutScaling() {
|
deba@2440
|
562 |
// Finding excess nodes
|
kpeter@2574
|
563 |
for (NodeIt n(_graph); n != INVALID; ++n)
|
kpeter@2574
|
564 |
if (_excess[n] > 0) _excess_nodes.push_back(n);
|
kpeter@2574
|
565 |
if (_excess_nodes.size() == 0) return true;
|
kpeter@2556
|
566 |
int next_node = 0;
|
deba@2440
|
567 |
|
deba@2457
|
568 |
// Finding shortest paths
|
kpeter@2535
|
569 |
Node s, t;
|
kpeter@2574
|
570 |
while ( _excess[_excess_nodes[next_node]] > 0 ||
|
kpeter@2574
|
571 |
++next_node < int(_excess_nodes.size()) )
|
deba@2440
|
572 |
{
|
kpeter@2535
|
573 |
// Running Dijkstra
|
kpeter@2574
|
574 |
s = _excess_nodes[next_node];
|
kpeter@2574
|
575 |
if ((t = _dijkstra.run(s, 1)) == INVALID)
|
kpeter@2535
|
576 |
return false;
|
deba@2440
|
577 |
|
kpeter@2535
|
578 |
// Augmenting along a shortest path from s to t
|
kpeter@2574
|
579 |
Capacity d = _excess[s] < -_excess[t] ? _excess[s] : -_excess[t];
|
kpeter@2535
|
580 |
Node u = t;
|
kpeter@2535
|
581 |
Edge e;
|
kpeter@2574
|
582 |
while ((e = _pred[u]) != INVALID) {
|
kpeter@2535
|
583 |
Capacity rc;
|
kpeter@2574
|
584 |
if (u == _graph.target(e)) {
|
kpeter@2574
|
585 |
rc = _res_cap[e];
|
kpeter@2574
|
586 |
u = _graph.source(e);
|
kpeter@2535
|
587 |
} else {
|
kpeter@2574
|
588 |
rc = _flow[e];
|
kpeter@2574
|
589 |
u = _graph.target(e);
|
kpeter@2535
|
590 |
}
|
kpeter@2535
|
591 |
if (rc < d) d = rc;
|
kpeter@2535
|
592 |
}
|
kpeter@2535
|
593 |
u = t;
|
kpeter@2574
|
594 |
while ((e = _pred[u]) != INVALID) {
|
kpeter@2574
|
595 |
if (u == _graph.target(e)) {
|
kpeter@2574
|
596 |
_flow[e] += d;
|
kpeter@2574
|
597 |
_res_cap[e] -= d;
|
kpeter@2574
|
598 |
u = _graph.source(e);
|
kpeter@2535
|
599 |
} else {
|
kpeter@2574
|
600 |
_flow[e] -= d;
|
kpeter@2574
|
601 |
_res_cap[e] += d;
|
kpeter@2574
|
602 |
u = _graph.target(e);
|
kpeter@2535
|
603 |
}
|
kpeter@2535
|
604 |
}
|
kpeter@2574
|
605 |
_excess[s] -= d;
|
kpeter@2574
|
606 |
_excess[t] += d;
|
deba@2440
|
607 |
}
|
deba@2440
|
608 |
|
kpeter@2556
|
609 |
// Handling non-zero lower bounds
|
kpeter@2574
|
610 |
if (_lower) {
|
kpeter@2574
|
611 |
for (EdgeIt e(_graph); e != INVALID; ++e)
|
kpeter@2574
|
612 |
_flow[e] += (*_lower)[e];
|
deba@2440
|
613 |
}
|
deba@2440
|
614 |
return true;
|
deba@2440
|
615 |
}
|
deba@2440
|
616 |
|
deba@2440
|
617 |
}; //class CapacityScaling
|
deba@2440
|
618 |
|
deba@2440
|
619 |
///@}
|
deba@2440
|
620 |
|
deba@2440
|
621 |
} //namespace lemon
|
deba@2440
|
622 |
|
deba@2440
|
623 |
#endif //LEMON_CAPACITY_SCALING_H
|