1 |
1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*-
|
2 |
2 |
*
|
3 |
3 |
* This file is a part of LEMON, a generic C++ optimization library.
|
4 |
4 |
*
|
5 |
5 |
* Copyright (C) 2003-2009
|
6 |
6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
|
7 |
7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES).
|
8 |
8 |
*
|
9 |
9 |
* Permission to use, modify and distribute this software is granted
|
10 |
10 |
* provided that this copyright notice appears in all copies. For
|
11 |
11 |
* precise terms see the accompanying LICENSE file.
|
12 |
12 |
*
|
13 |
13 |
* This software is provided "AS IS" with no warranty of any kind,
|
14 |
14 |
* express or implied, and with no claim as to its suitability for any
|
15 |
15 |
* purpose.
|
16 |
16 |
*
|
17 |
17 |
*/
|
18 |
18 |
|
19 |
19 |
#ifndef LEMON_NETWORK_SIMPLEX_H
|
20 |
20 |
#define LEMON_NETWORK_SIMPLEX_H
|
21 |
21 |
|
22 |
|
/// \ingroup min_cost_flow
|
|
22 |
/// \ingroup min_cost_flow_algs
|
23 |
23 |
///
|
24 |
24 |
/// \file
|
25 |
25 |
/// \brief Network Simplex algorithm for finding a minimum cost flow.
|
26 |
26 |
|
27 |
27 |
#include <vector>
|
28 |
28 |
#include <limits>
|
29 |
29 |
#include <algorithm>
|
30 |
30 |
|
31 |
31 |
#include <lemon/core.h>
|
32 |
32 |
#include <lemon/math.h>
|
33 |
33 |
|
34 |
34 |
namespace lemon {
|
35 |
35 |
|
36 |
|
/// \addtogroup min_cost_flow
|
|
36 |
/// \addtogroup min_cost_flow_algs
|
37 |
37 |
/// @{
|
38 |
38 |
|
39 |
39 |
/// \brief Implementation of the primal Network Simplex algorithm
|
40 |
40 |
/// for finding a \ref min_cost_flow "minimum cost flow".
|
41 |
41 |
///
|
42 |
42 |
/// \ref NetworkSimplex implements the primal Network Simplex algorithm
|
43 |
43 |
/// for finding a \ref min_cost_flow "minimum cost flow".
|
44 |
44 |
/// This algorithm is a specialized version of the linear programming
|
45 |
45 |
/// simplex method directly for the minimum cost flow problem.
|
46 |
46 |
/// It is one of the most efficient solution methods.
|
47 |
47 |
///
|
48 |
48 |
/// In general this class is the fastest implementation available
|
49 |
49 |
/// in LEMON for the minimum cost flow problem.
|
50 |
50 |
/// Moreover it supports both directions of the supply/demand inequality
|
51 |
51 |
/// constraints. For more information see \ref SupplyType.
|
52 |
52 |
///
|
53 |
53 |
/// Most of the parameters of the problem (except for the digraph)
|
54 |
54 |
/// can be given using separate functions, and the algorithm can be
|
55 |
55 |
/// executed using the \ref run() function. If some parameters are not
|
56 |
56 |
/// specified, then default values will be used.
|
57 |
57 |
///
|
58 |
58 |
/// \tparam GR The digraph type the algorithm runs on.
|
59 |
59 |
/// \tparam V The value type used for flow amounts, capacity bounds
|
60 |
60 |
/// and supply values in the algorithm. By default it is \c int.
|
61 |
61 |
/// \tparam C The value type used for costs and potentials in the
|
62 |
62 |
/// algorithm. By default it is the same as \c V.
|
63 |
63 |
///
|
64 |
64 |
/// \warning Both value types must be signed and all input data must
|
65 |
65 |
/// be integer.
|
66 |
66 |
///
|
67 |
67 |
/// \note %NetworkSimplex provides five different pivot rule
|
68 |
68 |
/// implementations, from which the most efficient one is used
|
69 |
69 |
/// by default. For more information see \ref PivotRule.
|
70 |
70 |
template <typename GR, typename V = int, typename C = V>
|
71 |
71 |
class NetworkSimplex
|
72 |
72 |
{
|
73 |
73 |
public:
|
74 |
74 |
|
75 |
75 |
/// The type of the flow amounts, capacity bounds and supply values
|
76 |
76 |
typedef V Value;
|
77 |
77 |
/// The type of the arc costs
|
78 |
78 |
typedef C Cost;
|
79 |
79 |
|
80 |
80 |
public:
|
81 |
81 |
|
82 |
82 |
/// \brief Problem type constants for the \c run() function.
|
83 |
83 |
///
|
84 |
84 |
/// Enum type containing the problem type constants that can be
|
85 |
85 |
/// returned by the \ref run() function of the algorithm.
|
86 |
86 |
enum ProblemType {
|
87 |
87 |
/// The problem has no feasible solution (flow).
|
88 |
88 |
INFEASIBLE,
|
89 |
89 |
/// The problem has optimal solution (i.e. it is feasible and
|
90 |
90 |
/// bounded), and the algorithm has found optimal flow and node
|
91 |
91 |
/// potentials (primal and dual solutions).
|
92 |
92 |
OPTIMAL,
|
93 |
93 |
/// The objective function of the problem is unbounded, i.e.
|
94 |
94 |
/// there is a directed cycle having negative total cost and
|
95 |
95 |
/// infinite upper bound.
|
96 |
96 |
UNBOUNDED
|
97 |
97 |
};
|
98 |
98 |
|
99 |
99 |
/// \brief Constants for selecting the type of the supply constraints.
|
100 |
100 |
///
|
101 |
101 |
/// Enum type containing constants for selecting the supply type,
|
102 |
102 |
/// i.e. the direction of the inequalities in the supply/demand
|
103 |
103 |
/// constraints of the \ref min_cost_flow "minimum cost flow problem".
|
104 |
104 |
///
|
105 |
|
/// The default supply type is \c GEQ, since this form is supported
|
106 |
|
/// by other minimum cost flow algorithms and the \ref Circulation
|
107 |
|
/// algorithm, as well.
|
108 |
|
/// The \c LEQ problem type can be selected using the \ref supplyType()
|
109 |
|
/// function.
|
110 |
|
///
|
111 |
|
/// Note that the equality form is a special case of both supply types.
|
|
105 |
/// The default supply type is \c GEQ, the \c LEQ type can be
|
|
106 |
/// selected using \ref supplyType().
|
|
107 |
/// The equality form is a special case of both supply types.
|
112 |
108 |
enum SupplyType {
|
113 |
|
|
114 |
109 |
/// This option means that there are <em>"greater or equal"</em>
|
115 |
|
/// supply/demand constraints in the definition, i.e. the exact
|
116 |
|
/// formulation of the problem is the following.
|
117 |
|
/**
|
118 |
|
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
|
119 |
|
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
|
120 |
|
sup(u) \quad \forall u\in V \f]
|
121 |
|
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
|
122 |
|
*/
|
123 |
|
/// It means that the total demand must be greater or equal to the
|
124 |
|
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
|
125 |
|
/// negative) and all the supplies have to be carried out from
|
126 |
|
/// the supply nodes, but there could be demands that are not
|
127 |
|
/// satisfied.
|
|
110 |
/// supply/demand constraints in the definition of the problem.
|
128 |
111 |
GEQ,
|
129 |
|
/// It is just an alias for the \c GEQ option.
|
130 |
|
CARRY_SUPPLIES = GEQ,
|
131 |
|
|
132 |
112 |
/// This option means that there are <em>"less or equal"</em>
|
133 |
|
/// supply/demand constraints in the definition, i.e. the exact
|
134 |
|
/// formulation of the problem is the following.
|
135 |
|
/**
|
136 |
|
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
|
137 |
|
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
|
138 |
|
sup(u) \quad \forall u\in V \f]
|
139 |
|
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
|
140 |
|
*/
|
141 |
|
/// It means that the total demand must be less or equal to the
|
142 |
|
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
|
143 |
|
/// positive) and all the demands have to be satisfied, but there
|
144 |
|
/// could be supplies that are not carried out from the supply
|
145 |
|
/// nodes.
|
146 |
|
LEQ,
|
147 |
|
/// It is just an alias for the \c LEQ option.
|
148 |
|
SATISFY_DEMANDS = LEQ
|
|
113 |
/// supply/demand constraints in the definition of the problem.
|
|
114 |
LEQ
|
149 |
115 |
};
|
150 |
116 |
|
151 |
117 |
/// \brief Constants for selecting the pivot rule.
|
152 |
118 |
///
|
153 |
119 |
/// Enum type containing constants for selecting the pivot rule for
|
154 |
120 |
/// the \ref run() function.
|
155 |
121 |
///
|
156 |
122 |
/// \ref NetworkSimplex provides five different pivot rule
|
157 |
123 |
/// implementations that significantly affect the running time
|
158 |
124 |
/// of the algorithm.
|
159 |
125 |
/// By default \ref BLOCK_SEARCH "Block Search" is used, which
|
160 |
126 |
/// proved to be the most efficient and the most robust on various
|
161 |
127 |
/// test inputs according to our benchmark tests.
|
162 |
128 |
/// However another pivot rule can be selected using the \ref run()
|
163 |
129 |
/// function with the proper parameter.
|
164 |
130 |
enum PivotRule {
|
165 |
131 |
|
166 |
132 |
/// The First Eligible pivot rule.
|
167 |
133 |
/// The next eligible arc is selected in a wraparound fashion
|
168 |
134 |
/// in every iteration.
|
169 |
135 |
FIRST_ELIGIBLE,
|
170 |
136 |
|
171 |
137 |
/// The Best Eligible pivot rule.
|
172 |
138 |
/// The best eligible arc is selected in every iteration.
|
173 |
139 |
BEST_ELIGIBLE,
|
174 |
140 |
|
175 |
141 |
/// The Block Search pivot rule.
|
176 |
142 |
/// A specified number of arcs are examined in every iteration
|
177 |
143 |
/// in a wraparound fashion and the best eligible arc is selected
|
178 |
144 |
/// from this block.
|
179 |
145 |
BLOCK_SEARCH,
|
180 |
146 |
|
181 |
147 |
/// The Candidate List pivot rule.
|
182 |
148 |
/// In a major iteration a candidate list is built from eligible arcs
|
183 |
149 |
/// in a wraparound fashion and in the following minor iterations
|
184 |
150 |
/// the best eligible arc is selected from this list.
|
185 |
151 |
CANDIDATE_LIST,
|
186 |
152 |
|
187 |
153 |
/// The Altering Candidate List pivot rule.
|
188 |
154 |
/// It is a modified version of the Candidate List method.
|
189 |
155 |
/// It keeps only the several best eligible arcs from the former
|
190 |
156 |
/// candidate list and extends this list in every iteration.
|
191 |
157 |
ALTERING_LIST
|
192 |
158 |
};
|
193 |
159 |
|
194 |
160 |
private:
|
195 |
161 |
|
196 |
162 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR);
|
197 |
163 |
|
198 |
164 |
typedef std::vector<Arc> ArcVector;
|
199 |
165 |
typedef std::vector<Node> NodeVector;
|
200 |
166 |
typedef std::vector<int> IntVector;
|
201 |
167 |
typedef std::vector<bool> BoolVector;
|
202 |
168 |
typedef std::vector<Value> ValueVector;
|
203 |
169 |
typedef std::vector<Cost> CostVector;
|
204 |
170 |
|
205 |
171 |
// State constants for arcs
|
206 |
172 |
enum ArcStateEnum {
|
207 |
173 |
STATE_UPPER = -1,
|
208 |
174 |
STATE_TREE = 0,
|
209 |
175 |
STATE_LOWER = 1
|
210 |
176 |
};
|
211 |
177 |
|
212 |
178 |
private:
|
213 |
179 |
|
214 |
180 |
// Data related to the underlying digraph
|
215 |
181 |
const GR &_graph;
|
216 |
182 |
int _node_num;
|
217 |
183 |
int _arc_num;
|
|
184 |
int _all_arc_num;
|
|
185 |
int _search_arc_num;
|
218 |
186 |
|
219 |
187 |
// Parameters of the problem
|
220 |
188 |
bool _have_lower;
|
221 |
189 |
SupplyType _stype;
|
222 |
190 |
Value _sum_supply;
|
223 |
191 |
|
224 |
192 |
// Data structures for storing the digraph
|
225 |
193 |
IntNodeMap _node_id;
|
226 |
194 |
IntArcMap _arc_id;
|
227 |
195 |
IntVector _source;
|
228 |
196 |
IntVector _target;
|
229 |
197 |
|
230 |
198 |
// Node and arc data
|
231 |
199 |
ValueVector _lower;
|
232 |
200 |
ValueVector _upper;
|
233 |
201 |
ValueVector _cap;
|
234 |
202 |
CostVector _cost;
|
235 |
203 |
ValueVector _supply;
|
236 |
204 |
ValueVector _flow;
|
237 |
205 |
CostVector _pi;
|
238 |
206 |
|
239 |
207 |
// Data for storing the spanning tree structure
|
240 |
208 |
IntVector _parent;
|
241 |
209 |
IntVector _pred;
|
242 |
210 |
IntVector _thread;
|
243 |
211 |
IntVector _rev_thread;
|
244 |
212 |
IntVector _succ_num;
|
245 |
213 |
IntVector _last_succ;
|
246 |
214 |
IntVector _dirty_revs;
|
247 |
215 |
BoolVector _forward;
|
248 |
216 |
IntVector _state;
|
249 |
217 |
int _root;
|
250 |
218 |
|
251 |
219 |
// Temporary data used in the current pivot iteration
|
252 |
220 |
int in_arc, join, u_in, v_in, u_out, v_out;
|
253 |
221 |
int first, second, right, last;
|
254 |
222 |
int stem, par_stem, new_stem;
|
255 |
223 |
Value delta;
|
256 |
224 |
|
257 |
225 |
public:
|
258 |
226 |
|
259 |
227 |
/// \brief Constant for infinite upper bounds (capacities).
|
260 |
228 |
///
|
261 |
229 |
/// Constant for infinite upper bounds (capacities).
|
262 |
230 |
/// It is \c std::numeric_limits<Value>::infinity() if available,
|
263 |
231 |
/// \c std::numeric_limits<Value>::max() otherwise.
|
264 |
232 |
const Value INF;
|
265 |
233 |
|
266 |
234 |
private:
|
267 |
235 |
|
268 |
236 |
// Implementation of the First Eligible pivot rule
|
269 |
237 |
class FirstEligiblePivotRule
|
270 |
238 |
{
|
271 |
239 |
private:
|
272 |
240 |
|
273 |
241 |
// References to the NetworkSimplex class
|
274 |
242 |
const IntVector &_source;
|
275 |
243 |
const IntVector &_target;
|
276 |
244 |
const CostVector &_cost;
|
277 |
245 |
const IntVector &_state;
|
278 |
246 |
const CostVector &_pi;
|
279 |
247 |
int &_in_arc;
|
280 |
|
int _arc_num;
|
|
248 |
int _search_arc_num;
|
281 |
249 |
|
282 |
250 |
// Pivot rule data
|
283 |
251 |
int _next_arc;
|
284 |
252 |
|
285 |
253 |
public:
|
286 |
254 |
|
287 |
255 |
// Constructor
|
288 |
256 |
FirstEligiblePivotRule(NetworkSimplex &ns) :
|
289 |
257 |
_source(ns._source), _target(ns._target),
|
290 |
258 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
291 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
|
259 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
|
|
260 |
_next_arc(0)
|
292 |
261 |
{}
|
293 |
262 |
|
294 |
263 |
// Find next entering arc
|
295 |
264 |
bool findEnteringArc() {
|
296 |
265 |
Cost c;
|
297 |
|
for (int e = _next_arc; e < _arc_num; ++e) {
|
|
266 |
for (int e = _next_arc; e < _search_arc_num; ++e) {
|
298 |
267 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
299 |
268 |
if (c < 0) {
|
300 |
269 |
_in_arc = e;
|
301 |
270 |
_next_arc = e + 1;
|
302 |
271 |
return true;
|
303 |
272 |
}
|
304 |
273 |
}
|
305 |
274 |
for (int e = 0; e < _next_arc; ++e) {
|
306 |
275 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
307 |
276 |
if (c < 0) {
|
308 |
277 |
_in_arc = e;
|
309 |
278 |
_next_arc = e + 1;
|
310 |
279 |
return true;
|
311 |
280 |
}
|
312 |
281 |
}
|
313 |
282 |
return false;
|
314 |
283 |
}
|
315 |
284 |
|
316 |
285 |
}; //class FirstEligiblePivotRule
|
317 |
286 |
|
318 |
287 |
|
319 |
288 |
// Implementation of the Best Eligible pivot rule
|
320 |
289 |
class BestEligiblePivotRule
|
321 |
290 |
{
|
322 |
291 |
private:
|
323 |
292 |
|
324 |
293 |
// References to the NetworkSimplex class
|
325 |
294 |
const IntVector &_source;
|
326 |
295 |
const IntVector &_target;
|
327 |
296 |
const CostVector &_cost;
|
328 |
297 |
const IntVector &_state;
|
329 |
298 |
const CostVector &_pi;
|
330 |
299 |
int &_in_arc;
|
331 |
|
int _arc_num;
|
|
300 |
int _search_arc_num;
|
332 |
301 |
|
333 |
302 |
public:
|
334 |
303 |
|
335 |
304 |
// Constructor
|
336 |
305 |
BestEligiblePivotRule(NetworkSimplex &ns) :
|
337 |
306 |
_source(ns._source), _target(ns._target),
|
338 |
307 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
339 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num)
|
|
308 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
|
340 |
309 |
{}
|
341 |
310 |
|
342 |
311 |
// Find next entering arc
|
343 |
312 |
bool findEnteringArc() {
|
344 |
313 |
Cost c, min = 0;
|
345 |
|
for (int e = 0; e < _arc_num; ++e) {
|
|
314 |
for (int e = 0; e < _search_arc_num; ++e) {
|
346 |
315 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
347 |
316 |
if (c < min) {
|
348 |
317 |
min = c;
|
349 |
318 |
_in_arc = e;
|
350 |
319 |
}
|
351 |
320 |
}
|
352 |
321 |
return min < 0;
|
353 |
322 |
}
|
354 |
323 |
|
355 |
324 |
}; //class BestEligiblePivotRule
|
356 |
325 |
|
357 |
326 |
|
358 |
327 |
// Implementation of the Block Search pivot rule
|
359 |
328 |
class BlockSearchPivotRule
|
360 |
329 |
{
|
361 |
330 |
private:
|
362 |
331 |
|
363 |
332 |
// References to the NetworkSimplex class
|
364 |
333 |
const IntVector &_source;
|
365 |
334 |
const IntVector &_target;
|
366 |
335 |
const CostVector &_cost;
|
367 |
336 |
const IntVector &_state;
|
368 |
337 |
const CostVector &_pi;
|
369 |
338 |
int &_in_arc;
|
370 |
|
int _arc_num;
|
|
339 |
int _search_arc_num;
|
371 |
340 |
|
372 |
341 |
// Pivot rule data
|
373 |
342 |
int _block_size;
|
374 |
343 |
int _next_arc;
|
375 |
344 |
|
376 |
345 |
public:
|
377 |
346 |
|
378 |
347 |
// Constructor
|
379 |
348 |
BlockSearchPivotRule(NetworkSimplex &ns) :
|
380 |
349 |
_source(ns._source), _target(ns._target),
|
381 |
350 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
382 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
|
351 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
|
|
352 |
_next_arc(0)
|
383 |
353 |
{
|
384 |
354 |
// The main parameters of the pivot rule
|
385 |
|
const double BLOCK_SIZE_FACTOR = 2.0;
|
|
355 |
const double BLOCK_SIZE_FACTOR = 0.5;
|
386 |
356 |
const int MIN_BLOCK_SIZE = 10;
|
387 |
357 |
|
388 |
358 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR *
|
389 |
|
std::sqrt(double(_arc_num))),
|
|
359 |
std::sqrt(double(_search_arc_num))),
|
390 |
360 |
MIN_BLOCK_SIZE );
|
391 |
361 |
}
|
392 |
362 |
|
393 |
363 |
// Find next entering arc
|
394 |
364 |
bool findEnteringArc() {
|
395 |
365 |
Cost c, min = 0;
|
396 |
366 |
int cnt = _block_size;
|
397 |
367 |
int e, min_arc = _next_arc;
|
398 |
|
for (e = _next_arc; e < _arc_num; ++e) {
|
|
368 |
for (e = _next_arc; e < _search_arc_num; ++e) {
|
399 |
369 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
400 |
370 |
if (c < min) {
|
401 |
371 |
min = c;
|
402 |
372 |
min_arc = e;
|
403 |
373 |
}
|
404 |
374 |
if (--cnt == 0) {
|
405 |
375 |
if (min < 0) break;
|
406 |
376 |
cnt = _block_size;
|
407 |
377 |
}
|
408 |
378 |
}
|
409 |
379 |
if (min == 0 || cnt > 0) {
|
410 |
380 |
for (e = 0; e < _next_arc; ++e) {
|
411 |
381 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
412 |
382 |
if (c < min) {
|
413 |
383 |
min = c;
|
414 |
384 |
min_arc = e;
|
415 |
385 |
}
|
416 |
386 |
if (--cnt == 0) {
|
417 |
387 |
if (min < 0) break;
|
418 |
388 |
cnt = _block_size;
|
419 |
389 |
}
|
420 |
390 |
}
|
421 |
391 |
}
|
422 |
392 |
if (min >= 0) return false;
|
423 |
393 |
_in_arc = min_arc;
|
424 |
394 |
_next_arc = e;
|
425 |
395 |
return true;
|
426 |
396 |
}
|
427 |
397 |
|
428 |
398 |
}; //class BlockSearchPivotRule
|
429 |
399 |
|
430 |
400 |
|
431 |
401 |
// Implementation of the Candidate List pivot rule
|
432 |
402 |
class CandidateListPivotRule
|
433 |
403 |
{
|
434 |
404 |
private:
|
435 |
405 |
|
436 |
406 |
// References to the NetworkSimplex class
|
437 |
407 |
const IntVector &_source;
|
438 |
408 |
const IntVector &_target;
|
439 |
409 |
const CostVector &_cost;
|
440 |
410 |
const IntVector &_state;
|
441 |
411 |
const CostVector &_pi;
|
442 |
412 |
int &_in_arc;
|
443 |
|
int _arc_num;
|
|
413 |
int _search_arc_num;
|
444 |
414 |
|
445 |
415 |
// Pivot rule data
|
446 |
416 |
IntVector _candidates;
|
447 |
417 |
int _list_length, _minor_limit;
|
448 |
418 |
int _curr_length, _minor_count;
|
449 |
419 |
int _next_arc;
|
450 |
420 |
|
451 |
421 |
public:
|
452 |
422 |
|
453 |
423 |
/// Constructor
|
454 |
424 |
CandidateListPivotRule(NetworkSimplex &ns) :
|
455 |
425 |
_source(ns._source), _target(ns._target),
|
456 |
426 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
457 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
|
427 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
|
|
428 |
_next_arc(0)
|
458 |
429 |
{
|
459 |
430 |
// The main parameters of the pivot rule
|
460 |
431 |
const double LIST_LENGTH_FACTOR = 1.0;
|
461 |
432 |
const int MIN_LIST_LENGTH = 10;
|
462 |
433 |
const double MINOR_LIMIT_FACTOR = 0.1;
|
463 |
434 |
const int MIN_MINOR_LIMIT = 3;
|
464 |
435 |
|
465 |
436 |
_list_length = std::max( int(LIST_LENGTH_FACTOR *
|
466 |
|
std::sqrt(double(_arc_num))),
|
|
437 |
std::sqrt(double(_search_arc_num))),
|
467 |
438 |
MIN_LIST_LENGTH );
|
468 |
439 |
_minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
|
469 |
440 |
MIN_MINOR_LIMIT );
|
470 |
441 |
_curr_length = _minor_count = 0;
|
471 |
442 |
_candidates.resize(_list_length);
|
472 |
443 |
}
|
473 |
444 |
|
474 |
445 |
/// Find next entering arc
|
475 |
446 |
bool findEnteringArc() {
|
476 |
447 |
Cost min, c;
|
477 |
448 |
int e, min_arc = _next_arc;
|
478 |
449 |
if (_curr_length > 0 && _minor_count < _minor_limit) {
|
479 |
450 |
// Minor iteration: select the best eligible arc from the
|
480 |
451 |
// current candidate list
|
481 |
452 |
++_minor_count;
|
482 |
453 |
min = 0;
|
483 |
454 |
for (int i = 0; i < _curr_length; ++i) {
|
484 |
455 |
e = _candidates[i];
|
485 |
456 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
486 |
457 |
if (c < min) {
|
487 |
458 |
min = c;
|
488 |
459 |
min_arc = e;
|
489 |
460 |
}
|
490 |
461 |
if (c >= 0) {
|
491 |
462 |
_candidates[i--] = _candidates[--_curr_length];
|
492 |
463 |
}
|
493 |
464 |
}
|
494 |
465 |
if (min < 0) {
|
495 |
466 |
_in_arc = min_arc;
|
496 |
467 |
return true;
|
497 |
468 |
}
|
498 |
469 |
}
|
499 |
470 |
|
500 |
471 |
// Major iteration: build a new candidate list
|
501 |
472 |
min = 0;
|
502 |
473 |
_curr_length = 0;
|
503 |
|
for (e = _next_arc; e < _arc_num; ++e) {
|
|
474 |
for (e = _next_arc; e < _search_arc_num; ++e) {
|
504 |
475 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
505 |
476 |
if (c < 0) {
|
506 |
477 |
_candidates[_curr_length++] = e;
|
507 |
478 |
if (c < min) {
|
508 |
479 |
min = c;
|
509 |
480 |
min_arc = e;
|
510 |
481 |
}
|
511 |
482 |
if (_curr_length == _list_length) break;
|
512 |
483 |
}
|
513 |
484 |
}
|
514 |
485 |
if (_curr_length < _list_length) {
|
515 |
486 |
for (e = 0; e < _next_arc; ++e) {
|
516 |
487 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
517 |
488 |
if (c < 0) {
|
518 |
489 |
_candidates[_curr_length++] = e;
|
519 |
490 |
if (c < min) {
|
520 |
491 |
min = c;
|
521 |
492 |
min_arc = e;
|
522 |
493 |
}
|
523 |
494 |
if (_curr_length == _list_length) break;
|
524 |
495 |
}
|
525 |
496 |
}
|
526 |
497 |
}
|
527 |
498 |
if (_curr_length == 0) return false;
|
528 |
499 |
_minor_count = 1;
|
529 |
500 |
_in_arc = min_arc;
|
530 |
501 |
_next_arc = e;
|
531 |
502 |
return true;
|
532 |
503 |
}
|
533 |
504 |
|
534 |
505 |
}; //class CandidateListPivotRule
|
535 |
506 |
|
536 |
507 |
|
537 |
508 |
// Implementation of the Altering Candidate List pivot rule
|
538 |
509 |
class AlteringListPivotRule
|
539 |
510 |
{
|
540 |
511 |
private:
|
541 |
512 |
|
542 |
513 |
// References to the NetworkSimplex class
|
543 |
514 |
const IntVector &_source;
|
544 |
515 |
const IntVector &_target;
|
545 |
516 |
const CostVector &_cost;
|
546 |
517 |
const IntVector &_state;
|
547 |
518 |
const CostVector &_pi;
|
548 |
519 |
int &_in_arc;
|
549 |
|
int _arc_num;
|
|
520 |
int _search_arc_num;
|
550 |
521 |
|
551 |
522 |
// Pivot rule data
|
552 |
523 |
int _block_size, _head_length, _curr_length;
|
553 |
524 |
int _next_arc;
|
554 |
525 |
IntVector _candidates;
|
555 |
526 |
CostVector _cand_cost;
|
556 |
527 |
|
557 |
528 |
// Functor class to compare arcs during sort of the candidate list
|
558 |
529 |
class SortFunc
|
559 |
530 |
{
|
560 |
531 |
private:
|
561 |
532 |
const CostVector &_map;
|
562 |
533 |
public:
|
563 |
534 |
SortFunc(const CostVector &map) : _map(map) {}
|
564 |
535 |
bool operator()(int left, int right) {
|
565 |
536 |
return _map[left] > _map[right];
|
566 |
537 |
}
|
567 |
538 |
};
|
568 |
539 |
|
569 |
540 |
SortFunc _sort_func;
|
570 |
541 |
|
571 |
542 |
public:
|
572 |
543 |
|
573 |
544 |
// Constructor
|
574 |
545 |
AlteringListPivotRule(NetworkSimplex &ns) :
|
575 |
546 |
_source(ns._source), _target(ns._target),
|
576 |
547 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
577 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num),
|
578 |
|
_next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
|
|
548 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
|
|
549 |
_next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
|
579 |
550 |
{
|
580 |
551 |
// The main parameters of the pivot rule
|
581 |
552 |
const double BLOCK_SIZE_FACTOR = 1.5;
|
582 |
553 |
const int MIN_BLOCK_SIZE = 10;
|
583 |
554 |
const double HEAD_LENGTH_FACTOR = 0.1;
|
584 |
555 |
const int MIN_HEAD_LENGTH = 3;
|
585 |
556 |
|
586 |
557 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR *
|
587 |
|
std::sqrt(double(_arc_num))),
|
|
558 |
std::sqrt(double(_search_arc_num))),
|
588 |
559 |
MIN_BLOCK_SIZE );
|
589 |
560 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
|
590 |
561 |
MIN_HEAD_LENGTH );
|
591 |
562 |
_candidates.resize(_head_length + _block_size);
|
592 |
563 |
_curr_length = 0;
|
593 |
564 |
}
|
594 |
565 |
|
595 |
566 |
// Find next entering arc
|
596 |
567 |
bool findEnteringArc() {
|
597 |
568 |
// Check the current candidate list
|
598 |
569 |
int e;
|
599 |
570 |
for (int i = 0; i < _curr_length; ++i) {
|
600 |
571 |
e = _candidates[i];
|
601 |
572 |
_cand_cost[e] = _state[e] *
|
602 |
573 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
603 |
574 |
if (_cand_cost[e] >= 0) {
|
604 |
575 |
_candidates[i--] = _candidates[--_curr_length];
|
605 |
576 |
}
|
606 |
577 |
}
|
607 |
578 |
|
608 |
579 |
// Extend the list
|
609 |
580 |
int cnt = _block_size;
|
610 |
581 |
int last_arc = 0;
|
611 |
582 |
int limit = _head_length;
|
612 |
583 |
|
613 |
|
for (int e = _next_arc; e < _arc_num; ++e) {
|
|
584 |
for (int e = _next_arc; e < _search_arc_num; ++e) {
|
614 |
585 |
_cand_cost[e] = _state[e] *
|
615 |
586 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
616 |
587 |
if (_cand_cost[e] < 0) {
|
617 |
588 |
_candidates[_curr_length++] = e;
|
618 |
589 |
last_arc = e;
|
619 |
590 |
}
|
620 |
591 |
if (--cnt == 0) {
|
621 |
592 |
if (_curr_length > limit) break;
|
622 |
593 |
limit = 0;
|
623 |
594 |
cnt = _block_size;
|
624 |
595 |
}
|
625 |
596 |
}
|
626 |
597 |
if (_curr_length <= limit) {
|
627 |
598 |
for (int e = 0; e < _next_arc; ++e) {
|
628 |
599 |
_cand_cost[e] = _state[e] *
|
629 |
600 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
630 |
601 |
if (_cand_cost[e] < 0) {
|
631 |
602 |
_candidates[_curr_length++] = e;
|
632 |
603 |
last_arc = e;
|
633 |
604 |
}
|
634 |
605 |
if (--cnt == 0) {
|
635 |
606 |
if (_curr_length > limit) break;
|
636 |
607 |
limit = 0;
|
637 |
608 |
cnt = _block_size;
|
638 |
609 |
}
|
639 |
610 |
}
|
640 |
611 |
}
|
641 |
612 |
if (_curr_length == 0) return false;
|
642 |
613 |
_next_arc = last_arc + 1;
|
643 |
614 |
|
644 |
615 |
// Make heap of the candidate list (approximating a partial sort)
|
645 |
616 |
make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
|
646 |
617 |
_sort_func );
|
647 |
618 |
|
648 |
619 |
// Pop the first element of the heap
|
649 |
620 |
_in_arc = _candidates[0];
|
650 |
621 |
pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
|
651 |
622 |
_sort_func );
|
652 |
623 |
_curr_length = std::min(_head_length, _curr_length - 1);
|
653 |
624 |
return true;
|
654 |
625 |
}
|
655 |
626 |
|
656 |
627 |
}; //class AlteringListPivotRule
|
657 |
628 |
|
658 |
629 |
public:
|
659 |
630 |
|
660 |
631 |
/// \brief Constructor.
|
661 |
632 |
///
|
662 |
633 |
/// The constructor of the class.
|
663 |
634 |
///
|
664 |
635 |
/// \param graph The digraph the algorithm runs on.
|
665 |
636 |
NetworkSimplex(const GR& graph) :
|
666 |
637 |
_graph(graph), _node_id(graph), _arc_id(graph),
|
667 |
638 |
INF(std::numeric_limits<Value>::has_infinity ?
|
668 |
639 |
std::numeric_limits<Value>::infinity() :
|
669 |
640 |
std::numeric_limits<Value>::max())
|
670 |
641 |
{
|
671 |
642 |
// Check the value types
|
672 |
643 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
|
673 |
644 |
"The flow type of NetworkSimplex must be signed");
|
674 |
645 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
|
675 |
646 |
"The cost type of NetworkSimplex must be signed");
|
676 |
647 |
|
677 |
648 |
// Resize vectors
|
678 |
649 |
_node_num = countNodes(_graph);
|
679 |
650 |
_arc_num = countArcs(_graph);
|
680 |
651 |
int all_node_num = _node_num + 1;
|
681 |
|
int all_arc_num = _arc_num + _node_num;
|
|
652 |
int max_arc_num = _arc_num + 2 * _node_num;
|
682 |
653 |
|
683 |
|
_source.resize(all_arc_num);
|
684 |
|
_target.resize(all_arc_num);
|
|
654 |
_source.resize(max_arc_num);
|
|
655 |
_target.resize(max_arc_num);
|
685 |
656 |
|
686 |
|
_lower.resize(all_arc_num);
|
687 |
|
_upper.resize(all_arc_num);
|
688 |
|
_cap.resize(all_arc_num);
|
689 |
|
_cost.resize(all_arc_num);
|
|
657 |
_lower.resize(_arc_num);
|
|
658 |
_upper.resize(_arc_num);
|
|
659 |
_cap.resize(max_arc_num);
|
|
660 |
_cost.resize(max_arc_num);
|
690 |
661 |
_supply.resize(all_node_num);
|
691 |
|
_flow.resize(all_arc_num);
|
|
662 |
_flow.resize(max_arc_num);
|
692 |
663 |
_pi.resize(all_node_num);
|
693 |
664 |
|
694 |
665 |
_parent.resize(all_node_num);
|
695 |
666 |
_pred.resize(all_node_num);
|
696 |
667 |
_forward.resize(all_node_num);
|
697 |
668 |
_thread.resize(all_node_num);
|
698 |
669 |
_rev_thread.resize(all_node_num);
|
699 |
670 |
_succ_num.resize(all_node_num);
|
700 |
671 |
_last_succ.resize(all_node_num);
|
701 |
|
_state.resize(all_arc_num);
|
|
672 |
_state.resize(max_arc_num);
|
702 |
673 |
|
703 |
674 |
// Copy the graph (store the arcs in a mixed order)
|
704 |
675 |
int i = 0;
|
705 |
676 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
|
706 |
677 |
_node_id[n] = i;
|
707 |
678 |
}
|
708 |
679 |
int k = std::max(int(std::sqrt(double(_arc_num))), 10);
|
709 |
680 |
i = 0;
|
710 |
681 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
711 |
682 |
_arc_id[a] = i;
|
712 |
683 |
_source[i] = _node_id[_graph.source(a)];
|
713 |
684 |
_target[i] = _node_id[_graph.target(a)];
|
714 |
685 |
if ((i += k) >= _arc_num) i = (i % k) + 1;
|
715 |
686 |
}
|
716 |
687 |
|
717 |
688 |
// Initialize maps
|
718 |
689 |
for (int i = 0; i != _node_num; ++i) {
|
719 |
690 |
_supply[i] = 0;
|
720 |
691 |
}
|
721 |
692 |
for (int i = 0; i != _arc_num; ++i) {
|
722 |
693 |
_lower[i] = 0;
|
723 |
694 |
_upper[i] = INF;
|
724 |
695 |
_cost[i] = 1;
|
725 |
696 |
}
|
726 |
697 |
_have_lower = false;
|
727 |
698 |
_stype = GEQ;
|
728 |
699 |
}
|
729 |
700 |
|
730 |
701 |
/// \name Parameters
|
731 |
702 |
/// The parameters of the algorithm can be specified using these
|
732 |
703 |
/// functions.
|
733 |
704 |
|
734 |
705 |
/// @{
|
735 |
706 |
|
736 |
707 |
/// \brief Set the lower bounds on the arcs.
|
737 |
708 |
///
|
738 |
709 |
/// This function sets the lower bounds on the arcs.
|
739 |
710 |
/// If it is not used before calling \ref run(), the lower bounds
|
740 |
711 |
/// will be set to zero on all arcs.
|
741 |
712 |
///
|
742 |
713 |
/// \param map An arc map storing the lower bounds.
|
743 |
714 |
/// Its \c Value type must be convertible to the \c Value type
|
744 |
715 |
/// of the algorithm.
|
745 |
716 |
///
|
746 |
717 |
/// \return <tt>(*this)</tt>
|
747 |
718 |
template <typename LowerMap>
|
748 |
719 |
NetworkSimplex& lowerMap(const LowerMap& map) {
|
749 |
720 |
_have_lower = true;
|
750 |
721 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
751 |
722 |
_lower[_arc_id[a]] = map[a];
|
752 |
723 |
}
|
753 |
724 |
return *this;
|
754 |
725 |
}
|
755 |
726 |
|
756 |
727 |
/// \brief Set the upper bounds (capacities) on the arcs.
|
757 |
728 |
///
|
758 |
729 |
/// This function sets the upper bounds (capacities) on the arcs.
|
759 |
730 |
/// If it is not used before calling \ref run(), the upper bounds
|
760 |
731 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be
|
761 |
732 |
/// unbounded from above on each arc).
|
762 |
733 |
///
|
763 |
734 |
/// \param map An arc map storing the upper bounds.
|
764 |
735 |
/// Its \c Value type must be convertible to the \c Value type
|
765 |
736 |
/// of the algorithm.
|
766 |
737 |
///
|
767 |
738 |
/// \return <tt>(*this)</tt>
|
768 |
739 |
template<typename UpperMap>
|
769 |
740 |
NetworkSimplex& upperMap(const UpperMap& map) {
|
770 |
741 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
771 |
742 |
_upper[_arc_id[a]] = map[a];
|
772 |
743 |
}
|
773 |
744 |
return *this;
|
774 |
745 |
}
|
775 |
746 |
|
776 |
747 |
/// \brief Set the costs of the arcs.
|
777 |
748 |
///
|
778 |
749 |
/// This function sets the costs of the arcs.
|
779 |
750 |
/// If it is not used before calling \ref run(), the costs
|
780 |
751 |
/// will be set to \c 1 on all arcs.
|
781 |
752 |
///
|
782 |
753 |
/// \param map An arc map storing the costs.
|
783 |
754 |
/// Its \c Value type must be convertible to the \c Cost type
|
784 |
755 |
/// of the algorithm.
|
785 |
756 |
///
|
786 |
757 |
/// \return <tt>(*this)</tt>
|
787 |
758 |
template<typename CostMap>
|
788 |
759 |
NetworkSimplex& costMap(const CostMap& map) {
|
789 |
760 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
790 |
761 |
_cost[_arc_id[a]] = map[a];
|
791 |
762 |
}
|
792 |
763 |
return *this;
|
793 |
764 |
}
|
794 |
765 |
|
795 |
766 |
/// \brief Set the supply values of the nodes.
|
796 |
767 |
///
|
797 |
768 |
/// This function sets the supply values of the nodes.
|
798 |
769 |
/// If neither this function nor \ref stSupply() is used before
|
799 |
770 |
/// calling \ref run(), the supply of each node will be set to zero.
|
800 |
771 |
/// (It makes sense only if non-zero lower bounds are given.)
|
801 |
772 |
///
|
802 |
773 |
/// \param map A node map storing the supply values.
|
803 |
774 |
/// Its \c Value type must be convertible to the \c Value type
|
804 |
775 |
/// of the algorithm.
|
805 |
776 |
///
|
806 |
777 |
/// \return <tt>(*this)</tt>
|
807 |
778 |
template<typename SupplyMap>
|
808 |
779 |
NetworkSimplex& supplyMap(const SupplyMap& map) {
|
809 |
780 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
810 |
781 |
_supply[_node_id[n]] = map[n];
|
811 |
782 |
}
|
812 |
783 |
return *this;
|
813 |
784 |
}
|
814 |
785 |
|
815 |
786 |
/// \brief Set single source and target nodes and a supply value.
|
816 |
787 |
///
|
817 |
788 |
/// This function sets a single source node and a single target node
|
818 |
789 |
/// and the required flow value.
|
819 |
790 |
/// If neither this function nor \ref supplyMap() is used before
|
820 |
791 |
/// calling \ref run(), the supply of each node will be set to zero.
|
821 |
792 |
/// (It makes sense only if non-zero lower bounds are given.)
|
822 |
793 |
///
|
823 |
794 |
/// Using this function has the same effect as using \ref supplyMap()
|
824 |
795 |
/// with such a map in which \c k is assigned to \c s, \c -k is
|
825 |
796 |
/// assigned to \c t and all other nodes have zero supply value.
|
826 |
797 |
///
|
827 |
798 |
/// \param s The source node.
|
828 |
799 |
/// \param t The target node.
|
829 |
800 |
/// \param k The required amount of flow from node \c s to node \c t
|
... |
... |
@@ -944,303 +915,395 @@
|
944 |
915 |
return *this;
|
945 |
916 |
}
|
946 |
917 |
|
947 |
918 |
/// @}
|
948 |
919 |
|
949 |
920 |
/// \name Query Functions
|
950 |
921 |
/// The results of the algorithm can be obtained using these
|
951 |
922 |
/// functions.\n
|
952 |
923 |
/// The \ref run() function must be called before using them.
|
953 |
924 |
|
954 |
925 |
/// @{
|
955 |
926 |
|
956 |
927 |
/// \brief Return the total cost of the found flow.
|
957 |
928 |
///
|
958 |
929 |
/// This function returns the total cost of the found flow.
|
959 |
930 |
/// Its complexity is O(e).
|
960 |
931 |
///
|
961 |
932 |
/// \note The return type of the function can be specified as a
|
962 |
933 |
/// template parameter. For example,
|
963 |
934 |
/// \code
|
964 |
935 |
/// ns.totalCost<double>();
|
965 |
936 |
/// \endcode
|
966 |
937 |
/// It is useful if the total cost cannot be stored in the \c Cost
|
967 |
938 |
/// type of the algorithm, which is the default return type of the
|
968 |
939 |
/// function.
|
969 |
940 |
///
|
970 |
941 |
/// \pre \ref run() must be called before using this function.
|
971 |
942 |
template <typename Number>
|
972 |
943 |
Number totalCost() const {
|
973 |
944 |
Number c = 0;
|
974 |
945 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
975 |
946 |
int i = _arc_id[a];
|
976 |
947 |
c += Number(_flow[i]) * Number(_cost[i]);
|
977 |
948 |
}
|
978 |
949 |
return c;
|
979 |
950 |
}
|
980 |
951 |
|
981 |
952 |
#ifndef DOXYGEN
|
982 |
953 |
Cost totalCost() const {
|
983 |
954 |
return totalCost<Cost>();
|
984 |
955 |
}
|
985 |
956 |
#endif
|
986 |
957 |
|
987 |
958 |
/// \brief Return the flow on the given arc.
|
988 |
959 |
///
|
989 |
960 |
/// This function returns the flow on the given arc.
|
990 |
961 |
///
|
991 |
962 |
/// \pre \ref run() must be called before using this function.
|
992 |
963 |
Value flow(const Arc& a) const {
|
993 |
964 |
return _flow[_arc_id[a]];
|
994 |
965 |
}
|
995 |
966 |
|
996 |
967 |
/// \brief Return the flow map (the primal solution).
|
997 |
968 |
///
|
998 |
969 |
/// This function copies the flow value on each arc into the given
|
999 |
970 |
/// map. The \c Value type of the algorithm must be convertible to
|
1000 |
971 |
/// the \c Value type of the map.
|
1001 |
972 |
///
|
1002 |
973 |
/// \pre \ref run() must be called before using this function.
|
1003 |
974 |
template <typename FlowMap>
|
1004 |
975 |
void flowMap(FlowMap &map) const {
|
1005 |
976 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
1006 |
977 |
map.set(a, _flow[_arc_id[a]]);
|
1007 |
978 |
}
|
1008 |
979 |
}
|
1009 |
980 |
|
1010 |
981 |
/// \brief Return the potential (dual value) of the given node.
|
1011 |
982 |
///
|
1012 |
983 |
/// This function returns the potential (dual value) of the
|
1013 |
984 |
/// given node.
|
1014 |
985 |
///
|
1015 |
986 |
/// \pre \ref run() must be called before using this function.
|
1016 |
987 |
Cost potential(const Node& n) const {
|
1017 |
988 |
return _pi[_node_id[n]];
|
1018 |
989 |
}
|
1019 |
990 |
|
1020 |
991 |
/// \brief Return the potential map (the dual solution).
|
1021 |
992 |
///
|
1022 |
993 |
/// This function copies the potential (dual value) of each node
|
1023 |
994 |
/// into the given map.
|
1024 |
995 |
/// The \c Cost type of the algorithm must be convertible to the
|
1025 |
996 |
/// \c Value type of the map.
|
1026 |
997 |
///
|
1027 |
998 |
/// \pre \ref run() must be called before using this function.
|
1028 |
999 |
template <typename PotentialMap>
|
1029 |
1000 |
void potentialMap(PotentialMap &map) const {
|
1030 |
1001 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
1031 |
1002 |
map.set(n, _pi[_node_id[n]]);
|
1032 |
1003 |
}
|
1033 |
1004 |
}
|
1034 |
1005 |
|
1035 |
1006 |
/// @}
|
1036 |
1007 |
|
1037 |
1008 |
private:
|
1038 |
1009 |
|
1039 |
1010 |
// Initialize internal data structures
|
1040 |
1011 |
bool init() {
|
1041 |
1012 |
if (_node_num == 0) return false;
|
1042 |
1013 |
|
1043 |
1014 |
// Check the sum of supply values
|
1044 |
1015 |
_sum_supply = 0;
|
1045 |
1016 |
for (int i = 0; i != _node_num; ++i) {
|
1046 |
1017 |
_sum_supply += _supply[i];
|
1047 |
1018 |
}
|
1048 |
1019 |
if ( !((_stype == GEQ && _sum_supply <= 0) ||
|
1049 |
1020 |
(_stype == LEQ && _sum_supply >= 0)) ) return false;
|
1050 |
1021 |
|
1051 |
1022 |
// Remove non-zero lower bounds
|
1052 |
1023 |
if (_have_lower) {
|
1053 |
1024 |
for (int i = 0; i != _arc_num; ++i) {
|
1054 |
1025 |
Value c = _lower[i];
|
1055 |
1026 |
if (c >= 0) {
|
1056 |
1027 |
_cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
|
1057 |
1028 |
} else {
|
1058 |
1029 |
_cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
|
1059 |
1030 |
}
|
1060 |
1031 |
_supply[_source[i]] -= c;
|
1061 |
1032 |
_supply[_target[i]] += c;
|
1062 |
1033 |
}
|
1063 |
1034 |
} else {
|
1064 |
1035 |
for (int i = 0; i != _arc_num; ++i) {
|
1065 |
1036 |
_cap[i] = _upper[i];
|
1066 |
1037 |
}
|
1067 |
1038 |
}
|
1068 |
1039 |
|
1069 |
1040 |
// Initialize artifical cost
|
1070 |
1041 |
Cost ART_COST;
|
1071 |
1042 |
if (std::numeric_limits<Cost>::is_exact) {
|
1072 |
|
ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
|
|
1043 |
ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
|
1073 |
1044 |
} else {
|
1074 |
1045 |
ART_COST = std::numeric_limits<Cost>::min();
|
1075 |
1046 |
for (int i = 0; i != _arc_num; ++i) {
|
1076 |
1047 |
if (_cost[i] > ART_COST) ART_COST = _cost[i];
|
1077 |
1048 |
}
|
1078 |
1049 |
ART_COST = (ART_COST + 1) * _node_num;
|
1079 |
1050 |
}
|
1080 |
1051 |
|
1081 |
1052 |
// Initialize arc maps
|
1082 |
1053 |
for (int i = 0; i != _arc_num; ++i) {
|
1083 |
1054 |
_flow[i] = 0;
|
1084 |
1055 |
_state[i] = STATE_LOWER;
|
1085 |
1056 |
}
|
1086 |
1057 |
|
1087 |
1058 |
// Set data for the artificial root node
|
1088 |
1059 |
_root = _node_num;
|
1089 |
1060 |
_parent[_root] = -1;
|
1090 |
1061 |
_pred[_root] = -1;
|
1091 |
1062 |
_thread[_root] = 0;
|
1092 |
1063 |
_rev_thread[0] = _root;
|
1093 |
1064 |
_succ_num[_root] = _node_num + 1;
|
1094 |
1065 |
_last_succ[_root] = _root - 1;
|
1095 |
1066 |
_supply[_root] = -_sum_supply;
|
1096 |
|
_pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
|
|
1067 |
_pi[_root] = 0;
|
1097 |
1068 |
|
1098 |
1069 |
// Add artificial arcs and initialize the spanning tree data structure
|
1099 |
|
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
1100 |
|
_parent[u] = _root;
|
1101 |
|
_pred[u] = e;
|
1102 |
|
_thread[u] = u + 1;
|
1103 |
|
_rev_thread[u + 1] = u;
|
1104 |
|
_succ_num[u] = 1;
|
1105 |
|
_last_succ[u] = u;
|
1106 |
|
_cost[e] = ART_COST;
|
1107 |
|
_cap[e] = INF;
|
1108 |
|
_state[e] = STATE_TREE;
|
1109 |
|
if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
|
1110 |
|
_flow[e] = _supply[u];
|
1111 |
|
_forward[u] = true;
|
1112 |
|
_pi[u] = -ART_COST + _pi[_root];
|
1113 |
|
} else {
|
1114 |
|
_flow[e] = -_supply[u];
|
1115 |
|
_forward[u] = false;
|
1116 |
|
_pi[u] = ART_COST + _pi[_root];
|
|
1070 |
if (_sum_supply == 0) {
|
|
1071 |
// EQ supply constraints
|
|
1072 |
_search_arc_num = _arc_num;
|
|
1073 |
_all_arc_num = _arc_num + _node_num;
|
|
1074 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
|
1075 |
_parent[u] = _root;
|
|
1076 |
_pred[u] = e;
|
|
1077 |
_thread[u] = u + 1;
|
|
1078 |
_rev_thread[u + 1] = u;
|
|
1079 |
_succ_num[u] = 1;
|
|
1080 |
_last_succ[u] = u;
|
|
1081 |
_cap[e] = INF;
|
|
1082 |
_state[e] = STATE_TREE;
|
|
1083 |
if (_supply[u] >= 0) {
|
|
1084 |
_forward[u] = true;
|
|
1085 |
_pi[u] = 0;
|
|
1086 |
_source[e] = u;
|
|
1087 |
_target[e] = _root;
|
|
1088 |
_flow[e] = _supply[u];
|
|
1089 |
_cost[e] = 0;
|
|
1090 |
} else {
|
|
1091 |
_forward[u] = false;
|
|
1092 |
_pi[u] = ART_COST;
|
|
1093 |
_source[e] = _root;
|
|
1094 |
_target[e] = u;
|
|
1095 |
_flow[e] = -_supply[u];
|
|
1096 |
_cost[e] = ART_COST;
|
|
1097 |
}
|
1117 |
1098 |
}
|
1118 |
1099 |
}
|
|
1100 |
else if (_sum_supply > 0) {
|
|
1101 |
// LEQ supply constraints
|
|
1102 |
_search_arc_num = _arc_num + _node_num;
|
|
1103 |
int f = _arc_num + _node_num;
|
|
1104 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
|
1105 |
_parent[u] = _root;
|
|
1106 |
_thread[u] = u + 1;
|
|
1107 |
_rev_thread[u + 1] = u;
|
|
1108 |
_succ_num[u] = 1;
|
|
1109 |
_last_succ[u] = u;
|
|
1110 |
if (_supply[u] >= 0) {
|
|
1111 |
_forward[u] = true;
|
|
1112 |
_pi[u] = 0;
|
|
1113 |
_pred[u] = e;
|
|
1114 |
_source[e] = u;
|
|
1115 |
_target[e] = _root;
|
|
1116 |
_cap[e] = INF;
|
|
1117 |
_flow[e] = _supply[u];
|
|
1118 |
_cost[e] = 0;
|
|
1119 |
_state[e] = STATE_TREE;
|
|
1120 |
} else {
|
|
1121 |
_forward[u] = false;
|
|
1122 |
_pi[u] = ART_COST;
|
|
1123 |
_pred[u] = f;
|
|
1124 |
_source[f] = _root;
|
|
1125 |
_target[f] = u;
|
|
1126 |
_cap[f] = INF;
|
|
1127 |
_flow[f] = -_supply[u];
|
|
1128 |
_cost[f] = ART_COST;
|
|
1129 |
_state[f] = STATE_TREE;
|
|
1130 |
_source[e] = u;
|
|
1131 |
_target[e] = _root;
|
|
1132 |
_cap[e] = INF;
|
|
1133 |
_flow[e] = 0;
|
|
1134 |
_cost[e] = 0;
|
|
1135 |
_state[e] = STATE_LOWER;
|
|
1136 |
++f;
|
|
1137 |
}
|
|
1138 |
}
|
|
1139 |
_all_arc_num = f;
|
|
1140 |
}
|
|
1141 |
else {
|
|
1142 |
// GEQ supply constraints
|
|
1143 |
_search_arc_num = _arc_num + _node_num;
|
|
1144 |
int f = _arc_num + _node_num;
|
|
1145 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
|
1146 |
_parent[u] = _root;
|
|
1147 |
_thread[u] = u + 1;
|
|
1148 |
_rev_thread[u + 1] = u;
|
|
1149 |
_succ_num[u] = 1;
|
|
1150 |
_last_succ[u] = u;
|
|
1151 |
if (_supply[u] <= 0) {
|
|
1152 |
_forward[u] = false;
|
|
1153 |
_pi[u] = 0;
|
|
1154 |
_pred[u] = e;
|
|
1155 |
_source[e] = _root;
|
|
1156 |
_target[e] = u;
|
|
1157 |
_cap[e] = INF;
|
|
1158 |
_flow[e] = -_supply[u];
|
|
1159 |
_cost[e] = 0;
|
|
1160 |
_state[e] = STATE_TREE;
|
|
1161 |
} else {
|
|
1162 |
_forward[u] = true;
|
|
1163 |
_pi[u] = -ART_COST;
|
|
1164 |
_pred[u] = f;
|
|
1165 |
_source[f] = u;
|
|
1166 |
_target[f] = _root;
|
|
1167 |
_cap[f] = INF;
|
|
1168 |
_flow[f] = _supply[u];
|
|
1169 |
_state[f] = STATE_TREE;
|
|
1170 |
_cost[f] = ART_COST;
|
|
1171 |
_source[e] = _root;
|
|
1172 |
_target[e] = u;
|
|
1173 |
_cap[e] = INF;
|
|
1174 |
_flow[e] = 0;
|
|
1175 |
_cost[e] = 0;
|
|
1176 |
_state[e] = STATE_LOWER;
|
|
1177 |
++f;
|
|
1178 |
}
|
|
1179 |
}
|
|
1180 |
_all_arc_num = f;
|
|
1181 |
}
|
1119 |
1182 |
|
1120 |
1183 |
return true;
|
1121 |
1184 |
}
|
1122 |
1185 |
|
1123 |
1186 |
// Find the join node
|
1124 |
1187 |
void findJoinNode() {
|
1125 |
1188 |
int u = _source[in_arc];
|
1126 |
1189 |
int v = _target[in_arc];
|
1127 |
1190 |
while (u != v) {
|
1128 |
1191 |
if (_succ_num[u] < _succ_num[v]) {
|
1129 |
1192 |
u = _parent[u];
|
1130 |
1193 |
} else {
|
1131 |
1194 |
v = _parent[v];
|
1132 |
1195 |
}
|
1133 |
1196 |
}
|
1134 |
1197 |
join = u;
|
1135 |
1198 |
}
|
1136 |
1199 |
|
1137 |
1200 |
// Find the leaving arc of the cycle and returns true if the
|
1138 |
1201 |
// leaving arc is not the same as the entering arc
|
1139 |
1202 |
bool findLeavingArc() {
|
1140 |
1203 |
// Initialize first and second nodes according to the direction
|
1141 |
1204 |
// of the cycle
|
1142 |
1205 |
if (_state[in_arc] == STATE_LOWER) {
|
1143 |
1206 |
first = _source[in_arc];
|
1144 |
1207 |
second = _target[in_arc];
|
1145 |
1208 |
} else {
|
1146 |
1209 |
first = _target[in_arc];
|
1147 |
1210 |
second = _source[in_arc];
|
1148 |
1211 |
}
|
1149 |
1212 |
delta = _cap[in_arc];
|
1150 |
1213 |
int result = 0;
|
1151 |
1214 |
Value d;
|
1152 |
1215 |
int e;
|
1153 |
1216 |
|
1154 |
1217 |
// Search the cycle along the path form the first node to the root
|
1155 |
1218 |
for (int u = first; u != join; u = _parent[u]) {
|
1156 |
1219 |
e = _pred[u];
|
1157 |
1220 |
d = _forward[u] ?
|
1158 |
1221 |
_flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
|
1159 |
1222 |
if (d < delta) {
|
1160 |
1223 |
delta = d;
|
1161 |
1224 |
u_out = u;
|
1162 |
1225 |
result = 1;
|
1163 |
1226 |
}
|
1164 |
1227 |
}
|
1165 |
1228 |
// Search the cycle along the path form the second node to the root
|
1166 |
1229 |
for (int u = second; u != join; u = _parent[u]) {
|
1167 |
1230 |
e = _pred[u];
|
1168 |
1231 |
d = _forward[u] ?
|
1169 |
1232 |
(_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
|
1170 |
1233 |
if (d <= delta) {
|
1171 |
1234 |
delta = d;
|
1172 |
1235 |
u_out = u;
|
1173 |
1236 |
result = 2;
|
1174 |
1237 |
}
|
1175 |
1238 |
}
|
1176 |
1239 |
|
1177 |
1240 |
if (result == 1) {
|
1178 |
1241 |
u_in = first;
|
1179 |
1242 |
v_in = second;
|
1180 |
1243 |
} else {
|
1181 |
1244 |
u_in = second;
|
1182 |
1245 |
v_in = first;
|
1183 |
1246 |
}
|
1184 |
1247 |
return result != 0;
|
1185 |
1248 |
}
|
1186 |
1249 |
|
1187 |
1250 |
// Change _flow and _state vectors
|
1188 |
1251 |
void changeFlow(bool change) {
|
1189 |
1252 |
// Augment along the cycle
|
1190 |
1253 |
if (delta > 0) {
|
1191 |
1254 |
Value val = _state[in_arc] * delta;
|
1192 |
1255 |
_flow[in_arc] += val;
|
1193 |
1256 |
for (int u = _source[in_arc]; u != join; u = _parent[u]) {
|
1194 |
1257 |
_flow[_pred[u]] += _forward[u] ? -val : val;
|
1195 |
1258 |
}
|
1196 |
1259 |
for (int u = _target[in_arc]; u != join; u = _parent[u]) {
|
1197 |
1260 |
_flow[_pred[u]] += _forward[u] ? val : -val;
|
1198 |
1261 |
}
|
1199 |
1262 |
}
|
1200 |
1263 |
// Update the state of the entering and leaving arcs
|
1201 |
1264 |
if (change) {
|
1202 |
1265 |
_state[in_arc] = STATE_TREE;
|
1203 |
1266 |
_state[_pred[u_out]] =
|
1204 |
1267 |
(_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
|
1205 |
1268 |
} else {
|
1206 |
1269 |
_state[in_arc] = -_state[in_arc];
|
1207 |
1270 |
}
|
1208 |
1271 |
}
|
1209 |
1272 |
|
1210 |
1273 |
// Update the tree structure
|
1211 |
1274 |
void updateTreeStructure() {
|
1212 |
1275 |
int u, w;
|
1213 |
1276 |
int old_rev_thread = _rev_thread[u_out];
|
1214 |
1277 |
int old_succ_num = _succ_num[u_out];
|
1215 |
1278 |
int old_last_succ = _last_succ[u_out];
|
1216 |
1279 |
v_out = _parent[u_out];
|
1217 |
1280 |
|
1218 |
1281 |
u = _last_succ[u_in]; // the last successor of u_in
|
1219 |
1282 |
right = _thread[u]; // the node after it
|
1220 |
1283 |
|
1221 |
1284 |
// Handle the case when old_rev_thread equals to v_in
|
1222 |
1285 |
// (it also means that join and v_out coincide)
|
1223 |
1286 |
if (old_rev_thread == v_in) {
|
1224 |
1287 |
last = _thread[_last_succ[u_out]];
|
1225 |
1288 |
} else {
|
1226 |
1289 |
last = _thread[v_in];
|
1227 |
1290 |
}
|
1228 |
1291 |
|
1229 |
1292 |
// Update _thread and _parent along the stem nodes (i.e. the nodes
|
1230 |
1293 |
// between u_in and u_out, whose parent have to be changed)
|
1231 |
1294 |
_thread[v_in] = stem = u_in;
|
1232 |
1295 |
_dirty_revs.clear();
|
1233 |
1296 |
_dirty_revs.push_back(v_in);
|
1234 |
1297 |
par_stem = v_in;
|
1235 |
1298 |
while (stem != u_out) {
|
1236 |
1299 |
// Insert the next stem node into the thread list
|
1237 |
1300 |
new_stem = _parent[stem];
|
1238 |
1301 |
_thread[u] = new_stem;
|
1239 |
1302 |
_dirty_revs.push_back(u);
|
1240 |
1303 |
|
1241 |
1304 |
// Remove the subtree of stem from the thread list
|
1242 |
1305 |
w = _rev_thread[stem];
|
1243 |
1306 |
_thread[w] = right;
|
1244 |
1307 |
_rev_thread[right] = w;
|
1245 |
1308 |
|
1246 |
1309 |
// Change the parent node and shift stem nodes
|
... |
... |
@@ -1249,166 +1312,178 @@
|
1249 |
1312 |
stem = new_stem;
|
1250 |
1313 |
|
1251 |
1314 |
// Update u and right
|
1252 |
1315 |
u = _last_succ[stem] == _last_succ[par_stem] ?
|
1253 |
1316 |
_rev_thread[par_stem] : _last_succ[stem];
|
1254 |
1317 |
right = _thread[u];
|
1255 |
1318 |
}
|
1256 |
1319 |
_parent[u_out] = par_stem;
|
1257 |
1320 |
_thread[u] = last;
|
1258 |
1321 |
_rev_thread[last] = u;
|
1259 |
1322 |
_last_succ[u_out] = u;
|
1260 |
1323 |
|
1261 |
1324 |
// Remove the subtree of u_out from the thread list except for
|
1262 |
1325 |
// the case when old_rev_thread equals to v_in
|
1263 |
1326 |
// (it also means that join and v_out coincide)
|
1264 |
1327 |
if (old_rev_thread != v_in) {
|
1265 |
1328 |
_thread[old_rev_thread] = right;
|
1266 |
1329 |
_rev_thread[right] = old_rev_thread;
|
1267 |
1330 |
}
|
1268 |
1331 |
|
1269 |
1332 |
// Update _rev_thread using the new _thread values
|
1270 |
1333 |
for (int i = 0; i < int(_dirty_revs.size()); ++i) {
|
1271 |
1334 |
u = _dirty_revs[i];
|
1272 |
1335 |
_rev_thread[_thread[u]] = u;
|
1273 |
1336 |
}
|
1274 |
1337 |
|
1275 |
1338 |
// Update _pred, _forward, _last_succ and _succ_num for the
|
1276 |
1339 |
// stem nodes from u_out to u_in
|
1277 |
1340 |
int tmp_sc = 0, tmp_ls = _last_succ[u_out];
|
1278 |
1341 |
u = u_out;
|
1279 |
1342 |
while (u != u_in) {
|
1280 |
1343 |
w = _parent[u];
|
1281 |
1344 |
_pred[u] = _pred[w];
|
1282 |
1345 |
_forward[u] = !_forward[w];
|
1283 |
1346 |
tmp_sc += _succ_num[u] - _succ_num[w];
|
1284 |
1347 |
_succ_num[u] = tmp_sc;
|
1285 |
1348 |
_last_succ[w] = tmp_ls;
|
1286 |
1349 |
u = w;
|
1287 |
1350 |
}
|
1288 |
1351 |
_pred[u_in] = in_arc;
|
1289 |
1352 |
_forward[u_in] = (u_in == _source[in_arc]);
|
1290 |
1353 |
_succ_num[u_in] = old_succ_num;
|
1291 |
1354 |
|
1292 |
1355 |
// Set limits for updating _last_succ form v_in and v_out
|
1293 |
1356 |
// towards the root
|
1294 |
1357 |
int up_limit_in = -1;
|
1295 |
1358 |
int up_limit_out = -1;
|
1296 |
1359 |
if (_last_succ[join] == v_in) {
|
1297 |
1360 |
up_limit_out = join;
|
1298 |
1361 |
} else {
|
1299 |
1362 |
up_limit_in = join;
|
1300 |
1363 |
}
|
1301 |
1364 |
|
1302 |
1365 |
// Update _last_succ from v_in towards the root
|
1303 |
1366 |
for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
|
1304 |
1367 |
u = _parent[u]) {
|
1305 |
1368 |
_last_succ[u] = _last_succ[u_out];
|
1306 |
1369 |
}
|
1307 |
1370 |
// Update _last_succ from v_out towards the root
|
1308 |
1371 |
if (join != old_rev_thread && v_in != old_rev_thread) {
|
1309 |
1372 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
|
1310 |
1373 |
u = _parent[u]) {
|
1311 |
1374 |
_last_succ[u] = old_rev_thread;
|
1312 |
1375 |
}
|
1313 |
1376 |
} else {
|
1314 |
1377 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
|
1315 |
1378 |
u = _parent[u]) {
|
1316 |
1379 |
_last_succ[u] = _last_succ[u_out];
|
1317 |
1380 |
}
|
1318 |
1381 |
}
|
1319 |
1382 |
|
1320 |
1383 |
// Update _succ_num from v_in to join
|
1321 |
1384 |
for (u = v_in; u != join; u = _parent[u]) {
|
1322 |
1385 |
_succ_num[u] += old_succ_num;
|
1323 |
1386 |
}
|
1324 |
1387 |
// Update _succ_num from v_out to join
|
1325 |
1388 |
for (u = v_out; u != join; u = _parent[u]) {
|
1326 |
1389 |
_succ_num[u] -= old_succ_num;
|
1327 |
1390 |
}
|
1328 |
1391 |
}
|
1329 |
1392 |
|
1330 |
1393 |
// Update potentials
|
1331 |
1394 |
void updatePotential() {
|
1332 |
1395 |
Cost sigma = _forward[u_in] ?
|
1333 |
1396 |
_pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
|
1334 |
1397 |
_pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
|
1335 |
1398 |
// Update potentials in the subtree, which has been moved
|
1336 |
1399 |
int end = _thread[_last_succ[u_in]];
|
1337 |
1400 |
for (int u = u_in; u != end; u = _thread[u]) {
|
1338 |
1401 |
_pi[u] += sigma;
|
1339 |
1402 |
}
|
1340 |
1403 |
}
|
1341 |
1404 |
|
1342 |
1405 |
// Execute the algorithm
|
1343 |
1406 |
ProblemType start(PivotRule pivot_rule) {
|
1344 |
1407 |
// Select the pivot rule implementation
|
1345 |
1408 |
switch (pivot_rule) {
|
1346 |
1409 |
case FIRST_ELIGIBLE:
|
1347 |
1410 |
return start<FirstEligiblePivotRule>();
|
1348 |
1411 |
case BEST_ELIGIBLE:
|
1349 |
1412 |
return start<BestEligiblePivotRule>();
|
1350 |
1413 |
case BLOCK_SEARCH:
|
1351 |
1414 |
return start<BlockSearchPivotRule>();
|
1352 |
1415 |
case CANDIDATE_LIST:
|
1353 |
1416 |
return start<CandidateListPivotRule>();
|
1354 |
1417 |
case ALTERING_LIST:
|
1355 |
1418 |
return start<AlteringListPivotRule>();
|
1356 |
1419 |
}
|
1357 |
1420 |
return INFEASIBLE; // avoid warning
|
1358 |
1421 |
}
|
1359 |
1422 |
|
1360 |
1423 |
template <typename PivotRuleImpl>
|
1361 |
1424 |
ProblemType start() {
|
1362 |
1425 |
PivotRuleImpl pivot(*this);
|
1363 |
1426 |
|
1364 |
1427 |
// Execute the Network Simplex algorithm
|
1365 |
1428 |
while (pivot.findEnteringArc()) {
|
1366 |
1429 |
findJoinNode();
|
1367 |
1430 |
bool change = findLeavingArc();
|
1368 |
1431 |
if (delta >= INF) return UNBOUNDED;
|
1369 |
1432 |
changeFlow(change);
|
1370 |
1433 |
if (change) {
|
1371 |
1434 |
updateTreeStructure();
|
1372 |
1435 |
updatePotential();
|
1373 |
1436 |
}
|
1374 |
1437 |
}
|
1375 |
1438 |
|
1376 |
1439 |
// Check feasibility
|
1377 |
|
if (_sum_supply < 0) {
|
1378 |
|
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
1379 |
|
if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
|
1380 |
|
}
|
1381 |
|
}
|
1382 |
|
else if (_sum_supply > 0) {
|
1383 |
|
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
1384 |
|
if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
|
1385 |
|
}
|
1386 |
|
}
|
1387 |
|
else {
|
1388 |
|
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
1389 |
|
if (_flow[e] != 0) return INFEASIBLE;
|
1390 |
|
}
|
|
1440 |
for (int e = _search_arc_num; e != _all_arc_num; ++e) {
|
|
1441 |
if (_flow[e] != 0) return INFEASIBLE;
|
1391 |
1442 |
}
|
1392 |
1443 |
|
1393 |
1444 |
// Transform the solution and the supply map to the original form
|
1394 |
1445 |
if (_have_lower) {
|
1395 |
1446 |
for (int i = 0; i != _arc_num; ++i) {
|
1396 |
1447 |
Value c = _lower[i];
|
1397 |
1448 |
if (c != 0) {
|
1398 |
1449 |
_flow[i] += c;
|
1399 |
1450 |
_supply[_source[i]] += c;
|
1400 |
1451 |
_supply[_target[i]] -= c;
|
1401 |
1452 |
}
|
1402 |
1453 |
}
|
1403 |
1454 |
}
|
|
1455 |
|
|
1456 |
// Shift potentials to meet the requirements of the GEQ/LEQ type
|
|
1457 |
// optimality conditions
|
|
1458 |
if (_sum_supply == 0) {
|
|
1459 |
if (_stype == GEQ) {
|
|
1460 |
Cost max_pot = std::numeric_limits<Cost>::min();
|
|
1461 |
for (int i = 0; i != _node_num; ++i) {
|
|
1462 |
if (_pi[i] > max_pot) max_pot = _pi[i];
|
|
1463 |
}
|
|
1464 |
if (max_pot > 0) {
|
|
1465 |
for (int i = 0; i != _node_num; ++i)
|
|
1466 |
_pi[i] -= max_pot;
|
|
1467 |
}
|
|
1468 |
} else {
|
|
1469 |
Cost min_pot = std::numeric_limits<Cost>::max();
|
|
1470 |
for (int i = 0; i != _node_num; ++i) {
|
|
1471 |
if (_pi[i] < min_pot) min_pot = _pi[i];
|
|
1472 |
}
|
|
1473 |
if (min_pot < 0) {
|
|
1474 |
for (int i = 0; i != _node_num; ++i)
|
|
1475 |
_pi[i] -= min_pot;
|
|
1476 |
}
|
|
1477 |
}
|
|
1478 |
}
|
1404 |
1479 |
|
1405 |
1480 |
return OPTIMAL;
|
1406 |
1481 |
}
|
1407 |
1482 |
|
1408 |
1483 |
}; //class NetworkSimplex
|
1409 |
1484 |
|
1410 |
1485 |
///@}
|
1411 |
1486 |
|
1412 |
1487 |
} //namespace lemon
|
1413 |
1488 |
|
1414 |
1489 |
#endif //LEMON_NETWORK_SIMPLEX_H
|