alpar@877
|
1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*-
|
kpeter@805
|
2 |
*
|
alpar@877
|
3 |
* This file is a part of LEMON, a generic C++ optimization library.
|
kpeter@805
|
4 |
*
|
alpar@877
|
5 |
* Copyright (C) 2003-2010
|
kpeter@805
|
6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
|
kpeter@805
|
7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES).
|
kpeter@805
|
8 |
*
|
kpeter@805
|
9 |
* Permission to use, modify and distribute this software is granted
|
kpeter@805
|
10 |
* provided that this copyright notice appears in all copies. For
|
kpeter@805
|
11 |
* precise terms see the accompanying LICENSE file.
|
kpeter@805
|
12 |
*
|
kpeter@805
|
13 |
* This software is provided "AS IS" with no warranty of any kind,
|
kpeter@805
|
14 |
* express or implied, and with no claim as to its suitability for any
|
kpeter@805
|
15 |
* purpose.
|
kpeter@805
|
16 |
*
|
kpeter@805
|
17 |
*/
|
kpeter@805
|
18 |
|
kpeter@805
|
19 |
#ifndef LEMON_CAPACITY_SCALING_H
|
kpeter@805
|
20 |
#define LEMON_CAPACITY_SCALING_H
|
kpeter@805
|
21 |
|
kpeter@806
|
22 |
/// \ingroup min_cost_flow_algs
|
kpeter@805
|
23 |
///
|
kpeter@805
|
24 |
/// \file
|
kpeter@806
|
25 |
/// \brief Capacity Scaling algorithm for finding a minimum cost flow.
|
kpeter@805
|
26 |
|
kpeter@805
|
27 |
#include <vector>
|
kpeter@806
|
28 |
#include <limits>
|
kpeter@806
|
29 |
#include <lemon/core.h>
|
kpeter@805
|
30 |
#include <lemon/bin_heap.h>
|
kpeter@805
|
31 |
|
kpeter@805
|
32 |
namespace lemon {
|
kpeter@805
|
33 |
|
kpeter@807
|
34 |
/// \brief Default traits class of CapacityScaling algorithm.
|
kpeter@807
|
35 |
///
|
kpeter@807
|
36 |
/// Default traits class of CapacityScaling algorithm.
|
kpeter@807
|
37 |
/// \tparam GR Digraph type.
|
kpeter@812
|
38 |
/// \tparam V The number type used for flow amounts, capacity bounds
|
kpeter@807
|
39 |
/// and supply values. By default it is \c int.
|
kpeter@812
|
40 |
/// \tparam C The number type used for costs and potentials.
|
kpeter@807
|
41 |
/// By default it is the same as \c V.
|
kpeter@807
|
42 |
template <typename GR, typename V = int, typename C = V>
|
kpeter@807
|
43 |
struct CapacityScalingDefaultTraits
|
kpeter@807
|
44 |
{
|
kpeter@807
|
45 |
/// The type of the digraph
|
kpeter@807
|
46 |
typedef GR Digraph;
|
kpeter@807
|
47 |
/// The type of the flow amounts, capacity bounds and supply values
|
kpeter@807
|
48 |
typedef V Value;
|
kpeter@807
|
49 |
/// The type of the arc costs
|
kpeter@807
|
50 |
typedef C Cost;
|
kpeter@807
|
51 |
|
kpeter@807
|
52 |
/// \brief The type of the heap used for internal Dijkstra computations.
|
kpeter@807
|
53 |
///
|
kpeter@807
|
54 |
/// The type of the heap used for internal Dijkstra computations.
|
kpeter@807
|
55 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
|
kpeter@807
|
56 |
/// its priority type must be \c Cost and its cross reference type
|
kpeter@807
|
57 |
/// must be \ref RangeMap "RangeMap<int>".
|
kpeter@807
|
58 |
typedef BinHeap<Cost, RangeMap<int> > Heap;
|
kpeter@807
|
59 |
};
|
kpeter@807
|
60 |
|
kpeter@806
|
61 |
/// \addtogroup min_cost_flow_algs
|
kpeter@805
|
62 |
/// @{
|
kpeter@805
|
63 |
|
kpeter@806
|
64 |
/// \brief Implementation of the Capacity Scaling algorithm for
|
kpeter@806
|
65 |
/// finding a \ref min_cost_flow "minimum cost flow".
|
kpeter@805
|
66 |
///
|
kpeter@805
|
67 |
/// \ref CapacityScaling implements the capacity scaling version
|
kpeter@806
|
68 |
/// of the successive shortest path algorithm for finding a
|
kpeter@813
|
69 |
/// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,
|
kpeter@813
|
70 |
/// \ref edmondskarp72theoretical. It is an efficient dual
|
kpeter@1049
|
71 |
/// solution method, which runs in polynomial time
|
kpeter@1049
|
72 |
/// \f$O(e\log U (n+e)\log n)\f$, where <i>U</i> denotes the maximum
|
kpeter@1049
|
73 |
/// of node supply and arc capacity values.
|
kpeter@805
|
74 |
///
|
kpeter@1003
|
75 |
/// This algorithm is typically slower than \ref CostScaling and
|
kpeter@1003
|
76 |
/// \ref NetworkSimplex, but in special cases, it can be more
|
kpeter@1003
|
77 |
/// efficient than them.
|
kpeter@1003
|
78 |
/// (For more information, see \ref min_cost_flow_algs "the module page".)
|
kpeter@1003
|
79 |
///
|
kpeter@806
|
80 |
/// Most of the parameters of the problem (except for the digraph)
|
kpeter@806
|
81 |
/// can be given using separate functions, and the algorithm can be
|
kpeter@806
|
82 |
/// executed using the \ref run() function. If some parameters are not
|
kpeter@806
|
83 |
/// specified, then default values will be used.
|
kpeter@805
|
84 |
///
|
kpeter@806
|
85 |
/// \tparam GR The digraph type the algorithm runs on.
|
kpeter@812
|
86 |
/// \tparam V The number type used for flow amounts, capacity bounds
|
kpeter@825
|
87 |
/// and supply values in the algorithm. By default, it is \c int.
|
kpeter@812
|
88 |
/// \tparam C The number type used for costs and potentials in the
|
kpeter@825
|
89 |
/// algorithm. By default, it is the same as \c V.
|
kpeter@825
|
90 |
/// \tparam TR The traits class that defines various types used by the
|
kpeter@825
|
91 |
/// algorithm. By default, it is \ref CapacityScalingDefaultTraits
|
kpeter@825
|
92 |
/// "CapacityScalingDefaultTraits<GR, V, C>".
|
kpeter@825
|
93 |
/// In most cases, this parameter should not be set directly,
|
kpeter@825
|
94 |
/// consider to use the named template parameters instead.
|
kpeter@805
|
95 |
///
|
kpeter@921
|
96 |
/// \warning Both \c V and \c C must be signed number types.
|
kpeter@985
|
97 |
/// \warning Capacity bounds and supply values must be integer, but
|
kpeter@985
|
98 |
/// arc costs can be arbitrary real numbers.
|
kpeter@919
|
99 |
/// \warning This algorithm does not support negative costs for
|
kpeter@919
|
100 |
/// arcs having infinite upper bound.
|
kpeter@807
|
101 |
#ifdef DOXYGEN
|
kpeter@807
|
102 |
template <typename GR, typename V, typename C, typename TR>
|
kpeter@807
|
103 |
#else
|
kpeter@807
|
104 |
template < typename GR, typename V = int, typename C = V,
|
kpeter@807
|
105 |
typename TR = CapacityScalingDefaultTraits<GR, V, C> >
|
kpeter@807
|
106 |
#endif
|
kpeter@805
|
107 |
class CapacityScaling
|
kpeter@805
|
108 |
{
|
kpeter@806
|
109 |
public:
|
kpeter@805
|
110 |
|
kpeter@807
|
111 |
/// The type of the digraph
|
kpeter@807
|
112 |
typedef typename TR::Digraph Digraph;
|
kpeter@806
|
113 |
/// The type of the flow amounts, capacity bounds and supply values
|
kpeter@807
|
114 |
typedef typename TR::Value Value;
|
kpeter@806
|
115 |
/// The type of the arc costs
|
kpeter@807
|
116 |
typedef typename TR::Cost Cost;
|
kpeter@807
|
117 |
|
kpeter@807
|
118 |
/// The type of the heap used for internal Dijkstra computations
|
kpeter@807
|
119 |
typedef typename TR::Heap Heap;
|
kpeter@807
|
120 |
|
kpeter@807
|
121 |
/// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
|
kpeter@807
|
122 |
typedef TR Traits;
|
kpeter@805
|
123 |
|
kpeter@805
|
124 |
public:
|
kpeter@805
|
125 |
|
kpeter@806
|
126 |
/// \brief Problem type constants for the \c run() function.
|
kpeter@806
|
127 |
///
|
kpeter@806
|
128 |
/// Enum type containing the problem type constants that can be
|
kpeter@806
|
129 |
/// returned by the \ref run() function of the algorithm.
|
kpeter@806
|
130 |
enum ProblemType {
|
kpeter@806
|
131 |
/// The problem has no feasible solution (flow).
|
kpeter@806
|
132 |
INFEASIBLE,
|
kpeter@806
|
133 |
/// The problem has optimal solution (i.e. it is feasible and
|
kpeter@806
|
134 |
/// bounded), and the algorithm has found optimal flow and node
|
kpeter@806
|
135 |
/// potentials (primal and dual solutions).
|
kpeter@806
|
136 |
OPTIMAL,
|
kpeter@806
|
137 |
/// The digraph contains an arc of negative cost and infinite
|
kpeter@806
|
138 |
/// upper bound. It means that the objective function is unbounded
|
kpeter@812
|
139 |
/// on that arc, however, note that it could actually be bounded
|
kpeter@806
|
140 |
/// over the feasible flows, but this algroithm cannot handle
|
kpeter@806
|
141 |
/// these cases.
|
kpeter@806
|
142 |
UNBOUNDED
|
kpeter@806
|
143 |
};
|
alpar@877
|
144 |
|
kpeter@806
|
145 |
private:
|
kpeter@806
|
146 |
|
kpeter@806
|
147 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR);
|
kpeter@806
|
148 |
|
kpeter@806
|
149 |
typedef std::vector<int> IntVector;
|
kpeter@806
|
150 |
typedef std::vector<Value> ValueVector;
|
kpeter@806
|
151 |
typedef std::vector<Cost> CostVector;
|
kpeter@839
|
152 |
typedef std::vector<char> BoolVector;
|
kpeter@839
|
153 |
// Note: vector<char> is used instead of vector<bool> for efficiency reasons
|
kpeter@805
|
154 |
|
kpeter@805
|
155 |
private:
|
kpeter@805
|
156 |
|
kpeter@806
|
157 |
// Data related to the underlying digraph
|
kpeter@806
|
158 |
const GR &_graph;
|
kpeter@806
|
159 |
int _node_num;
|
kpeter@806
|
160 |
int _arc_num;
|
kpeter@806
|
161 |
int _res_arc_num;
|
kpeter@806
|
162 |
int _root;
|
kpeter@806
|
163 |
|
kpeter@806
|
164 |
// Parameters of the problem
|
kpeter@806
|
165 |
bool _have_lower;
|
kpeter@806
|
166 |
Value _sum_supply;
|
kpeter@806
|
167 |
|
kpeter@806
|
168 |
// Data structures for storing the digraph
|
kpeter@806
|
169 |
IntNodeMap _node_id;
|
kpeter@806
|
170 |
IntArcMap _arc_idf;
|
kpeter@806
|
171 |
IntArcMap _arc_idb;
|
kpeter@806
|
172 |
IntVector _first_out;
|
kpeter@806
|
173 |
BoolVector _forward;
|
kpeter@806
|
174 |
IntVector _source;
|
kpeter@806
|
175 |
IntVector _target;
|
kpeter@806
|
176 |
IntVector _reverse;
|
kpeter@806
|
177 |
|
kpeter@806
|
178 |
// Node and arc data
|
kpeter@806
|
179 |
ValueVector _lower;
|
kpeter@806
|
180 |
ValueVector _upper;
|
kpeter@806
|
181 |
CostVector _cost;
|
kpeter@806
|
182 |
ValueVector _supply;
|
kpeter@806
|
183 |
|
kpeter@806
|
184 |
ValueVector _res_cap;
|
kpeter@806
|
185 |
CostVector _pi;
|
kpeter@806
|
186 |
ValueVector _excess;
|
kpeter@806
|
187 |
IntVector _excess_nodes;
|
kpeter@806
|
188 |
IntVector _deficit_nodes;
|
kpeter@806
|
189 |
|
kpeter@806
|
190 |
Value _delta;
|
kpeter@810
|
191 |
int _factor;
|
kpeter@806
|
192 |
IntVector _pred;
|
kpeter@806
|
193 |
|
kpeter@806
|
194 |
public:
|
alpar@877
|
195 |
|
kpeter@806
|
196 |
/// \brief Constant for infinite upper bounds (capacities).
|
kpeter@805
|
197 |
///
|
kpeter@806
|
198 |
/// Constant for infinite upper bounds (capacities).
|
kpeter@806
|
199 |
/// It is \c std::numeric_limits<Value>::infinity() if available,
|
kpeter@806
|
200 |
/// \c std::numeric_limits<Value>::max() otherwise.
|
kpeter@806
|
201 |
const Value INF;
|
kpeter@806
|
202 |
|
kpeter@806
|
203 |
private:
|
kpeter@806
|
204 |
|
kpeter@806
|
205 |
// Special implementation of the Dijkstra algorithm for finding
|
kpeter@806
|
206 |
// shortest paths in the residual network of the digraph with
|
kpeter@806
|
207 |
// respect to the reduced arc costs and modifying the node
|
kpeter@806
|
208 |
// potentials according to the found distance labels.
|
kpeter@805
|
209 |
class ResidualDijkstra
|
kpeter@805
|
210 |
{
|
kpeter@805
|
211 |
private:
|
kpeter@805
|
212 |
|
kpeter@806
|
213 |
int _node_num;
|
kpeter@811
|
214 |
bool _geq;
|
kpeter@806
|
215 |
const IntVector &_first_out;
|
kpeter@806
|
216 |
const IntVector &_target;
|
kpeter@806
|
217 |
const CostVector &_cost;
|
kpeter@806
|
218 |
const ValueVector &_res_cap;
|
kpeter@806
|
219 |
const ValueVector &_excess;
|
kpeter@806
|
220 |
CostVector &_pi;
|
kpeter@806
|
221 |
IntVector &_pred;
|
alpar@877
|
222 |
|
kpeter@806
|
223 |
IntVector _proc_nodes;
|
kpeter@806
|
224 |
CostVector _dist;
|
alpar@877
|
225 |
|
kpeter@805
|
226 |
public:
|
kpeter@805
|
227 |
|
kpeter@806
|
228 |
ResidualDijkstra(CapacityScaling& cs) :
|
kpeter@811
|
229 |
_node_num(cs._node_num), _geq(cs._sum_supply < 0),
|
kpeter@811
|
230 |
_first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
|
kpeter@811
|
231 |
_res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
|
kpeter@811
|
232 |
_pred(cs._pred), _dist(cs._node_num)
|
kpeter@805
|
233 |
{}
|
kpeter@805
|
234 |
|
kpeter@806
|
235 |
int run(int s, Value delta = 1) {
|
kpeter@807
|
236 |
RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
|
kpeter@805
|
237 |
Heap heap(heap_cross_ref);
|
kpeter@805
|
238 |
heap.push(s, 0);
|
kpeter@806
|
239 |
_pred[s] = -1;
|
kpeter@805
|
240 |
_proc_nodes.clear();
|
kpeter@805
|
241 |
|
kpeter@806
|
242 |
// Process nodes
|
kpeter@805
|
243 |
while (!heap.empty() && _excess[heap.top()] > -delta) {
|
kpeter@806
|
244 |
int u = heap.top(), v;
|
kpeter@806
|
245 |
Cost d = heap.prio() + _pi[u], dn;
|
kpeter@805
|
246 |
_dist[u] = heap.prio();
|
kpeter@806
|
247 |
_proc_nodes.push_back(u);
|
kpeter@805
|
248 |
heap.pop();
|
kpeter@805
|
249 |
|
kpeter@806
|
250 |
// Traverse outgoing residual arcs
|
kpeter@811
|
251 |
int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
|
kpeter@811
|
252 |
for (int a = _first_out[u]; a != last_out; ++a) {
|
kpeter@806
|
253 |
if (_res_cap[a] < delta) continue;
|
kpeter@806
|
254 |
v = _target[a];
|
kpeter@806
|
255 |
switch (heap.state(v)) {
|
kpeter@805
|
256 |
case Heap::PRE_HEAP:
|
kpeter@806
|
257 |
heap.push(v, d + _cost[a] - _pi[v]);
|
kpeter@806
|
258 |
_pred[v] = a;
|
kpeter@805
|
259 |
break;
|
kpeter@805
|
260 |
case Heap::IN_HEAP:
|
kpeter@806
|
261 |
dn = d + _cost[a] - _pi[v];
|
kpeter@806
|
262 |
if (dn < heap[v]) {
|
kpeter@806
|
263 |
heap.decrease(v, dn);
|
kpeter@806
|
264 |
_pred[v] = a;
|
kpeter@805
|
265 |
}
|
kpeter@805
|
266 |
break;
|
kpeter@805
|
267 |
case Heap::POST_HEAP:
|
kpeter@805
|
268 |
break;
|
kpeter@805
|
269 |
}
|
kpeter@805
|
270 |
}
|
kpeter@805
|
271 |
}
|
kpeter@806
|
272 |
if (heap.empty()) return -1;
|
kpeter@805
|
273 |
|
kpeter@806
|
274 |
// Update potentials of processed nodes
|
kpeter@806
|
275 |
int t = heap.top();
|
kpeter@806
|
276 |
Cost dt = heap.prio();
|
kpeter@806
|
277 |
for (int i = 0; i < int(_proc_nodes.size()); ++i) {
|
kpeter@806
|
278 |
_pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
|
kpeter@806
|
279 |
}
|
kpeter@805
|
280 |
|
kpeter@805
|
281 |
return t;
|
kpeter@805
|
282 |
}
|
kpeter@805
|
283 |
|
kpeter@805
|
284 |
}; //class ResidualDijkstra
|
kpeter@805
|
285 |
|
kpeter@805
|
286 |
public:
|
kpeter@805
|
287 |
|
kpeter@807
|
288 |
/// \name Named Template Parameters
|
kpeter@807
|
289 |
/// @{
|
kpeter@807
|
290 |
|
kpeter@807
|
291 |
template <typename T>
|
kpeter@807
|
292 |
struct SetHeapTraits : public Traits {
|
kpeter@807
|
293 |
typedef T Heap;
|
kpeter@807
|
294 |
};
|
kpeter@807
|
295 |
|
kpeter@807
|
296 |
/// \brief \ref named-templ-param "Named parameter" for setting
|
kpeter@807
|
297 |
/// \c Heap type.
|
kpeter@807
|
298 |
///
|
kpeter@807
|
299 |
/// \ref named-templ-param "Named parameter" for setting \c Heap
|
kpeter@807
|
300 |
/// type, which is used for internal Dijkstra computations.
|
kpeter@807
|
301 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
|
kpeter@807
|
302 |
/// its priority type must be \c Cost and its cross reference type
|
kpeter@807
|
303 |
/// must be \ref RangeMap "RangeMap<int>".
|
kpeter@807
|
304 |
template <typename T>
|
kpeter@807
|
305 |
struct SetHeap
|
kpeter@807
|
306 |
: public CapacityScaling<GR, V, C, SetHeapTraits<T> > {
|
kpeter@807
|
307 |
typedef CapacityScaling<GR, V, C, SetHeapTraits<T> > Create;
|
kpeter@807
|
308 |
};
|
kpeter@807
|
309 |
|
kpeter@807
|
310 |
/// @}
|
kpeter@807
|
311 |
|
kpeter@863
|
312 |
protected:
|
kpeter@863
|
313 |
|
kpeter@863
|
314 |
CapacityScaling() {}
|
kpeter@863
|
315 |
|
kpeter@807
|
316 |
public:
|
kpeter@807
|
317 |
|
kpeter@806
|
318 |
/// \brief Constructor.
|
kpeter@805
|
319 |
///
|
kpeter@806
|
320 |
/// The constructor of the class.
|
kpeter@805
|
321 |
///
|
kpeter@806
|
322 |
/// \param graph The digraph the algorithm runs on.
|
kpeter@806
|
323 |
CapacityScaling(const GR& graph) :
|
kpeter@806
|
324 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
|
kpeter@806
|
325 |
INF(std::numeric_limits<Value>::has_infinity ?
|
kpeter@806
|
326 |
std::numeric_limits<Value>::infinity() :
|
kpeter@806
|
327 |
std::numeric_limits<Value>::max())
|
kpeter@805
|
328 |
{
|
kpeter@812
|
329 |
// Check the number types
|
kpeter@806
|
330 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
|
kpeter@806
|
331 |
"The flow type of CapacityScaling must be signed");
|
kpeter@806
|
332 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
|
kpeter@806
|
333 |
"The cost type of CapacityScaling must be signed");
|
kpeter@806
|
334 |
|
kpeter@830
|
335 |
// Reset data structures
|
kpeter@806
|
336 |
reset();
|
kpeter@805
|
337 |
}
|
kpeter@805
|
338 |
|
kpeter@806
|
339 |
/// \name Parameters
|
kpeter@806
|
340 |
/// The parameters of the algorithm can be specified using these
|
kpeter@806
|
341 |
/// functions.
|
kpeter@806
|
342 |
|
kpeter@806
|
343 |
/// @{
|
kpeter@806
|
344 |
|
kpeter@806
|
345 |
/// \brief Set the lower bounds on the arcs.
|
kpeter@805
|
346 |
///
|
kpeter@806
|
347 |
/// This function sets the lower bounds on the arcs.
|
kpeter@806
|
348 |
/// If it is not used before calling \ref run(), the lower bounds
|
kpeter@806
|
349 |
/// will be set to zero on all arcs.
|
kpeter@805
|
350 |
///
|
kpeter@806
|
351 |
/// \param map An arc map storing the lower bounds.
|
kpeter@806
|
352 |
/// Its \c Value type must be convertible to the \c Value type
|
kpeter@806
|
353 |
/// of the algorithm.
|
kpeter@806
|
354 |
///
|
kpeter@806
|
355 |
/// \return <tt>(*this)</tt>
|
kpeter@806
|
356 |
template <typename LowerMap>
|
kpeter@806
|
357 |
CapacityScaling& lowerMap(const LowerMap& map) {
|
kpeter@806
|
358 |
_have_lower = true;
|
kpeter@806
|
359 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
kpeter@806
|
360 |
_lower[_arc_idf[a]] = map[a];
|
kpeter@806
|
361 |
_lower[_arc_idb[a]] = map[a];
|
kpeter@805
|
362 |
}
|
kpeter@805
|
363 |
return *this;
|
kpeter@805
|
364 |
}
|
kpeter@805
|
365 |
|
kpeter@806
|
366 |
/// \brief Set the upper bounds (capacities) on the arcs.
|
kpeter@805
|
367 |
///
|
kpeter@806
|
368 |
/// This function sets the upper bounds (capacities) on the arcs.
|
kpeter@806
|
369 |
/// If it is not used before calling \ref run(), the upper bounds
|
kpeter@806
|
370 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be
|
kpeter@812
|
371 |
/// unbounded from above).
|
kpeter@805
|
372 |
///
|
kpeter@806
|
373 |
/// \param map An arc map storing the upper bounds.
|
kpeter@806
|
374 |
/// Its \c Value type must be convertible to the \c Value type
|
kpeter@806
|
375 |
/// of the algorithm.
|
kpeter@806
|
376 |
///
|
kpeter@806
|
377 |
/// \return <tt>(*this)</tt>
|
kpeter@806
|
378 |
template<typename UpperMap>
|
kpeter@806
|
379 |
CapacityScaling& upperMap(const UpperMap& map) {
|
kpeter@806
|
380 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
kpeter@806
|
381 |
_upper[_arc_idf[a]] = map[a];
|
kpeter@805
|
382 |
}
|
kpeter@805
|
383 |
return *this;
|
kpeter@805
|
384 |
}
|
kpeter@805
|
385 |
|
kpeter@806
|
386 |
/// \brief Set the costs of the arcs.
|
kpeter@806
|
387 |
///
|
kpeter@806
|
388 |
/// This function sets the costs of the arcs.
|
kpeter@806
|
389 |
/// If it is not used before calling \ref run(), the costs
|
kpeter@806
|
390 |
/// will be set to \c 1 on all arcs.
|
kpeter@806
|
391 |
///
|
kpeter@806
|
392 |
/// \param map An arc map storing the costs.
|
kpeter@806
|
393 |
/// Its \c Value type must be convertible to the \c Cost type
|
kpeter@806
|
394 |
/// of the algorithm.
|
kpeter@806
|
395 |
///
|
kpeter@806
|
396 |
/// \return <tt>(*this)</tt>
|
kpeter@806
|
397 |
template<typename CostMap>
|
kpeter@806
|
398 |
CapacityScaling& costMap(const CostMap& map) {
|
kpeter@806
|
399 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
kpeter@806
|
400 |
_cost[_arc_idf[a]] = map[a];
|
kpeter@806
|
401 |
_cost[_arc_idb[a]] = -map[a];
|
kpeter@806
|
402 |
}
|
kpeter@806
|
403 |
return *this;
|
kpeter@806
|
404 |
}
|
kpeter@806
|
405 |
|
kpeter@806
|
406 |
/// \brief Set the supply values of the nodes.
|
kpeter@806
|
407 |
///
|
kpeter@806
|
408 |
/// This function sets the supply values of the nodes.
|
kpeter@806
|
409 |
/// If neither this function nor \ref stSupply() is used before
|
kpeter@806
|
410 |
/// calling \ref run(), the supply of each node will be set to zero.
|
kpeter@806
|
411 |
///
|
kpeter@806
|
412 |
/// \param map A node map storing the supply values.
|
kpeter@806
|
413 |
/// Its \c Value type must be convertible to the \c Value type
|
kpeter@806
|
414 |
/// of the algorithm.
|
kpeter@806
|
415 |
///
|
kpeter@806
|
416 |
/// \return <tt>(*this)</tt>
|
kpeter@806
|
417 |
template<typename SupplyMap>
|
kpeter@806
|
418 |
CapacityScaling& supplyMap(const SupplyMap& map) {
|
kpeter@806
|
419 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@806
|
420 |
_supply[_node_id[n]] = map[n];
|
kpeter@806
|
421 |
}
|
kpeter@806
|
422 |
return *this;
|
kpeter@806
|
423 |
}
|
kpeter@806
|
424 |
|
kpeter@806
|
425 |
/// \brief Set single source and target nodes and a supply value.
|
kpeter@806
|
426 |
///
|
kpeter@806
|
427 |
/// This function sets a single source node and a single target node
|
kpeter@806
|
428 |
/// and the required flow value.
|
kpeter@806
|
429 |
/// If neither this function nor \ref supplyMap() is used before
|
kpeter@806
|
430 |
/// calling \ref run(), the supply of each node will be set to zero.
|
kpeter@806
|
431 |
///
|
kpeter@806
|
432 |
/// Using this function has the same effect as using \ref supplyMap()
|
kpeter@919
|
433 |
/// with a map in which \c k is assigned to \c s, \c -k is
|
kpeter@806
|
434 |
/// assigned to \c t and all other nodes have zero supply value.
|
kpeter@806
|
435 |
///
|
kpeter@806
|
436 |
/// \param s The source node.
|
kpeter@806
|
437 |
/// \param t The target node.
|
kpeter@806
|
438 |
/// \param k The required amount of flow from node \c s to node \c t
|
kpeter@806
|
439 |
/// (i.e. the supply of \c s and the demand of \c t).
|
kpeter@806
|
440 |
///
|
kpeter@806
|
441 |
/// \return <tt>(*this)</tt>
|
kpeter@806
|
442 |
CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
|
kpeter@806
|
443 |
for (int i = 0; i != _node_num; ++i) {
|
kpeter@806
|
444 |
_supply[i] = 0;
|
kpeter@806
|
445 |
}
|
kpeter@806
|
446 |
_supply[_node_id[s]] = k;
|
kpeter@806
|
447 |
_supply[_node_id[t]] = -k;
|
kpeter@806
|
448 |
return *this;
|
kpeter@806
|
449 |
}
|
alpar@877
|
450 |
|
kpeter@806
|
451 |
/// @}
|
kpeter@806
|
452 |
|
kpeter@805
|
453 |
/// \name Execution control
|
kpeter@807
|
454 |
/// The algorithm can be executed using \ref run().
|
kpeter@805
|
455 |
|
kpeter@805
|
456 |
/// @{
|
kpeter@805
|
457 |
|
kpeter@805
|
458 |
/// \brief Run the algorithm.
|
kpeter@805
|
459 |
///
|
kpeter@805
|
460 |
/// This function runs the algorithm.
|
kpeter@806
|
461 |
/// The paramters can be specified using functions \ref lowerMap(),
|
kpeter@806
|
462 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
|
kpeter@806
|
463 |
/// For example,
|
kpeter@806
|
464 |
/// \code
|
kpeter@806
|
465 |
/// CapacityScaling<ListDigraph> cs(graph);
|
kpeter@806
|
466 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost)
|
kpeter@806
|
467 |
/// .supplyMap(sup).run();
|
kpeter@806
|
468 |
/// \endcode
|
kpeter@806
|
469 |
///
|
kpeter@830
|
470 |
/// This function can be called more than once. All the given parameters
|
kpeter@830
|
471 |
/// are kept for the next call, unless \ref resetParams() or \ref reset()
|
kpeter@830
|
472 |
/// is used, thus only the modified parameters have to be set again.
|
kpeter@830
|
473 |
/// If the underlying digraph was also modified after the construction
|
kpeter@830
|
474 |
/// of the class (or the last \ref reset() call), then the \ref reset()
|
kpeter@830
|
475 |
/// function must be called.
|
kpeter@805
|
476 |
///
|
kpeter@810
|
477 |
/// \param factor The capacity scaling factor. It must be larger than
|
kpeter@810
|
478 |
/// one to use scaling. If it is less or equal to one, then scaling
|
kpeter@810
|
479 |
/// will be disabled.
|
kpeter@805
|
480 |
///
|
kpeter@806
|
481 |
/// \return \c INFEASIBLE if no feasible flow exists,
|
kpeter@806
|
482 |
/// \n \c OPTIMAL if the problem has optimal solution
|
kpeter@806
|
483 |
/// (i.e. it is feasible and bounded), and the algorithm has found
|
kpeter@806
|
484 |
/// optimal flow and node potentials (primal and dual solutions),
|
kpeter@806
|
485 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost
|
kpeter@806
|
486 |
/// and infinite upper bound. It means that the objective function
|
kpeter@812
|
487 |
/// is unbounded on that arc, however, note that it could actually be
|
kpeter@806
|
488 |
/// bounded over the feasible flows, but this algroithm cannot handle
|
kpeter@806
|
489 |
/// these cases.
|
kpeter@806
|
490 |
///
|
kpeter@806
|
491 |
/// \see ProblemType
|
kpeter@830
|
492 |
/// \see resetParams(), reset()
|
kpeter@810
|
493 |
ProblemType run(int factor = 4) {
|
kpeter@810
|
494 |
_factor = factor;
|
kpeter@810
|
495 |
ProblemType pt = init();
|
kpeter@806
|
496 |
if (pt != OPTIMAL) return pt;
|
kpeter@806
|
497 |
return start();
|
kpeter@806
|
498 |
}
|
kpeter@806
|
499 |
|
kpeter@806
|
500 |
/// \brief Reset all the parameters that have been given before.
|
kpeter@806
|
501 |
///
|
kpeter@806
|
502 |
/// This function resets all the paramaters that have been given
|
kpeter@806
|
503 |
/// before using functions \ref lowerMap(), \ref upperMap(),
|
kpeter@806
|
504 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply().
|
kpeter@806
|
505 |
///
|
kpeter@830
|
506 |
/// It is useful for multiple \ref run() calls. Basically, all the given
|
kpeter@830
|
507 |
/// parameters are kept for the next \ref run() call, unless
|
kpeter@830
|
508 |
/// \ref resetParams() or \ref reset() is used.
|
kpeter@830
|
509 |
/// If the underlying digraph was also modified after the construction
|
kpeter@830
|
510 |
/// of the class or the last \ref reset() call, then the \ref reset()
|
kpeter@830
|
511 |
/// function must be used, otherwise \ref resetParams() is sufficient.
|
kpeter@806
|
512 |
///
|
kpeter@806
|
513 |
/// For example,
|
kpeter@806
|
514 |
/// \code
|
kpeter@806
|
515 |
/// CapacityScaling<ListDigraph> cs(graph);
|
kpeter@806
|
516 |
///
|
kpeter@806
|
517 |
/// // First run
|
kpeter@806
|
518 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost)
|
kpeter@806
|
519 |
/// .supplyMap(sup).run();
|
kpeter@806
|
520 |
///
|
kpeter@830
|
521 |
/// // Run again with modified cost map (resetParams() is not called,
|
kpeter@806
|
522 |
/// // so only the cost map have to be set again)
|
kpeter@806
|
523 |
/// cost[e] += 100;
|
kpeter@806
|
524 |
/// cs.costMap(cost).run();
|
kpeter@806
|
525 |
///
|
kpeter@830
|
526 |
/// // Run again from scratch using resetParams()
|
kpeter@806
|
527 |
/// // (the lower bounds will be set to zero on all arcs)
|
kpeter@830
|
528 |
/// cs.resetParams();
|
kpeter@806
|
529 |
/// cs.upperMap(capacity).costMap(cost)
|
kpeter@806
|
530 |
/// .supplyMap(sup).run();
|
kpeter@806
|
531 |
/// \endcode
|
kpeter@806
|
532 |
///
|
kpeter@806
|
533 |
/// \return <tt>(*this)</tt>
|
kpeter@830
|
534 |
///
|
kpeter@830
|
535 |
/// \see reset(), run()
|
kpeter@830
|
536 |
CapacityScaling& resetParams() {
|
kpeter@806
|
537 |
for (int i = 0; i != _node_num; ++i) {
|
kpeter@806
|
538 |
_supply[i] = 0;
|
kpeter@806
|
539 |
}
|
kpeter@806
|
540 |
for (int j = 0; j != _res_arc_num; ++j) {
|
kpeter@806
|
541 |
_lower[j] = 0;
|
kpeter@806
|
542 |
_upper[j] = INF;
|
kpeter@806
|
543 |
_cost[j] = _forward[j] ? 1 : -1;
|
kpeter@806
|
544 |
}
|
kpeter@806
|
545 |
_have_lower = false;
|
kpeter@806
|
546 |
return *this;
|
kpeter@805
|
547 |
}
|
kpeter@805
|
548 |
|
kpeter@830
|
549 |
/// \brief Reset the internal data structures and all the parameters
|
kpeter@830
|
550 |
/// that have been given before.
|
kpeter@830
|
551 |
///
|
kpeter@830
|
552 |
/// This function resets the internal data structures and all the
|
kpeter@830
|
553 |
/// paramaters that have been given before using functions \ref lowerMap(),
|
kpeter@830
|
554 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
|
kpeter@830
|
555 |
///
|
kpeter@830
|
556 |
/// It is useful for multiple \ref run() calls. Basically, all the given
|
kpeter@830
|
557 |
/// parameters are kept for the next \ref run() call, unless
|
kpeter@830
|
558 |
/// \ref resetParams() or \ref reset() is used.
|
kpeter@830
|
559 |
/// If the underlying digraph was also modified after the construction
|
kpeter@830
|
560 |
/// of the class or the last \ref reset() call, then the \ref reset()
|
kpeter@830
|
561 |
/// function must be used, otherwise \ref resetParams() is sufficient.
|
kpeter@830
|
562 |
///
|
kpeter@830
|
563 |
/// See \ref resetParams() for examples.
|
kpeter@830
|
564 |
///
|
kpeter@830
|
565 |
/// \return <tt>(*this)</tt>
|
kpeter@830
|
566 |
///
|
kpeter@830
|
567 |
/// \see resetParams(), run()
|
kpeter@830
|
568 |
CapacityScaling& reset() {
|
kpeter@830
|
569 |
// Resize vectors
|
kpeter@830
|
570 |
_node_num = countNodes(_graph);
|
kpeter@830
|
571 |
_arc_num = countArcs(_graph);
|
kpeter@830
|
572 |
_res_arc_num = 2 * (_arc_num + _node_num);
|
kpeter@830
|
573 |
_root = _node_num;
|
kpeter@830
|
574 |
++_node_num;
|
kpeter@830
|
575 |
|
kpeter@830
|
576 |
_first_out.resize(_node_num + 1);
|
kpeter@830
|
577 |
_forward.resize(_res_arc_num);
|
kpeter@830
|
578 |
_source.resize(_res_arc_num);
|
kpeter@830
|
579 |
_target.resize(_res_arc_num);
|
kpeter@830
|
580 |
_reverse.resize(_res_arc_num);
|
kpeter@830
|
581 |
|
kpeter@830
|
582 |
_lower.resize(_res_arc_num);
|
kpeter@830
|
583 |
_upper.resize(_res_arc_num);
|
kpeter@830
|
584 |
_cost.resize(_res_arc_num);
|
kpeter@830
|
585 |
_supply.resize(_node_num);
|
alpar@877
|
586 |
|
kpeter@830
|
587 |
_res_cap.resize(_res_arc_num);
|
kpeter@830
|
588 |
_pi.resize(_node_num);
|
kpeter@830
|
589 |
_excess.resize(_node_num);
|
kpeter@830
|
590 |
_pred.resize(_node_num);
|
kpeter@830
|
591 |
|
kpeter@830
|
592 |
// Copy the graph
|
kpeter@830
|
593 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
|
kpeter@830
|
594 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
|
kpeter@830
|
595 |
_node_id[n] = i;
|
kpeter@830
|
596 |
}
|
kpeter@830
|
597 |
i = 0;
|
kpeter@830
|
598 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
|
kpeter@830
|
599 |
_first_out[i] = j;
|
kpeter@830
|
600 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
|
kpeter@830
|
601 |
_arc_idf[a] = j;
|
kpeter@830
|
602 |
_forward[j] = true;
|
kpeter@830
|
603 |
_source[j] = i;
|
kpeter@830
|
604 |
_target[j] = _node_id[_graph.runningNode(a)];
|
kpeter@830
|
605 |
}
|
kpeter@830
|
606 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
|
kpeter@830
|
607 |
_arc_idb[a] = j;
|
kpeter@830
|
608 |
_forward[j] = false;
|
kpeter@830
|
609 |
_source[j] = i;
|
kpeter@830
|
610 |
_target[j] = _node_id[_graph.runningNode(a)];
|
kpeter@830
|
611 |
}
|
kpeter@830
|
612 |
_forward[j] = false;
|
kpeter@830
|
613 |
_source[j] = i;
|
kpeter@830
|
614 |
_target[j] = _root;
|
kpeter@830
|
615 |
_reverse[j] = k;
|
kpeter@830
|
616 |
_forward[k] = true;
|
kpeter@830
|
617 |
_source[k] = _root;
|
kpeter@830
|
618 |
_target[k] = i;
|
kpeter@830
|
619 |
_reverse[k] = j;
|
kpeter@830
|
620 |
++j; ++k;
|
kpeter@830
|
621 |
}
|
kpeter@830
|
622 |
_first_out[i] = j;
|
kpeter@830
|
623 |
_first_out[_node_num] = k;
|
kpeter@830
|
624 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
kpeter@830
|
625 |
int fi = _arc_idf[a];
|
kpeter@830
|
626 |
int bi = _arc_idb[a];
|
kpeter@830
|
627 |
_reverse[fi] = bi;
|
kpeter@830
|
628 |
_reverse[bi] = fi;
|
kpeter@830
|
629 |
}
|
alpar@877
|
630 |
|
kpeter@830
|
631 |
// Reset parameters
|
kpeter@830
|
632 |
resetParams();
|
kpeter@830
|
633 |
return *this;
|
kpeter@830
|
634 |
}
|
kpeter@830
|
635 |
|
kpeter@805
|
636 |
/// @}
|
kpeter@805
|
637 |
|
kpeter@805
|
638 |
/// \name Query Functions
|
kpeter@805
|
639 |
/// The results of the algorithm can be obtained using these
|
kpeter@805
|
640 |
/// functions.\n
|
kpeter@806
|
641 |
/// The \ref run() function must be called before using them.
|
kpeter@805
|
642 |
|
kpeter@805
|
643 |
/// @{
|
kpeter@805
|
644 |
|
kpeter@806
|
645 |
/// \brief Return the total cost of the found flow.
|
kpeter@805
|
646 |
///
|
kpeter@806
|
647 |
/// This function returns the total cost of the found flow.
|
kpeter@806
|
648 |
/// Its complexity is O(e).
|
kpeter@806
|
649 |
///
|
kpeter@806
|
650 |
/// \note The return type of the function can be specified as a
|
kpeter@806
|
651 |
/// template parameter. For example,
|
kpeter@806
|
652 |
/// \code
|
kpeter@806
|
653 |
/// cs.totalCost<double>();
|
kpeter@806
|
654 |
/// \endcode
|
kpeter@806
|
655 |
/// It is useful if the total cost cannot be stored in the \c Cost
|
kpeter@806
|
656 |
/// type of the algorithm, which is the default return type of the
|
kpeter@806
|
657 |
/// function.
|
kpeter@805
|
658 |
///
|
kpeter@805
|
659 |
/// \pre \ref run() must be called before using this function.
|
kpeter@806
|
660 |
template <typename Number>
|
kpeter@806
|
661 |
Number totalCost() const {
|
kpeter@806
|
662 |
Number c = 0;
|
kpeter@806
|
663 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
kpeter@806
|
664 |
int i = _arc_idb[a];
|
kpeter@806
|
665 |
c += static_cast<Number>(_res_cap[i]) *
|
kpeter@806
|
666 |
(-static_cast<Number>(_cost[i]));
|
kpeter@806
|
667 |
}
|
kpeter@806
|
668 |
return c;
|
kpeter@805
|
669 |
}
|
kpeter@805
|
670 |
|
kpeter@806
|
671 |
#ifndef DOXYGEN
|
kpeter@806
|
672 |
Cost totalCost() const {
|
kpeter@806
|
673 |
return totalCost<Cost>();
|
kpeter@805
|
674 |
}
|
kpeter@806
|
675 |
#endif
|
kpeter@805
|
676 |
|
kpeter@805
|
677 |
/// \brief Return the flow on the given arc.
|
kpeter@805
|
678 |
///
|
kpeter@806
|
679 |
/// This function returns the flow on the given arc.
|
kpeter@805
|
680 |
///
|
kpeter@805
|
681 |
/// \pre \ref run() must be called before using this function.
|
kpeter@806
|
682 |
Value flow(const Arc& a) const {
|
kpeter@806
|
683 |
return _res_cap[_arc_idb[a]];
|
kpeter@805
|
684 |
}
|
kpeter@805
|
685 |
|
kpeter@1003
|
686 |
/// \brief Copy the flow values (the primal solution) into the
|
kpeter@1003
|
687 |
/// given map.
|
kpeter@805
|
688 |
///
|
kpeter@806
|
689 |
/// This function copies the flow value on each arc into the given
|
kpeter@806
|
690 |
/// map. The \c Value type of the algorithm must be convertible to
|
kpeter@806
|
691 |
/// the \c Value type of the map.
|
kpeter@805
|
692 |
///
|
kpeter@805
|
693 |
/// \pre \ref run() must be called before using this function.
|
kpeter@806
|
694 |
template <typename FlowMap>
|
kpeter@806
|
695 |
void flowMap(FlowMap &map) const {
|
kpeter@806
|
696 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
kpeter@806
|
697 |
map.set(a, _res_cap[_arc_idb[a]]);
|
kpeter@806
|
698 |
}
|
kpeter@805
|
699 |
}
|
kpeter@805
|
700 |
|
kpeter@806
|
701 |
/// \brief Return the potential (dual value) of the given node.
|
kpeter@805
|
702 |
///
|
kpeter@806
|
703 |
/// This function returns the potential (dual value) of the
|
kpeter@806
|
704 |
/// given node.
|
kpeter@805
|
705 |
///
|
kpeter@805
|
706 |
/// \pre \ref run() must be called before using this function.
|
kpeter@806
|
707 |
Cost potential(const Node& n) const {
|
kpeter@806
|
708 |
return _pi[_node_id[n]];
|
kpeter@806
|
709 |
}
|
kpeter@806
|
710 |
|
kpeter@1003
|
711 |
/// \brief Copy the potential values (the dual solution) into the
|
kpeter@1003
|
712 |
/// given map.
|
kpeter@806
|
713 |
///
|
kpeter@806
|
714 |
/// This function copies the potential (dual value) of each node
|
kpeter@806
|
715 |
/// into the given map.
|
kpeter@806
|
716 |
/// The \c Cost type of the algorithm must be convertible to the
|
kpeter@806
|
717 |
/// \c Value type of the map.
|
kpeter@806
|
718 |
///
|
kpeter@806
|
719 |
/// \pre \ref run() must be called before using this function.
|
kpeter@806
|
720 |
template <typename PotentialMap>
|
kpeter@806
|
721 |
void potentialMap(PotentialMap &map) const {
|
kpeter@806
|
722 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@806
|
723 |
map.set(n, _pi[_node_id[n]]);
|
kpeter@806
|
724 |
}
|
kpeter@805
|
725 |
}
|
kpeter@805
|
726 |
|
kpeter@805
|
727 |
/// @}
|
kpeter@805
|
728 |
|
kpeter@805
|
729 |
private:
|
kpeter@805
|
730 |
|
kpeter@806
|
731 |
// Initialize the algorithm
|
kpeter@810
|
732 |
ProblemType init() {
|
kpeter@821
|
733 |
if (_node_num <= 1) return INFEASIBLE;
|
kpeter@805
|
734 |
|
kpeter@806
|
735 |
// Check the sum of supply values
|
kpeter@806
|
736 |
_sum_supply = 0;
|
kpeter@806
|
737 |
for (int i = 0; i != _root; ++i) {
|
kpeter@806
|
738 |
_sum_supply += _supply[i];
|
kpeter@805
|
739 |
}
|
kpeter@806
|
740 |
if (_sum_supply > 0) return INFEASIBLE;
|
alpar@877
|
741 |
|
kpeter@811
|
742 |
// Initialize vectors
|
kpeter@806
|
743 |
for (int i = 0; i != _root; ++i) {
|
kpeter@806
|
744 |
_pi[i] = 0;
|
kpeter@806
|
745 |
_excess[i] = _supply[i];
|
kpeter@805
|
746 |
}
|
kpeter@805
|
747 |
|
kpeter@806
|
748 |
// Remove non-zero lower bounds
|
kpeter@811
|
749 |
const Value MAX = std::numeric_limits<Value>::max();
|
kpeter@811
|
750 |
int last_out;
|
kpeter@806
|
751 |
if (_have_lower) {
|
kpeter@806
|
752 |
for (int i = 0; i != _root; ++i) {
|
kpeter@811
|
753 |
last_out = _first_out[i+1];
|
kpeter@811
|
754 |
for (int j = _first_out[i]; j != last_out; ++j) {
|
kpeter@806
|
755 |
if (_forward[j]) {
|
kpeter@806
|
756 |
Value c = _lower[j];
|
kpeter@806
|
757 |
if (c >= 0) {
|
kpeter@811
|
758 |
_res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
|
kpeter@806
|
759 |
} else {
|
kpeter@811
|
760 |
_res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
|
kpeter@806
|
761 |
}
|
kpeter@806
|
762 |
_excess[i] -= c;
|
kpeter@806
|
763 |
_excess[_target[j]] += c;
|
kpeter@806
|
764 |
} else {
|
kpeter@806
|
765 |
_res_cap[j] = 0;
|
kpeter@806
|
766 |
}
|
kpeter@806
|
767 |
}
|
kpeter@806
|
768 |
}
|
kpeter@806
|
769 |
} else {
|
kpeter@806
|
770 |
for (int j = 0; j != _res_arc_num; ++j) {
|
kpeter@806
|
771 |
_res_cap[j] = _forward[j] ? _upper[j] : 0;
|
kpeter@806
|
772 |
}
|
kpeter@806
|
773 |
}
|
kpeter@805
|
774 |
|
kpeter@806
|
775 |
// Handle negative costs
|
kpeter@811
|
776 |
for (int i = 0; i != _root; ++i) {
|
kpeter@811
|
777 |
last_out = _first_out[i+1] - 1;
|
kpeter@811
|
778 |
for (int j = _first_out[i]; j != last_out; ++j) {
|
kpeter@811
|
779 |
Value rc = _res_cap[j];
|
kpeter@811
|
780 |
if (_cost[j] < 0 && rc > 0) {
|
kpeter@811
|
781 |
if (rc >= MAX) return UNBOUNDED;
|
kpeter@811
|
782 |
_excess[i] -= rc;
|
kpeter@811
|
783 |
_excess[_target[j]] += rc;
|
kpeter@811
|
784 |
_res_cap[j] = 0;
|
kpeter@811
|
785 |
_res_cap[_reverse[j]] += rc;
|
kpeter@806
|
786 |
}
|
kpeter@806
|
787 |
}
|
kpeter@806
|
788 |
}
|
alpar@877
|
789 |
|
kpeter@806
|
790 |
// Handle GEQ supply type
|
kpeter@806
|
791 |
if (_sum_supply < 0) {
|
kpeter@806
|
792 |
_pi[_root] = 0;
|
kpeter@806
|
793 |
_excess[_root] = -_sum_supply;
|
kpeter@806
|
794 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
|
kpeter@811
|
795 |
int ra = _reverse[a];
|
kpeter@811
|
796 |
_res_cap[a] = -_sum_supply + 1;
|
kpeter@811
|
797 |
_res_cap[ra] = 0;
|
kpeter@806
|
798 |
_cost[a] = 0;
|
kpeter@811
|
799 |
_cost[ra] = 0;
|
kpeter@806
|
800 |
}
|
kpeter@806
|
801 |
} else {
|
kpeter@806
|
802 |
_pi[_root] = 0;
|
kpeter@806
|
803 |
_excess[_root] = 0;
|
kpeter@806
|
804 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
|
kpeter@811
|
805 |
int ra = _reverse[a];
|
kpeter@806
|
806 |
_res_cap[a] = 1;
|
kpeter@811
|
807 |
_res_cap[ra] = 0;
|
kpeter@806
|
808 |
_cost[a] = 0;
|
kpeter@811
|
809 |
_cost[ra] = 0;
|
kpeter@806
|
810 |
}
|
kpeter@806
|
811 |
}
|
kpeter@806
|
812 |
|
kpeter@806
|
813 |
// Initialize delta value
|
kpeter@810
|
814 |
if (_factor > 1) {
|
kpeter@805
|
815 |
// With scaling
|
kpeter@839
|
816 |
Value max_sup = 0, max_dem = 0, max_cap = 0;
|
kpeter@839
|
817 |
for (int i = 0; i != _root; ++i) {
|
kpeter@811
|
818 |
Value ex = _excess[i];
|
kpeter@811
|
819 |
if ( ex > max_sup) max_sup = ex;
|
kpeter@811
|
820 |
if (-ex > max_dem) max_dem = -ex;
|
kpeter@839
|
821 |
int last_out = _first_out[i+1] - 1;
|
kpeter@839
|
822 |
for (int j = _first_out[i]; j != last_out; ++j) {
|
kpeter@839
|
823 |
if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
|
kpeter@839
|
824 |
}
|
kpeter@805
|
825 |
}
|
kpeter@805
|
826 |
max_sup = std::min(std::min(max_sup, max_dem), max_cap);
|
kpeter@810
|
827 |
for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
|
kpeter@805
|
828 |
} else {
|
kpeter@805
|
829 |
// Without scaling
|
kpeter@805
|
830 |
_delta = 1;
|
kpeter@805
|
831 |
}
|
kpeter@805
|
832 |
|
kpeter@806
|
833 |
return OPTIMAL;
|
kpeter@805
|
834 |
}
|
kpeter@805
|
835 |
|
kpeter@806
|
836 |
ProblemType start() {
|
kpeter@806
|
837 |
// Execute the algorithm
|
kpeter@806
|
838 |
ProblemType pt;
|
kpeter@805
|
839 |
if (_delta > 1)
|
kpeter@806
|
840 |
pt = startWithScaling();
|
kpeter@805
|
841 |
else
|
kpeter@806
|
842 |
pt = startWithoutScaling();
|
kpeter@806
|
843 |
|
kpeter@806
|
844 |
// Handle non-zero lower bounds
|
kpeter@806
|
845 |
if (_have_lower) {
|
kpeter@811
|
846 |
int limit = _first_out[_root];
|
kpeter@811
|
847 |
for (int j = 0; j != limit; ++j) {
|
kpeter@806
|
848 |
if (!_forward[j]) _res_cap[j] += _lower[j];
|
kpeter@806
|
849 |
}
|
kpeter@806
|
850 |
}
|
kpeter@806
|
851 |
|
kpeter@806
|
852 |
// Shift potentials if necessary
|
kpeter@806
|
853 |
Cost pr = _pi[_root];
|
kpeter@806
|
854 |
if (_sum_supply < 0 || pr > 0) {
|
kpeter@806
|
855 |
for (int i = 0; i != _node_num; ++i) {
|
kpeter@806
|
856 |
_pi[i] -= pr;
|
alpar@877
|
857 |
}
|
kpeter@806
|
858 |
}
|
alpar@877
|
859 |
|
kpeter@806
|
860 |
return pt;
|
kpeter@805
|
861 |
}
|
kpeter@805
|
862 |
|
kpeter@806
|
863 |
// Execute the capacity scaling algorithm
|
kpeter@806
|
864 |
ProblemType startWithScaling() {
|
kpeter@807
|
865 |
// Perform capacity scaling phases
|
kpeter@806
|
866 |
int s, t;
|
kpeter@806
|
867 |
ResidualDijkstra _dijkstra(*this);
|
kpeter@805
|
868 |
while (true) {
|
kpeter@806
|
869 |
// Saturate all arcs not satisfying the optimality condition
|
kpeter@811
|
870 |
int last_out;
|
kpeter@806
|
871 |
for (int u = 0; u != _node_num; ++u) {
|
kpeter@811
|
872 |
last_out = _sum_supply < 0 ?
|
kpeter@811
|
873 |
_first_out[u+1] : _first_out[u+1] - 1;
|
kpeter@811
|
874 |
for (int a = _first_out[u]; a != last_out; ++a) {
|
kpeter@806
|
875 |
int v = _target[a];
|
kpeter@806
|
876 |
Cost c = _cost[a] + _pi[u] - _pi[v];
|
kpeter@806
|
877 |
Value rc = _res_cap[a];
|
kpeter@806
|
878 |
if (c < 0 && rc >= _delta) {
|
kpeter@806
|
879 |
_excess[u] -= rc;
|
kpeter@806
|
880 |
_excess[v] += rc;
|
kpeter@806
|
881 |
_res_cap[a] = 0;
|
kpeter@806
|
882 |
_res_cap[_reverse[a]] += rc;
|
kpeter@806
|
883 |
}
|
kpeter@805
|
884 |
}
|
kpeter@805
|
885 |
}
|
kpeter@805
|
886 |
|
kpeter@806
|
887 |
// Find excess nodes and deficit nodes
|
kpeter@805
|
888 |
_excess_nodes.clear();
|
kpeter@805
|
889 |
_deficit_nodes.clear();
|
kpeter@806
|
890 |
for (int u = 0; u != _node_num; ++u) {
|
kpeter@811
|
891 |
Value ex = _excess[u];
|
kpeter@811
|
892 |
if (ex >= _delta) _excess_nodes.push_back(u);
|
kpeter@811
|
893 |
if (ex <= -_delta) _deficit_nodes.push_back(u);
|
kpeter@805
|
894 |
}
|
kpeter@805
|
895 |
int next_node = 0, next_def_node = 0;
|
kpeter@805
|
896 |
|
kpeter@806
|
897 |
// Find augmenting shortest paths
|
kpeter@805
|
898 |
while (next_node < int(_excess_nodes.size())) {
|
kpeter@806
|
899 |
// Check deficit nodes
|
kpeter@805
|
900 |
if (_delta > 1) {
|
kpeter@805
|
901 |
bool delta_deficit = false;
|
kpeter@805
|
902 |
for ( ; next_def_node < int(_deficit_nodes.size());
|
kpeter@805
|
903 |
++next_def_node ) {
|
kpeter@805
|
904 |
if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
|
kpeter@805
|
905 |
delta_deficit = true;
|
kpeter@805
|
906 |
break;
|
kpeter@805
|
907 |
}
|
kpeter@805
|
908 |
}
|
kpeter@805
|
909 |
if (!delta_deficit) break;
|
kpeter@805
|
910 |
}
|
kpeter@805
|
911 |
|
kpeter@806
|
912 |
// Run Dijkstra in the residual network
|
kpeter@805
|
913 |
s = _excess_nodes[next_node];
|
kpeter@806
|
914 |
if ((t = _dijkstra.run(s, _delta)) == -1) {
|
kpeter@805
|
915 |
if (_delta > 1) {
|
kpeter@805
|
916 |
++next_node;
|
kpeter@805
|
917 |
continue;
|
kpeter@805
|
918 |
}
|
kpeter@806
|
919 |
return INFEASIBLE;
|
kpeter@805
|
920 |
}
|
kpeter@805
|
921 |
|
kpeter@806
|
922 |
// Augment along a shortest path from s to t
|
kpeter@806
|
923 |
Value d = std::min(_excess[s], -_excess[t]);
|
kpeter@806
|
924 |
int u = t;
|
kpeter@806
|
925 |
int a;
|
kpeter@805
|
926 |
if (d > _delta) {
|
kpeter@806
|
927 |
while ((a = _pred[u]) != -1) {
|
kpeter@806
|
928 |
if (_res_cap[a] < d) d = _res_cap[a];
|
kpeter@806
|
929 |
u = _source[a];
|
kpeter@805
|
930 |
}
|
kpeter@805
|
931 |
}
|
kpeter@805
|
932 |
u = t;
|
kpeter@806
|
933 |
while ((a = _pred[u]) != -1) {
|
kpeter@806
|
934 |
_res_cap[a] -= d;
|
kpeter@806
|
935 |
_res_cap[_reverse[a]] += d;
|
kpeter@806
|
936 |
u = _source[a];
|
kpeter@805
|
937 |
}
|
kpeter@805
|
938 |
_excess[s] -= d;
|
kpeter@805
|
939 |
_excess[t] += d;
|
kpeter@805
|
940 |
|
kpeter@805
|
941 |
if (_excess[s] < _delta) ++next_node;
|
kpeter@805
|
942 |
}
|
kpeter@805
|
943 |
|
kpeter@805
|
944 |
if (_delta == 1) break;
|
kpeter@810
|
945 |
_delta = _delta <= _factor ? 1 : _delta / _factor;
|
kpeter@805
|
946 |
}
|
kpeter@805
|
947 |
|
kpeter@806
|
948 |
return OPTIMAL;
|
kpeter@805
|
949 |
}
|
kpeter@805
|
950 |
|
kpeter@806
|
951 |
// Execute the successive shortest path algorithm
|
kpeter@806
|
952 |
ProblemType startWithoutScaling() {
|
kpeter@806
|
953 |
// Find excess nodes
|
kpeter@806
|
954 |
_excess_nodes.clear();
|
kpeter@806
|
955 |
for (int i = 0; i != _node_num; ++i) {
|
kpeter@806
|
956 |
if (_excess[i] > 0) _excess_nodes.push_back(i);
|
kpeter@806
|
957 |
}
|
kpeter@806
|
958 |
if (_excess_nodes.size() == 0) return OPTIMAL;
|
kpeter@805
|
959 |
int next_node = 0;
|
kpeter@805
|
960 |
|
kpeter@806
|
961 |
// Find shortest paths
|
kpeter@806
|
962 |
int s, t;
|
kpeter@806
|
963 |
ResidualDijkstra _dijkstra(*this);
|
kpeter@805
|
964 |
while ( _excess[_excess_nodes[next_node]] > 0 ||
|
kpeter@805
|
965 |
++next_node < int(_excess_nodes.size()) )
|
kpeter@805
|
966 |
{
|
kpeter@806
|
967 |
// Run Dijkstra in the residual network
|
kpeter@805
|
968 |
s = _excess_nodes[next_node];
|
kpeter@806
|
969 |
if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
|
kpeter@805
|
970 |
|
kpeter@806
|
971 |
// Augment along a shortest path from s to t
|
kpeter@806
|
972 |
Value d = std::min(_excess[s], -_excess[t]);
|
kpeter@806
|
973 |
int u = t;
|
kpeter@806
|
974 |
int a;
|
kpeter@805
|
975 |
if (d > 1) {
|
kpeter@806
|
976 |
while ((a = _pred[u]) != -1) {
|
kpeter@806
|
977 |
if (_res_cap[a] < d) d = _res_cap[a];
|
kpeter@806
|
978 |
u = _source[a];
|
kpeter@805
|
979 |
}
|
kpeter@805
|
980 |
}
|
kpeter@805
|
981 |
u = t;
|
kpeter@806
|
982 |
while ((a = _pred[u]) != -1) {
|
kpeter@806
|
983 |
_res_cap[a] -= d;
|
kpeter@806
|
984 |
_res_cap[_reverse[a]] += d;
|
kpeter@806
|
985 |
u = _source[a];
|
kpeter@805
|
986 |
}
|
kpeter@805
|
987 |
_excess[s] -= d;
|
kpeter@805
|
988 |
_excess[t] += d;
|
kpeter@805
|
989 |
}
|
kpeter@805
|
990 |
|
kpeter@806
|
991 |
return OPTIMAL;
|
kpeter@805
|
992 |
}
|
kpeter@805
|
993 |
|
kpeter@805
|
994 |
}; //class CapacityScaling
|
kpeter@805
|
995 |
|
kpeter@805
|
996 |
///@}
|
kpeter@805
|
997 |
|
kpeter@805
|
998 |
} //namespace lemon
|
kpeter@805
|
999 |
|
kpeter@805
|
1000 |
#endif //LEMON_CAPACITY_SCALING_H
|