| 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 |
22 |
/// \ingroup min_cost_flow
|
| 23 |
23 |
///
|
| 24 |
24 |
/// \file
|
| 25 |
|
/// \brief Network simplex algorithm for finding a minimum cost flow.
|
|
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 |
36 |
/// \addtogroup min_cost_flow
|
| 37 |
37 |
/// @{
|
| 38 |
38 |
|
| 39 |
|
/// \brief Implementation of the primal network simplex algorithm
|
|
39 |
/// \brief Implementation of the primal Network Simplex algorithm
|
| 40 |
40 |
/// for finding a \ref min_cost_flow "minimum cost flow".
|
| 41 |
41 |
///
|
| 42 |
|
/// \ref NetworkSimplex implements the primal network simplex algorithm
|
|
42 |
/// \ref NetworkSimplex implements the primal Network Simplex algorithm
|
| 43 |
43 |
/// for finding a \ref min_cost_flow "minimum cost flow".
|
| 44 |
44 |
///
|
| 45 |
|
/// \tparam Digraph The digraph type the algorithm runs on.
|
| 46 |
|
/// \tparam LowerMap The type of the lower bound map.
|
| 47 |
|
/// \tparam CapacityMap The type of the capacity (upper bound) map.
|
| 48 |
|
/// \tparam CostMap The type of the cost (length) map.
|
| 49 |
|
/// \tparam SupplyMap The type of the supply map.
|
|
45 |
/// \tparam GR The digraph type the algorithm runs on.
|
|
46 |
/// \tparam V The value type used in the algorithm.
|
|
47 |
/// By default it is \c int.
|
| 50 |
48 |
///
|
| 51 |
|
/// \warning
|
| 52 |
|
/// - Arc capacities and costs should be \e non-negative \e integers.
|
| 53 |
|
/// - Supply values should be \e signed \e integers.
|
| 54 |
|
/// - The value types of the maps should be convertible to each other.
|
| 55 |
|
/// - \c CostMap::Value must be signed type.
|
|
49 |
/// \warning \c V must be a signed integer type.
|
| 56 |
50 |
///
|
| 57 |
|
/// \note \ref NetworkSimplex provides five different pivot rule
|
| 58 |
|
/// implementations that significantly affect the efficiency of the
|
| 59 |
|
/// algorithm.
|
| 60 |
|
/// By default "Block Search" pivot rule is used, which proved to be
|
| 61 |
|
/// by far the most efficient according to our benchmark tests.
|
| 62 |
|
/// However another pivot rule can be selected using \ref run()
|
| 63 |
|
/// function with the proper parameter.
|
| 64 |
|
#ifdef DOXYGEN
|
| 65 |
|
template < typename Digraph,
|
| 66 |
|
typename LowerMap,
|
| 67 |
|
typename CapacityMap,
|
| 68 |
|
typename CostMap,
|
| 69 |
|
typename SupplyMap >
|
| 70 |
|
|
| 71 |
|
#else
|
| 72 |
|
template < typename Digraph,
|
| 73 |
|
typename LowerMap = typename Digraph::template ArcMap<int>,
|
| 74 |
|
typename CapacityMap = typename Digraph::template ArcMap<int>,
|
| 75 |
|
typename CostMap = typename Digraph::template ArcMap<int>,
|
| 76 |
|
typename SupplyMap = typename Digraph::template NodeMap<int> >
|
| 77 |
|
#endif
|
|
51 |
/// \note %NetworkSimplex provides five different pivot rule
|
|
52 |
/// implementations. For more information see \ref PivotRule.
|
|
53 |
template <typename GR, typename V = int>
|
| 78 |
54 |
class NetworkSimplex
|
| 79 |
55 |
{
|
| 80 |
|
TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
|
|
56 |
public:
|
| 81 |
57 |
|
| 82 |
|
typedef typename CapacityMap::Value Capacity;
|
| 83 |
|
typedef typename CostMap::Value Cost;
|
| 84 |
|
typedef typename SupplyMap::Value Supply;
|
|
58 |
/// The value type of the algorithm
|
|
59 |
typedef V Value;
|
|
60 |
/// The type of the flow map
|
|
61 |
typedef typename GR::template ArcMap<Value> FlowMap;
|
|
62 |
/// The type of the potential map
|
|
63 |
typedef typename GR::template NodeMap<Value> PotentialMap;
|
|
64 |
|
|
65 |
public:
|
|
66 |
|
|
67 |
/// \brief Enum type for selecting the pivot rule.
|
|
68 |
///
|
|
69 |
/// Enum type for selecting the pivot rule for the \ref run()
|
|
70 |
/// function.
|
|
71 |
///
|
|
72 |
/// \ref NetworkSimplex provides five different pivot rule
|
|
73 |
/// implementations that significantly affect the running time
|
|
74 |
/// of the algorithm.
|
|
75 |
/// By default \ref BLOCK_SEARCH "Block Search" is used, which
|
|
76 |
/// proved to be the most efficient and the most robust on various
|
|
77 |
/// test inputs according to our benchmark tests.
|
|
78 |
/// However another pivot rule can be selected using the \ref run()
|
|
79 |
/// function with the proper parameter.
|
|
80 |
enum PivotRule {
|
|
81 |
|
|
82 |
/// The First Eligible pivot rule.
|
|
83 |
/// The next eligible arc is selected in a wraparound fashion
|
|
84 |
/// in every iteration.
|
|
85 |
FIRST_ELIGIBLE,
|
|
86 |
|
|
87 |
/// The Best Eligible pivot rule.
|
|
88 |
/// The best eligible arc is selected in every iteration.
|
|
89 |
BEST_ELIGIBLE,
|
|
90 |
|
|
91 |
/// The Block Search pivot rule.
|
|
92 |
/// A specified number of arcs are examined in every iteration
|
|
93 |
/// in a wraparound fashion and the best eligible arc is selected
|
|
94 |
/// from this block.
|
|
95 |
BLOCK_SEARCH,
|
|
96 |
|
|
97 |
/// The Candidate List pivot rule.
|
|
98 |
/// In a major iteration a candidate list is built from eligible arcs
|
|
99 |
/// in a wraparound fashion and in the following minor iterations
|
|
100 |
/// the best eligible arc is selected from this list.
|
|
101 |
CANDIDATE_LIST,
|
|
102 |
|
|
103 |
/// The Altering Candidate List pivot rule.
|
|
104 |
/// It is a modified version of the Candidate List method.
|
|
105 |
/// It keeps only the several best eligible arcs from the former
|
|
106 |
/// candidate list and extends this list in every iteration.
|
|
107 |
ALTERING_LIST
|
|
108 |
};
|
|
109 |
|
|
110 |
private:
|
|
111 |
|
|
112 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR);
|
|
113 |
|
|
114 |
typedef typename GR::template ArcMap<Value> ValueArcMap;
|
|
115 |
typedef typename GR::template NodeMap<Value> ValueNodeMap;
|
| 85 |
116 |
|
| 86 |
117 |
typedef std::vector<Arc> ArcVector;
|
| 87 |
118 |
typedef std::vector<Node> NodeVector;
|
| 88 |
119 |
typedef std::vector<int> IntVector;
|
| 89 |
120 |
typedef std::vector<bool> BoolVector;
|
| 90 |
|
typedef std::vector<Capacity> CapacityVector;
|
| 91 |
|
typedef std::vector<Cost> CostVector;
|
| 92 |
|
typedef std::vector<Supply> SupplyVector;
|
| 93 |
|
|
| 94 |
|
public:
|
| 95 |
|
|
| 96 |
|
/// The type of the flow map
|
| 97 |
|
typedef typename Digraph::template ArcMap<Capacity> FlowMap;
|
| 98 |
|
/// The type of the potential map
|
| 99 |
|
typedef typename Digraph::template NodeMap<Cost> PotentialMap;
|
| 100 |
|
|
| 101 |
|
public:
|
| 102 |
|
|
| 103 |
|
/// Enum type for selecting the pivot rule used by \ref run()
|
| 104 |
|
enum PivotRuleEnum {
|
| 105 |
|
FIRST_ELIGIBLE_PIVOT,
|
| 106 |
|
BEST_ELIGIBLE_PIVOT,
|
| 107 |
|
BLOCK_SEARCH_PIVOT,
|
| 108 |
|
CANDIDATE_LIST_PIVOT,
|
| 109 |
|
ALTERING_LIST_PIVOT
|
| 110 |
|
};
|
| 111 |
|
|
| 112 |
|
private:
|
|
121 |
typedef std::vector<Value> ValueVector;
|
| 113 |
122 |
|
| 114 |
123 |
// State constants for arcs
|
| 115 |
124 |
enum ArcStateEnum {
|
| 116 |
125 |
STATE_UPPER = -1,
|
| 117 |
126 |
STATE_TREE = 0,
|
| 118 |
127 |
STATE_LOWER = 1
|
| 119 |
128 |
};
|
| 120 |
129 |
|
| 121 |
130 |
private:
|
| 122 |
131 |
|
| 123 |
|
// References for the original data
|
| 124 |
|
const Digraph &_graph;
|
| 125 |
|
const LowerMap *_orig_lower;
|
| 126 |
|
const CapacityMap &_orig_cap;
|
| 127 |
|
const CostMap &_orig_cost;
|
| 128 |
|
const SupplyMap *_orig_supply;
|
| 129 |
|
Node _orig_source;
|
| 130 |
|
Node _orig_target;
|
| 131 |
|
Capacity _orig_flow_value;
|
|
132 |
// Data related to the underlying digraph
|
|
133 |
const GR &_graph;
|
|
134 |
int _node_num;
|
|
135 |
int _arc_num;
|
|
136 |
|
|
137 |
// Parameters of the problem
|
|
138 |
ValueArcMap *_plower;
|
|
139 |
ValueArcMap *_pupper;
|
|
140 |
ValueArcMap *_pcost;
|
|
141 |
ValueNodeMap *_psupply;
|
|
142 |
bool _pstsup;
|
|
143 |
Node _psource, _ptarget;
|
|
144 |
Value _pstflow;
|
| 132 |
145 |
|
| 133 |
146 |
// Result maps
|
| 134 |
147 |
FlowMap *_flow_map;
|
| 135 |
148 |
PotentialMap *_potential_map;
|
| 136 |
149 |
bool _local_flow;
|
| 137 |
150 |
bool _local_potential;
|
| 138 |
151 |
|
| 139 |
|
// The number of nodes and arcs in the original graph
|
| 140 |
|
int _node_num;
|
| 141 |
|
int _arc_num;
|
| 142 |
|
|
| 143 |
|
// Data structures for storing the graph
|
|
152 |
// Data structures for storing the digraph
|
| 144 |
153 |
IntNodeMap _node_id;
|
| 145 |
154 |
ArcVector _arc_ref;
|
| 146 |
155 |
IntVector _source;
|
| 147 |
156 |
IntVector _target;
|
| 148 |
157 |
|
| 149 |
|
// Node and arc maps
|
| 150 |
|
CapacityVector _cap;
|
| 151 |
|
CostVector _cost;
|
| 152 |
|
CostVector _supply;
|
| 153 |
|
CapacityVector _flow;
|
| 154 |
|
CostVector _pi;
|
|
158 |
// Node and arc data
|
|
159 |
ValueVector _cap;
|
|
160 |
ValueVector _cost;
|
|
161 |
ValueVector _supply;
|
|
162 |
ValueVector _flow;
|
|
163 |
ValueVector _pi;
|
| 155 |
164 |
|
| 156 |
165 |
// Data for storing the spanning tree structure
|
| 157 |
166 |
IntVector _parent;
|
| 158 |
167 |
IntVector _pred;
|
| 159 |
168 |
IntVector _thread;
|
| 160 |
169 |
IntVector _rev_thread;
|
| 161 |
170 |
IntVector _succ_num;
|
| 162 |
171 |
IntVector _last_succ;
|
| 163 |
172 |
IntVector _dirty_revs;
|
| 164 |
173 |
BoolVector _forward;
|
| 165 |
174 |
IntVector _state;
|
| 166 |
175 |
int _root;
|
| 167 |
176 |
|
| 168 |
177 |
// Temporary data used in the current pivot iteration
|
| 169 |
178 |
int in_arc, join, u_in, v_in, u_out, v_out;
|
| 170 |
179 |
int first, second, right, last;
|
| 171 |
180 |
int stem, par_stem, new_stem;
|
| 172 |
|
Capacity delta;
|
|
181 |
Value delta;
|
| 173 |
182 |
|
| 174 |
183 |
private:
|
| 175 |
184 |
|
| 176 |
|
/// \brief Implementation of the "First Eligible" pivot rule for the
|
| 177 |
|
/// \ref NetworkSimplex "network simplex" algorithm.
|
| 178 |
|
///
|
| 179 |
|
/// This class implements the "First Eligible" pivot rule
|
| 180 |
|
/// for the \ref NetworkSimplex "network simplex" algorithm.
|
| 181 |
|
///
|
| 182 |
|
/// For more information see \ref NetworkSimplex::run().
|
|
185 |
// Implementation of the First Eligible pivot rule
|
| 183 |
186 |
class FirstEligiblePivotRule
|
| 184 |
187 |
{
|
| 185 |
188 |
private:
|
| 186 |
189 |
|
| 187 |
190 |
// References to the NetworkSimplex class
|
| 188 |
191 |
const IntVector &_source;
|
| 189 |
192 |
const IntVector &_target;
|
| 190 |
|
const CostVector &_cost;
|
|
193 |
const ValueVector &_cost;
|
| 191 |
194 |
const IntVector &_state;
|
| 192 |
|
const CostVector &_pi;
|
|
195 |
const ValueVector &_pi;
|
| 193 |
196 |
int &_in_arc;
|
| 194 |
197 |
int _arc_num;
|
| 195 |
198 |
|
| 196 |
199 |
// Pivot rule data
|
| 197 |
200 |
int _next_arc;
|
| 198 |
201 |
|
| 199 |
202 |
public:
|
| 200 |
203 |
|
| 201 |
|
/// Constructor
|
|
204 |
// Constructor
|
| 202 |
205 |
FirstEligiblePivotRule(NetworkSimplex &ns) :
|
| 203 |
206 |
_source(ns._source), _target(ns._target),
|
| 204 |
207 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 205 |
208 |
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
| 206 |
209 |
{}
|
| 207 |
210 |
|
| 208 |
|
/// Find next entering arc
|
|
211 |
// Find next entering arc
|
| 209 |
212 |
bool findEnteringArc() {
|
| 210 |
|
Cost c;
|
|
213 |
Value c;
|
| 211 |
214 |
for (int e = _next_arc; e < _arc_num; ++e) {
|
| 212 |
215 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 213 |
216 |
if (c < 0) {
|
| 214 |
217 |
_in_arc = e;
|
| 215 |
218 |
_next_arc = e + 1;
|
| 216 |
219 |
return true;
|
| 217 |
220 |
}
|
| 218 |
221 |
}
|
| 219 |
222 |
for (int e = 0; e < _next_arc; ++e) {
|
| 220 |
223 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 221 |
224 |
if (c < 0) {
|
| 222 |
225 |
_in_arc = e;
|
| 223 |
226 |
_next_arc = e + 1;
|
| 224 |
227 |
return true;
|
| 225 |
228 |
}
|
| 226 |
229 |
}
|
| 227 |
230 |
return false;
|
| 228 |
231 |
}
|
| 229 |
232 |
|
| 230 |
233 |
}; //class FirstEligiblePivotRule
|
| 231 |
234 |
|
| 232 |
235 |
|
| 233 |
|
/// \brief Implementation of the "Best Eligible" pivot rule for the
|
| 234 |
|
/// \ref NetworkSimplex "network simplex" algorithm.
|
| 235 |
|
///
|
| 236 |
|
/// This class implements the "Best Eligible" pivot rule
|
| 237 |
|
/// for the \ref NetworkSimplex "network simplex" algorithm.
|
| 238 |
|
///
|
| 239 |
|
/// For more information see \ref NetworkSimplex::run().
|
|
236 |
// Implementation of the Best Eligible pivot rule
|
| 240 |
237 |
class BestEligiblePivotRule
|
| 241 |
238 |
{
|
| 242 |
239 |
private:
|
| 243 |
240 |
|
| 244 |
241 |
// References to the NetworkSimplex class
|
| 245 |
242 |
const IntVector &_source;
|
| 246 |
243 |
const IntVector &_target;
|
| 247 |
|
const CostVector &_cost;
|
|
244 |
const ValueVector &_cost;
|
| 248 |
245 |
const IntVector &_state;
|
| 249 |
|
const CostVector &_pi;
|
|
246 |
const ValueVector &_pi;
|
| 250 |
247 |
int &_in_arc;
|
| 251 |
248 |
int _arc_num;
|
| 252 |
249 |
|
| 253 |
250 |
public:
|
| 254 |
251 |
|
| 255 |
|
/// Constructor
|
|
252 |
// Constructor
|
| 256 |
253 |
BestEligiblePivotRule(NetworkSimplex &ns) :
|
| 257 |
254 |
_source(ns._source), _target(ns._target),
|
| 258 |
255 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 259 |
256 |
_in_arc(ns.in_arc), _arc_num(ns._arc_num)
|
| 260 |
257 |
{}
|
| 261 |
258 |
|
| 262 |
|
/// Find next entering arc
|
|
259 |
// Find next entering arc
|
| 263 |
260 |
bool findEnteringArc() {
|
| 264 |
|
Cost c, min = 0;
|
|
261 |
Value c, min = 0;
|
| 265 |
262 |
for (int e = 0; e < _arc_num; ++e) {
|
| 266 |
263 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 267 |
264 |
if (c < min) {
|
| 268 |
265 |
min = c;
|
| 269 |
266 |
_in_arc = e;
|
| 270 |
267 |
}
|
| 271 |
268 |
}
|
| 272 |
269 |
return min < 0;
|
| 273 |
270 |
}
|
| 274 |
271 |
|
| 275 |
272 |
}; //class BestEligiblePivotRule
|
| 276 |
273 |
|
| 277 |
274 |
|
| 278 |
|
/// \brief Implementation of the "Block Search" pivot rule for the
|
| 279 |
|
/// \ref NetworkSimplex "network simplex" algorithm.
|
| 280 |
|
///
|
| 281 |
|
/// This class implements the "Block Search" pivot rule
|
| 282 |
|
/// for the \ref NetworkSimplex "network simplex" algorithm.
|
| 283 |
|
///
|
| 284 |
|
/// For more information see \ref NetworkSimplex::run().
|
|
275 |
// Implementation of the Block Search pivot rule
|
| 285 |
276 |
class BlockSearchPivotRule
|
| 286 |
277 |
{
|
| 287 |
278 |
private:
|
| 288 |
279 |
|
| 289 |
280 |
// References to the NetworkSimplex class
|
| 290 |
281 |
const IntVector &_source;
|
| 291 |
282 |
const IntVector &_target;
|
| 292 |
|
const CostVector &_cost;
|
|
283 |
const ValueVector &_cost;
|
| 293 |
284 |
const IntVector &_state;
|
| 294 |
|
const CostVector &_pi;
|
|
285 |
const ValueVector &_pi;
|
| 295 |
286 |
int &_in_arc;
|
| 296 |
287 |
int _arc_num;
|
| 297 |
288 |
|
| 298 |
289 |
// Pivot rule data
|
| 299 |
290 |
int _block_size;
|
| 300 |
291 |
int _next_arc;
|
| 301 |
292 |
|
| 302 |
293 |
public:
|
| 303 |
294 |
|
| 304 |
|
/// Constructor
|
|
295 |
// Constructor
|
| 305 |
296 |
BlockSearchPivotRule(NetworkSimplex &ns) :
|
| 306 |
297 |
_source(ns._source), _target(ns._target),
|
| 307 |
298 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 308 |
299 |
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
| 309 |
300 |
{
|
| 310 |
301 |
// The main parameters of the pivot rule
|
| 311 |
302 |
const double BLOCK_SIZE_FACTOR = 2.0;
|
| 312 |
303 |
const int MIN_BLOCK_SIZE = 10;
|
| 313 |
304 |
|
| 314 |
305 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
|
| 315 |
306 |
MIN_BLOCK_SIZE );
|
| 316 |
307 |
}
|
| 317 |
308 |
|
| 318 |
|
/// Find next entering arc
|
|
309 |
// Find next entering arc
|
| 319 |
310 |
bool findEnteringArc() {
|
| 320 |
|
Cost c, min = 0;
|
|
311 |
Value c, min = 0;
|
| 321 |
312 |
int cnt = _block_size;
|
| 322 |
313 |
int e, min_arc = _next_arc;
|
| 323 |
314 |
for (e = _next_arc; e < _arc_num; ++e) {
|
| 324 |
315 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 325 |
316 |
if (c < min) {
|
| 326 |
317 |
min = c;
|
| 327 |
318 |
min_arc = e;
|
| 328 |
319 |
}
|
| 329 |
320 |
if (--cnt == 0) {
|
| 330 |
321 |
if (min < 0) break;
|
| 331 |
322 |
cnt = _block_size;
|
| 332 |
323 |
}
|
| 333 |
324 |
}
|
| 334 |
325 |
if (min == 0 || cnt > 0) {
|
| 335 |
326 |
for (e = 0; e < _next_arc; ++e) {
|
| 336 |
327 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 337 |
328 |
if (c < min) {
|
| 338 |
329 |
min = c;
|
| 339 |
330 |
min_arc = e;
|
| 340 |
331 |
}
|
| 341 |
332 |
if (--cnt == 0) {
|
| 342 |
333 |
if (min < 0) break;
|
| 343 |
334 |
cnt = _block_size;
|
| 344 |
335 |
}
|
| 345 |
336 |
}
|
| 346 |
337 |
}
|
| 347 |
338 |
if (min >= 0) return false;
|
| 348 |
339 |
_in_arc = min_arc;
|
| 349 |
340 |
_next_arc = e;
|
| 350 |
341 |
return true;
|
| 351 |
342 |
}
|
| 352 |
343 |
|
| 353 |
344 |
}; //class BlockSearchPivotRule
|
| 354 |
345 |
|
| 355 |
346 |
|
| 356 |
|
/// \brief Implementation of the "Candidate List" pivot rule for the
|
| 357 |
|
/// \ref NetworkSimplex "network simplex" algorithm.
|
| 358 |
|
///
|
| 359 |
|
/// This class implements the "Candidate List" pivot rule
|
| 360 |
|
/// for the \ref NetworkSimplex "network simplex" algorithm.
|
| 361 |
|
///
|
| 362 |
|
/// For more information see \ref NetworkSimplex::run().
|
|
347 |
// Implementation of the Candidate List pivot rule
|
| 363 |
348 |
class CandidateListPivotRule
|
| 364 |
349 |
{
|
| 365 |
350 |
private:
|
| 366 |
351 |
|
| 367 |
352 |
// References to the NetworkSimplex class
|
| 368 |
353 |
const IntVector &_source;
|
| 369 |
354 |
const IntVector &_target;
|
| 370 |
|
const CostVector &_cost;
|
|
355 |
const ValueVector &_cost;
|
| 371 |
356 |
const IntVector &_state;
|
| 372 |
|
const CostVector &_pi;
|
|
357 |
const ValueVector &_pi;
|
| 373 |
358 |
int &_in_arc;
|
| 374 |
359 |
int _arc_num;
|
| 375 |
360 |
|
| 376 |
361 |
// Pivot rule data
|
| 377 |
362 |
IntVector _candidates;
|
| 378 |
363 |
int _list_length, _minor_limit;
|
| 379 |
364 |
int _curr_length, _minor_count;
|
| 380 |
365 |
int _next_arc;
|
| 381 |
366 |
|
| 382 |
367 |
public:
|
| 383 |
368 |
|
| 384 |
369 |
/// Constructor
|
| 385 |
370 |
CandidateListPivotRule(NetworkSimplex &ns) :
|
| 386 |
371 |
_source(ns._source), _target(ns._target),
|
| 387 |
372 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 388 |
373 |
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
| 389 |
374 |
{
|
| 390 |
375 |
// The main parameters of the pivot rule
|
| 391 |
376 |
const double LIST_LENGTH_FACTOR = 1.0;
|
| 392 |
377 |
const int MIN_LIST_LENGTH = 10;
|
| 393 |
378 |
const double MINOR_LIMIT_FACTOR = 0.1;
|
| 394 |
379 |
const int MIN_MINOR_LIMIT = 3;
|
| 395 |
380 |
|
| 396 |
381 |
_list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
|
| 397 |
382 |
MIN_LIST_LENGTH );
|
| 398 |
383 |
_minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
|
| 399 |
384 |
MIN_MINOR_LIMIT );
|
| 400 |
385 |
_curr_length = _minor_count = 0;
|
| 401 |
386 |
_candidates.resize(_list_length);
|
| 402 |
387 |
}
|
| 403 |
388 |
|
| 404 |
389 |
/// Find next entering arc
|
| 405 |
390 |
bool findEnteringArc() {
|
| 406 |
|
Cost min, c;
|
|
391 |
Value min, c;
|
| 407 |
392 |
int e, min_arc = _next_arc;
|
| 408 |
393 |
if (_curr_length > 0 && _minor_count < _minor_limit) {
|
| 409 |
394 |
// Minor iteration: select the best eligible arc from the
|
| 410 |
395 |
// current candidate list
|
| 411 |
396 |
++_minor_count;
|
| 412 |
397 |
min = 0;
|
| 413 |
398 |
for (int i = 0; i < _curr_length; ++i) {
|
| 414 |
399 |
e = _candidates[i];
|
| 415 |
400 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 416 |
401 |
if (c < min) {
|
| 417 |
402 |
min = c;
|
| 418 |
403 |
min_arc = e;
|
| 419 |
404 |
}
|
| 420 |
405 |
if (c >= 0) {
|
| 421 |
406 |
_candidates[i--] = _candidates[--_curr_length];
|
| 422 |
407 |
}
|
| 423 |
408 |
}
|
| 424 |
409 |
if (min < 0) {
|
| 425 |
410 |
_in_arc = min_arc;
|
| 426 |
411 |
return true;
|
| 427 |
412 |
}
|
| 428 |
413 |
}
|
| 429 |
414 |
|
| 430 |
415 |
// Major iteration: build a new candidate list
|
| ... |
... |
@@ -443,112 +428,106 @@
|
| 443 |
428 |
}
|
| 444 |
429 |
if (_curr_length < _list_length) {
|
| 445 |
430 |
for (e = 0; e < _next_arc; ++e) {
|
| 446 |
431 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 447 |
432 |
if (c < 0) {
|
| 448 |
433 |
_candidates[_curr_length++] = e;
|
| 449 |
434 |
if (c < min) {
|
| 450 |
435 |
min = c;
|
| 451 |
436 |
min_arc = e;
|
| 452 |
437 |
}
|
| 453 |
438 |
if (_curr_length == _list_length) break;
|
| 454 |
439 |
}
|
| 455 |
440 |
}
|
| 456 |
441 |
}
|
| 457 |
442 |
if (_curr_length == 0) return false;
|
| 458 |
443 |
_minor_count = 1;
|
| 459 |
444 |
_in_arc = min_arc;
|
| 460 |
445 |
_next_arc = e;
|
| 461 |
446 |
return true;
|
| 462 |
447 |
}
|
| 463 |
448 |
|
| 464 |
449 |
}; //class CandidateListPivotRule
|
| 465 |
450 |
|
| 466 |
451 |
|
| 467 |
|
/// \brief Implementation of the "Altering Candidate List" pivot rule
|
| 468 |
|
/// for the \ref NetworkSimplex "network simplex" algorithm.
|
| 469 |
|
///
|
| 470 |
|
/// This class implements the "Altering Candidate List" pivot rule
|
| 471 |
|
/// for the \ref NetworkSimplex "network simplex" algorithm.
|
| 472 |
|
///
|
| 473 |
|
/// For more information see \ref NetworkSimplex::run().
|
|
452 |
// Implementation of the Altering Candidate List pivot rule
|
| 474 |
453 |
class AlteringListPivotRule
|
| 475 |
454 |
{
|
| 476 |
455 |
private:
|
| 477 |
456 |
|
| 478 |
457 |
// References to the NetworkSimplex class
|
| 479 |
458 |
const IntVector &_source;
|
| 480 |
459 |
const IntVector &_target;
|
| 481 |
|
const CostVector &_cost;
|
|
460 |
const ValueVector &_cost;
|
| 482 |
461 |
const IntVector &_state;
|
| 483 |
|
const CostVector &_pi;
|
|
462 |
const ValueVector &_pi;
|
| 484 |
463 |
int &_in_arc;
|
| 485 |
464 |
int _arc_num;
|
| 486 |
465 |
|
| 487 |
466 |
// Pivot rule data
|
| 488 |
467 |
int _block_size, _head_length, _curr_length;
|
| 489 |
468 |
int _next_arc;
|
| 490 |
469 |
IntVector _candidates;
|
| 491 |
|
CostVector _cand_cost;
|
|
470 |
ValueVector _cand_cost;
|
| 492 |
471 |
|
| 493 |
472 |
// Functor class to compare arcs during sort of the candidate list
|
| 494 |
473 |
class SortFunc
|
| 495 |
474 |
{
|
| 496 |
475 |
private:
|
| 497 |
|
const CostVector &_map;
|
|
476 |
const ValueVector &_map;
|
| 498 |
477 |
public:
|
| 499 |
|
SortFunc(const CostVector &map) : _map(map) {}
|
|
478 |
SortFunc(const ValueVector &map) : _map(map) {}
|
| 500 |
479 |
bool operator()(int left, int right) {
|
| 501 |
480 |
return _map[left] > _map[right];
|
| 502 |
481 |
}
|
| 503 |
482 |
};
|
| 504 |
483 |
|
| 505 |
484 |
SortFunc _sort_func;
|
| 506 |
485 |
|
| 507 |
486 |
public:
|
| 508 |
487 |
|
| 509 |
|
/// Constructor
|
|
488 |
// Constructor
|
| 510 |
489 |
AlteringListPivotRule(NetworkSimplex &ns) :
|
| 511 |
490 |
_source(ns._source), _target(ns._target),
|
| 512 |
491 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 513 |
492 |
_in_arc(ns.in_arc), _arc_num(ns._arc_num),
|
| 514 |
493 |
_next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
|
| 515 |
494 |
{
|
| 516 |
495 |
// The main parameters of the pivot rule
|
| 517 |
496 |
const double BLOCK_SIZE_FACTOR = 1.5;
|
| 518 |
497 |
const int MIN_BLOCK_SIZE = 10;
|
| 519 |
498 |
const double HEAD_LENGTH_FACTOR = 0.1;
|
| 520 |
499 |
const int MIN_HEAD_LENGTH = 3;
|
| 521 |
500 |
|
| 522 |
501 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
|
| 523 |
502 |
MIN_BLOCK_SIZE );
|
| 524 |
503 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
|
| 525 |
504 |
MIN_HEAD_LENGTH );
|
| 526 |
505 |
_candidates.resize(_head_length + _block_size);
|
| 527 |
506 |
_curr_length = 0;
|
| 528 |
507 |
}
|
| 529 |
508 |
|
| 530 |
|
/// Find next entering arc
|
|
509 |
// Find next entering arc
|
| 531 |
510 |
bool findEnteringArc() {
|
| 532 |
511 |
// Check the current candidate list
|
| 533 |
512 |
int e;
|
| 534 |
513 |
for (int i = 0; i < _curr_length; ++i) {
|
| 535 |
514 |
e = _candidates[i];
|
| 536 |
515 |
_cand_cost[e] = _state[e] *
|
| 537 |
516 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 538 |
517 |
if (_cand_cost[e] >= 0) {
|
| 539 |
518 |
_candidates[i--] = _candidates[--_curr_length];
|
| 540 |
519 |
}
|
| 541 |
520 |
}
|
| 542 |
521 |
|
| 543 |
522 |
// Extend the list
|
| 544 |
523 |
int cnt = _block_size;
|
| 545 |
524 |
int last_arc = 0;
|
| 546 |
525 |
int limit = _head_length;
|
| 547 |
526 |
|
| 548 |
527 |
for (int e = _next_arc; e < _arc_num; ++e) {
|
| 549 |
528 |
_cand_cost[e] = _state[e] *
|
| 550 |
529 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 551 |
530 |
if (_cand_cost[e] < 0) {
|
| 552 |
531 |
_candidates[_curr_length++] = e;
|
| 553 |
532 |
last_arc = e;
|
| 554 |
533 |
}
|
| ... |
... |
@@ -571,383 +550,502 @@
|
| 571 |
550 |
limit = 0;
|
| 572 |
551 |
cnt = _block_size;
|
| 573 |
552 |
}
|
| 574 |
553 |
}
|
| 575 |
554 |
}
|
| 576 |
555 |
if (_curr_length == 0) return false;
|
| 577 |
556 |
_next_arc = last_arc + 1;
|
| 578 |
557 |
|
| 579 |
558 |
// Make heap of the candidate list (approximating a partial sort)
|
| 580 |
559 |
make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
|
| 581 |
560 |
_sort_func );
|
| 582 |
561 |
|
| 583 |
562 |
// Pop the first element of the heap
|
| 584 |
563 |
_in_arc = _candidates[0];
|
| 585 |
564 |
pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
|
| 586 |
565 |
_sort_func );
|
| 587 |
566 |
_curr_length = std::min(_head_length, _curr_length - 1);
|
| 588 |
567 |
return true;
|
| 589 |
568 |
}
|
| 590 |
569 |
|
| 591 |
570 |
}; //class AlteringListPivotRule
|
| 592 |
571 |
|
| 593 |
572 |
public:
|
| 594 |
573 |
|
| 595 |
|
/// \brief General constructor (with lower bounds).
|
|
574 |
/// \brief Constructor.
|
| 596 |
575 |
///
|
| 597 |
|
/// General constructor (with lower bounds).
|
|
576 |
/// Constructor.
|
| 598 |
577 |
///
|
| 599 |
578 |
/// \param graph The digraph the algorithm runs on.
|
| 600 |
|
/// \param lower The lower bounds of the arcs.
|
| 601 |
|
/// \param capacity The capacities (upper bounds) of the arcs.
|
| 602 |
|
/// \param cost The cost (length) values of the arcs.
|
| 603 |
|
/// \param supply The supply values of the nodes (signed).
|
| 604 |
|
NetworkSimplex( const Digraph &graph,
|
| 605 |
|
const LowerMap &lower,
|
| 606 |
|
const CapacityMap &capacity,
|
| 607 |
|
const CostMap &cost,
|
| 608 |
|
const SupplyMap &supply ) :
|
| 609 |
|
_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
|
| 610 |
|
_orig_cost(cost), _orig_supply(&supply),
|
|
579 |
NetworkSimplex(const GR& graph) :
|
|
580 |
_graph(graph),
|
|
581 |
_plower(NULL), _pupper(NULL), _pcost(NULL),
|
|
582 |
_psupply(NULL), _pstsup(false),
|
| 611 |
583 |
_flow_map(NULL), _potential_map(NULL),
|
| 612 |
584 |
_local_flow(false), _local_potential(false),
|
| 613 |
585 |
_node_id(graph)
|
| 614 |
|
{}
|
| 615 |
|
|
| 616 |
|
/// \brief General constructor (without lower bounds).
|
| 617 |
|
///
|
| 618 |
|
/// General constructor (without lower bounds).
|
| 619 |
|
///
|
| 620 |
|
/// \param graph The digraph the algorithm runs on.
|
| 621 |
|
/// \param capacity The capacities (upper bounds) of the arcs.
|
| 622 |
|
/// \param cost The cost (length) values of the arcs.
|
| 623 |
|
/// \param supply The supply values of the nodes (signed).
|
| 624 |
|
NetworkSimplex( const Digraph &graph,
|
| 625 |
|
const CapacityMap &capacity,
|
| 626 |
|
const CostMap &cost,
|
| 627 |
|
const SupplyMap &supply ) :
|
| 628 |
|
_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
|
| 629 |
|
_orig_cost(cost), _orig_supply(&supply),
|
| 630 |
|
_flow_map(NULL), _potential_map(NULL),
|
| 631 |
|
_local_flow(false), _local_potential(false),
|
| 632 |
|
_node_id(graph)
|
| 633 |
|
{}
|
| 634 |
|
|
| 635 |
|
/// \brief Simple constructor (with lower bounds).
|
| 636 |
|
///
|
| 637 |
|
/// Simple constructor (with lower bounds).
|
| 638 |
|
///
|
| 639 |
|
/// \param graph The digraph the algorithm runs on.
|
| 640 |
|
/// \param lower The lower bounds of the arcs.
|
| 641 |
|
/// \param capacity The capacities (upper bounds) of the arcs.
|
| 642 |
|
/// \param cost The cost (length) values of the arcs.
|
| 643 |
|
/// \param s The source node.
|
| 644 |
|
/// \param t The target node.
|
| 645 |
|
/// \param flow_value The required amount of flow from node \c s
|
| 646 |
|
/// to node \c t (i.e. the supply of \c s and the demand of \c t).
|
| 647 |
|
NetworkSimplex( const Digraph &graph,
|
| 648 |
|
const LowerMap &lower,
|
| 649 |
|
const CapacityMap &capacity,
|
| 650 |
|
const CostMap &cost,
|
| 651 |
|
Node s, Node t,
|
| 652 |
|
Capacity flow_value ) :
|
| 653 |
|
_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
|
| 654 |
|
_orig_cost(cost), _orig_supply(NULL),
|
| 655 |
|
_orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
|
| 656 |
|
_flow_map(NULL), _potential_map(NULL),
|
| 657 |
|
_local_flow(false), _local_potential(false),
|
| 658 |
|
_node_id(graph)
|
| 659 |
|
{}
|
| 660 |
|
|
| 661 |
|
/// \brief Simple constructor (without lower bounds).
|
| 662 |
|
///
|
| 663 |
|
/// Simple constructor (without lower bounds).
|
| 664 |
|
///
|
| 665 |
|
/// \param graph The digraph the algorithm runs on.
|
| 666 |
|
/// \param capacity The capacities (upper bounds) of the arcs.
|
| 667 |
|
/// \param cost The cost (length) values of the arcs.
|
| 668 |
|
/// \param s The source node.
|
| 669 |
|
/// \param t The target node.
|
| 670 |
|
/// \param flow_value The required amount of flow from node \c s
|
| 671 |
|
/// to node \c t (i.e. the supply of \c s and the demand of \c t).
|
| 672 |
|
NetworkSimplex( const Digraph &graph,
|
| 673 |
|
const CapacityMap &capacity,
|
| 674 |
|
const CostMap &cost,
|
| 675 |
|
Node s, Node t,
|
| 676 |
|
Capacity flow_value ) :
|
| 677 |
|
_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
|
| 678 |
|
_orig_cost(cost), _orig_supply(NULL),
|
| 679 |
|
_orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
|
| 680 |
|
_flow_map(NULL), _potential_map(NULL),
|
| 681 |
|
_local_flow(false), _local_potential(false),
|
| 682 |
|
_node_id(graph)
|
| 683 |
|
{}
|
|
586 |
{
|
|
587 |
LEMON_ASSERT(std::numeric_limits<Value>::is_integer &&
|
|
588 |
std::numeric_limits<Value>::is_signed,
|
|
589 |
"The value type of NetworkSimplex must be a signed integer");
|
|
590 |
}
|
| 684 |
591 |
|
| 685 |
592 |
/// Destructor.
|
| 686 |
593 |
~NetworkSimplex() {
|
| 687 |
594 |
if (_local_flow) delete _flow_map;
|
| 688 |
595 |
if (_local_potential) delete _potential_map;
|
| 689 |
596 |
}
|
| 690 |
597 |
|
|
598 |
/// \brief Set the lower bounds on the arcs.
|
|
599 |
///
|
|
600 |
/// This function sets the lower bounds on the arcs.
|
|
601 |
/// If neither this function nor \ref boundMaps() is used before
|
|
602 |
/// calling \ref run(), the lower bounds will be set to zero
|
|
603 |
/// on all arcs.
|
|
604 |
///
|
|
605 |
/// \param map An arc map storing the lower bounds.
|
|
606 |
/// Its \c Value type must be convertible to the \c Value type
|
|
607 |
/// of the algorithm.
|
|
608 |
///
|
|
609 |
/// \return <tt>(*this)</tt>
|
|
610 |
template <typename LOWER>
|
|
611 |
NetworkSimplex& lowerMap(const LOWER& map) {
|
|
612 |
delete _plower;
|
|
613 |
_plower = new ValueArcMap(_graph);
|
|
614 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
|
615 |
(*_plower)[a] = map[a];
|
|
616 |
}
|
|
617 |
return *this;
|
|
618 |
}
|
|
619 |
|
|
620 |
/// \brief Set the upper bounds (capacities) on the arcs.
|
|
621 |
///
|
|
622 |
/// This function sets the upper bounds (capacities) on the arcs.
|
|
623 |
/// If none of the functions \ref upperMap(), \ref capacityMap()
|
|
624 |
/// and \ref boundMaps() is used before calling \ref run(),
|
|
625 |
/// the upper bounds (capacities) will be set to
|
|
626 |
/// \c std::numeric_limits<Value>::max() on all arcs.
|
|
627 |
///
|
|
628 |
/// \param map An arc map storing the upper bounds.
|
|
629 |
/// Its \c Value type must be convertible to the \c Value type
|
|
630 |
/// of the algorithm.
|
|
631 |
///
|
|
632 |
/// \return <tt>(*this)</tt>
|
|
633 |
template<typename UPPER>
|
|
634 |
NetworkSimplex& upperMap(const UPPER& map) {
|
|
635 |
delete _pupper;
|
|
636 |
_pupper = new ValueArcMap(_graph);
|
|
637 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
|
638 |
(*_pupper)[a] = map[a];
|
|
639 |
}
|
|
640 |
return *this;
|
|
641 |
}
|
|
642 |
|
|
643 |
/// \brief Set the upper bounds (capacities) on the arcs.
|
|
644 |
///
|
|
645 |
/// This function sets the upper bounds (capacities) on the arcs.
|
|
646 |
/// It is just an alias for \ref upperMap().
|
|
647 |
///
|
|
648 |
/// \return <tt>(*this)</tt>
|
|
649 |
template<typename CAP>
|
|
650 |
NetworkSimplex& capacityMap(const CAP& map) {
|
|
651 |
return upperMap(map);
|
|
652 |
}
|
|
653 |
|
|
654 |
/// \brief Set the lower and upper bounds on the arcs.
|
|
655 |
///
|
|
656 |
/// This function sets the lower and upper bounds on the arcs.
|
|
657 |
/// If neither this function nor \ref lowerMap() is used before
|
|
658 |
/// calling \ref run(), the lower bounds will be set to zero
|
|
659 |
/// on all arcs.
|
|
660 |
/// If none of the functions \ref upperMap(), \ref capacityMap()
|
|
661 |
/// and \ref boundMaps() is used before calling \ref run(),
|
|
662 |
/// the upper bounds (capacities) will be set to
|
|
663 |
/// \c std::numeric_limits<Value>::max() on all arcs.
|
|
664 |
///
|
|
665 |
/// \param lower An arc map storing the lower bounds.
|
|
666 |
/// \param upper An arc map storing the upper bounds.
|
|
667 |
///
|
|
668 |
/// The \c Value type of the maps must be convertible to the
|
|
669 |
/// \c Value type of the algorithm.
|
|
670 |
///
|
|
671 |
/// \note This function is just a shortcut of calling \ref lowerMap()
|
|
672 |
/// and \ref upperMap() separately.
|
|
673 |
///
|
|
674 |
/// \return <tt>(*this)</tt>
|
|
675 |
template <typename LOWER, typename UPPER>
|
|
676 |
NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
|
|
677 |
return lowerMap(lower).upperMap(upper);
|
|
678 |
}
|
|
679 |
|
|
680 |
/// \brief Set the costs of the arcs.
|
|
681 |
///
|
|
682 |
/// This function sets the costs of the arcs.
|
|
683 |
/// If it is not used before calling \ref run(), the costs
|
|
684 |
/// will be set to \c 1 on all arcs.
|
|
685 |
///
|
|
686 |
/// \param map An arc map storing the costs.
|
|
687 |
/// Its \c Value type must be convertible to the \c Value type
|
|
688 |
/// of the algorithm.
|
|
689 |
///
|
|
690 |
/// \return <tt>(*this)</tt>
|
|
691 |
template<typename COST>
|
|
692 |
NetworkSimplex& costMap(const COST& map) {
|
|
693 |
delete _pcost;
|
|
694 |
_pcost = new ValueArcMap(_graph);
|
|
695 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
|
696 |
(*_pcost)[a] = map[a];
|
|
697 |
}
|
|
698 |
return *this;
|
|
699 |
}
|
|
700 |
|
|
701 |
/// \brief Set the supply values of the nodes.
|
|
702 |
///
|
|
703 |
/// This function sets the supply values of the nodes.
|
|
704 |
/// If neither this function nor \ref stSupply() is used before
|
|
705 |
/// calling \ref run(), the supply of each node will be set to zero.
|
|
706 |
/// (It makes sense only if non-zero lower bounds are given.)
|
|
707 |
///
|
|
708 |
/// \param map A node map storing the supply values.
|
|
709 |
/// Its \c Value type must be convertible to the \c Value type
|
|
710 |
/// of the algorithm.
|
|
711 |
///
|
|
712 |
/// \return <tt>(*this)</tt>
|
|
713 |
template<typename SUP>
|
|
714 |
NetworkSimplex& supplyMap(const SUP& map) {
|
|
715 |
delete _psupply;
|
|
716 |
_pstsup = false;
|
|
717 |
_psupply = new ValueNodeMap(_graph);
|
|
718 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
719 |
(*_psupply)[n] = map[n];
|
|
720 |
}
|
|
721 |
return *this;
|
|
722 |
}
|
|
723 |
|
|
724 |
/// \brief Set single source and target nodes and a supply value.
|
|
725 |
///
|
|
726 |
/// This function sets a single source node and a single target node
|
|
727 |
/// and the required flow value.
|
|
728 |
/// If neither this function nor \ref supplyMap() is used before
|
|
729 |
/// calling \ref run(), the supply of each node will be set to zero.
|
|
730 |
/// (It makes sense only if non-zero lower bounds are given.)
|
|
731 |
///
|
|
732 |
/// \param s The source node.
|
|
733 |
/// \param t The target node.
|
|
734 |
/// \param k The required amount of flow from node \c s to node \c t
|
|
735 |
/// (i.e. the supply of \c s and the demand of \c t).
|
|
736 |
///
|
|
737 |
/// \return <tt>(*this)</tt>
|
|
738 |
NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
|
|
739 |
delete _psupply;
|
|
740 |
_psupply = NULL;
|
|
741 |
_pstsup = true;
|
|
742 |
_psource = s;
|
|
743 |
_ptarget = t;
|
|
744 |
_pstflow = k;
|
|
745 |
return *this;
|
|
746 |
}
|
|
747 |
|
| 691 |
748 |
/// \brief Set the flow map.
|
| 692 |
749 |
///
|
| 693 |
750 |
/// This function sets the flow map.
|
|
751 |
/// If it is not used before calling \ref run(), an instance will
|
|
752 |
/// be allocated automatically. The destructor deallocates this
|
|
753 |
/// automatically allocated map, of course.
|
| 694 |
754 |
///
|
| 695 |
755 |
/// \return <tt>(*this)</tt>
|
| 696 |
756 |
NetworkSimplex& flowMap(FlowMap &map) {
|
| 697 |
757 |
if (_local_flow) {
|
| 698 |
758 |
delete _flow_map;
|
| 699 |
759 |
_local_flow = false;
|
| 700 |
760 |
}
|
| 701 |
761 |
_flow_map = ↦
|
| 702 |
762 |
return *this;
|
| 703 |
763 |
}
|
| 704 |
764 |
|
| 705 |
765 |
/// \brief Set the potential map.
|
| 706 |
766 |
///
|
| 707 |
|
/// This function sets the potential map.
|
|
767 |
/// This function sets the potential map, which is used for storing
|
|
768 |
/// the dual solution.
|
|
769 |
/// If it is not used before calling \ref run(), an instance will
|
|
770 |
/// be allocated automatically. The destructor deallocates this
|
|
771 |
/// automatically allocated map, of course.
|
| 708 |
772 |
///
|
| 709 |
773 |
/// \return <tt>(*this)</tt>
|
| 710 |
774 |
NetworkSimplex& potentialMap(PotentialMap &map) {
|
| 711 |
775 |
if (_local_potential) {
|
| 712 |
776 |
delete _potential_map;
|
| 713 |
777 |
_local_potential = false;
|
| 714 |
778 |
}
|
| 715 |
779 |
_potential_map = ↦
|
| 716 |
780 |
return *this;
|
| 717 |
781 |
}
|
| 718 |
782 |
|
| 719 |
|
/// \name Execution control
|
| 720 |
|
/// The algorithm can be executed using the
|
| 721 |
|
/// \ref lemon::NetworkSimplex::run() "run()" function.
|
|
783 |
/// \name Execution Control
|
|
784 |
/// The algorithm can be executed using \ref run().
|
|
785 |
|
| 722 |
786 |
/// @{
|
| 723 |
787 |
|
| 724 |
788 |
/// \brief Run the algorithm.
|
| 725 |
789 |
///
|
| 726 |
790 |
/// This function runs the algorithm.
|
|
791 |
/// The paramters can be specified using \ref lowerMap(),
|
|
792 |
/// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
|
|
793 |
/// \ref costMap(), \ref supplyMap() and \ref stSupply()
|
|
794 |
/// functions. For example,
|
|
795 |
/// \code
|
|
796 |
/// NetworkSimplex<ListDigraph> ns(graph);
|
|
797 |
/// ns.boundMaps(lower, upper).costMap(cost)
|
|
798 |
/// .supplyMap(sup).run();
|
|
799 |
/// \endcode
|
| 727 |
800 |
///
|
| 728 |
|
/// \param pivot_rule The pivot rule that is used during the
|
| 729 |
|
/// algorithm.
|
| 730 |
|
///
|
| 731 |
|
/// The available pivot rules:
|
| 732 |
|
///
|
| 733 |
|
/// - FIRST_ELIGIBLE_PIVOT The next eligible arc is selected in
|
| 734 |
|
/// a wraparound fashion in every iteration
|
| 735 |
|
/// (\ref FirstEligiblePivotRule).
|
| 736 |
|
///
|
| 737 |
|
/// - BEST_ELIGIBLE_PIVOT The best eligible arc is selected in
|
| 738 |
|
/// every iteration (\ref BestEligiblePivotRule).
|
| 739 |
|
///
|
| 740 |
|
/// - BLOCK_SEARCH_PIVOT A specified number of arcs are examined in
|
| 741 |
|
/// every iteration in a wraparound fashion and the best eligible
|
| 742 |
|
/// arc is selected from this block (\ref BlockSearchPivotRule).
|
| 743 |
|
///
|
| 744 |
|
/// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
|
| 745 |
|
/// built from eligible arcs in a wraparound fashion and in the
|
| 746 |
|
/// following minor iterations the best eligible arc is selected
|
| 747 |
|
/// from this list (\ref CandidateListPivotRule).
|
| 748 |
|
///
|
| 749 |
|
/// - ALTERING_LIST_PIVOT It is a modified version of the
|
| 750 |
|
/// "Candidate List" pivot rule. It keeps only the several best
|
| 751 |
|
/// eligible arcs from the former candidate list and extends this
|
| 752 |
|
/// list in every iteration (\ref AlteringListPivotRule).
|
| 753 |
|
///
|
| 754 |
|
/// According to our comprehensive benchmark tests the "Block Search"
|
| 755 |
|
/// pivot rule proved to be the fastest and the most robust on
|
| 756 |
|
/// various test inputs. Thus it is the default option.
|
|
801 |
/// \param pivot_rule The pivot rule that will be used during the
|
|
802 |
/// algorithm. For more information see \ref PivotRule.
|
| 757 |
803 |
///
|
| 758 |
804 |
/// \return \c true if a feasible flow can be found.
|
| 759 |
|
bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
|
|
805 |
bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
|
| 760 |
806 |
return init() && start(pivot_rule);
|
| 761 |
807 |
}
|
| 762 |
808 |
|
| 763 |
809 |
/// @}
|
| 764 |
810 |
|
| 765 |
811 |
/// \name Query Functions
|
| 766 |
812 |
/// The results of the algorithm can be obtained using these
|
| 767 |
813 |
/// functions.\n
|
| 768 |
|
/// \ref lemon::NetworkSimplex::run() "run()" must be called before
|
| 769 |
|
/// using them.
|
|
814 |
/// The \ref run() function must be called before using them.
|
|
815 |
|
| 770 |
816 |
/// @{
|
| 771 |
817 |
|
|
818 |
/// \brief Return the total cost of the found flow.
|
|
819 |
///
|
|
820 |
/// This function returns the total cost of the found flow.
|
|
821 |
/// The complexity of the function is \f$ O(e) \f$.
|
|
822 |
///
|
|
823 |
/// \note The return type of the function can be specified as a
|
|
824 |
/// template parameter. For example,
|
|
825 |
/// \code
|
|
826 |
/// ns.totalCost<double>();
|
|
827 |
/// \endcode
|
|
828 |
/// It is useful if the total cost cannot be stored in the \c Value
|
|
829 |
/// type of the algorithm, which is the default return type of the
|
|
830 |
/// function.
|
|
831 |
///
|
|
832 |
/// \pre \ref run() must be called before using this function.
|
|
833 |
template <typename Num>
|
|
834 |
Num totalCost() const {
|
|
835 |
Num c = 0;
|
|
836 |
if (_pcost) {
|
|
837 |
for (ArcIt e(_graph); e != INVALID; ++e)
|
|
838 |
c += (*_flow_map)[e] * (*_pcost)[e];
|
|
839 |
} else {
|
|
840 |
for (ArcIt e(_graph); e != INVALID; ++e)
|
|
841 |
c += (*_flow_map)[e];
|
|
842 |
}
|
|
843 |
return c;
|
|
844 |
}
|
|
845 |
|
|
846 |
#ifndef DOXYGEN
|
|
847 |
Value totalCost() const {
|
|
848 |
return totalCost<Value>();
|
|
849 |
}
|
|
850 |
#endif
|
|
851 |
|
|
852 |
/// \brief Return the flow on the given arc.
|
|
853 |
///
|
|
854 |
/// This function returns the flow on the given arc.
|
|
855 |
///
|
|
856 |
/// \pre \ref run() must be called before using this function.
|
|
857 |
Value flow(const Arc& a) const {
|
|
858 |
return (*_flow_map)[a];
|
|
859 |
}
|
|
860 |
|
| 772 |
861 |
/// \brief Return a const reference to the flow map.
|
| 773 |
862 |
///
|
| 774 |
863 |
/// This function returns a const reference to an arc map storing
|
| 775 |
864 |
/// the found flow.
|
| 776 |
865 |
///
|
| 777 |
866 |
/// \pre \ref run() must be called before using this function.
|
| 778 |
867 |
const FlowMap& flowMap() const {
|
| 779 |
868 |
return *_flow_map;
|
| 780 |
869 |
}
|
| 781 |
870 |
|
|
871 |
/// \brief Return the potential (dual value) of the given node.
|
|
872 |
///
|
|
873 |
/// This function returns the potential (dual value) of the
|
|
874 |
/// given node.
|
|
875 |
///
|
|
876 |
/// \pre \ref run() must be called before using this function.
|
|
877 |
Value potential(const Node& n) const {
|
|
878 |
return (*_potential_map)[n];
|
|
879 |
}
|
|
880 |
|
| 782 |
881 |
/// \brief Return a const reference to the potential map
|
| 783 |
882 |
/// (the dual solution).
|
| 784 |
883 |
///
|
| 785 |
884 |
/// This function returns a const reference to a node map storing
|
| 786 |
|
/// the found potentials (the dual solution).
|
|
885 |
/// the found potentials, which form the dual solution of the
|
|
886 |
/// \ref min_cost_flow "minimum cost flow" problem.
|
| 787 |
887 |
///
|
| 788 |
888 |
/// \pre \ref run() must be called before using this function.
|
| 789 |
889 |
const PotentialMap& potentialMap() const {
|
| 790 |
890 |
return *_potential_map;
|
| 791 |
891 |
}
|
| 792 |
892 |
|
| 793 |
|
/// \brief Return the flow on the given arc.
|
| 794 |
|
///
|
| 795 |
|
/// This function returns the flow on the given arc.
|
| 796 |
|
///
|
| 797 |
|
/// \pre \ref run() must be called before using this function.
|
| 798 |
|
Capacity flow(const Arc& arc) const {
|
| 799 |
|
return (*_flow_map)[arc];
|
| 800 |
|
}
|
| 801 |
|
|
| 802 |
|
/// \brief Return the potential of the given node.
|
| 803 |
|
///
|
| 804 |
|
/// This function returns the potential of the given node.
|
| 805 |
|
///
|
| 806 |
|
/// \pre \ref run() must be called before using this function.
|
| 807 |
|
Cost potential(const Node& node) const {
|
| 808 |
|
return (*_potential_map)[node];
|
| 809 |
|
}
|
| 810 |
|
|
| 811 |
|
/// \brief Return the total cost of the found flow.
|
| 812 |
|
///
|
| 813 |
|
/// This function returns the total cost of the found flow.
|
| 814 |
|
/// The complexity of the function is \f$ O(e) \f$.
|
| 815 |
|
///
|
| 816 |
|
/// \pre \ref run() must be called before using this function.
|
| 817 |
|
Cost totalCost() const {
|
| 818 |
|
Cost c = 0;
|
| 819 |
|
for (ArcIt e(_graph); e != INVALID; ++e)
|
| 820 |
|
c += (*_flow_map)[e] * _orig_cost[e];
|
| 821 |
|
return c;
|
| 822 |
|
}
|
| 823 |
|
|
| 824 |
893 |
/// @}
|
| 825 |
894 |
|
| 826 |
895 |
private:
|
| 827 |
896 |
|
| 828 |
897 |
// Initialize internal data structures
|
| 829 |
898 |
bool init() {
|
| 830 |
899 |
// Initialize result maps
|
| 831 |
900 |
if (!_flow_map) {
|
| 832 |
901 |
_flow_map = new FlowMap(_graph);
|
| 833 |
902 |
_local_flow = true;
|
| 834 |
903 |
}
|
| 835 |
904 |
if (!_potential_map) {
|
| 836 |
905 |
_potential_map = new PotentialMap(_graph);
|
| 837 |
906 |
_local_potential = true;
|
| 838 |
907 |
}
|
| 839 |
908 |
|
| 840 |
909 |
// Initialize vectors
|
| 841 |
910 |
_node_num = countNodes(_graph);
|
| 842 |
911 |
_arc_num = countArcs(_graph);
|
| 843 |
912 |
int all_node_num = _node_num + 1;
|
| 844 |
913 |
int all_arc_num = _arc_num + _node_num;
|
|
914 |
if (_node_num == 0) return false;
|
| 845 |
915 |
|
| 846 |
916 |
_arc_ref.resize(_arc_num);
|
| 847 |
917 |
_source.resize(all_arc_num);
|
| 848 |
918 |
_target.resize(all_arc_num);
|
| 849 |
919 |
|
| 850 |
920 |
_cap.resize(all_arc_num);
|
| 851 |
921 |
_cost.resize(all_arc_num);
|
| 852 |
922 |
_supply.resize(all_node_num);
|
| 853 |
923 |
_flow.resize(all_arc_num, 0);
|
| 854 |
924 |
_pi.resize(all_node_num, 0);
|
| 855 |
925 |
|
| 856 |
926 |
_parent.resize(all_node_num);
|
| 857 |
927 |
_pred.resize(all_node_num);
|
| 858 |
928 |
_forward.resize(all_node_num);
|
| 859 |
929 |
_thread.resize(all_node_num);
|
| 860 |
930 |
_rev_thread.resize(all_node_num);
|
| 861 |
931 |
_succ_num.resize(all_node_num);
|
| 862 |
932 |
_last_succ.resize(all_node_num);
|
| 863 |
933 |
_state.resize(all_arc_num, STATE_LOWER);
|
| 864 |
934 |
|
| 865 |
935 |
// Initialize node related data
|
| 866 |
936 |
bool valid_supply = true;
|
| 867 |
|
if (_orig_supply) {
|
| 868 |
|
Supply sum = 0;
|
|
937 |
if (!_pstsup && !_psupply) {
|
|
938 |
_pstsup = true;
|
|
939 |
_psource = _ptarget = NodeIt(_graph);
|
|
940 |
_pstflow = 0;
|
|
941 |
}
|
|
942 |
if (_psupply) {
|
|
943 |
Value sum = 0;
|
| 869 |
944 |
int i = 0;
|
| 870 |
945 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
|
| 871 |
946 |
_node_id[n] = i;
|
| 872 |
|
_supply[i] = (*_orig_supply)[n];
|
|
947 |
_supply[i] = (*_psupply)[n];
|
| 873 |
948 |
sum += _supply[i];
|
| 874 |
949 |
}
|
| 875 |
950 |
valid_supply = (sum == 0);
|
| 876 |
951 |
} else {
|
| 877 |
952 |
int i = 0;
|
| 878 |
953 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
|
| 879 |
954 |
_node_id[n] = i;
|
| 880 |
955 |
_supply[i] = 0;
|
| 881 |
956 |
}
|
| 882 |
|
_supply[_node_id[_orig_source]] = _orig_flow_value;
|
| 883 |
|
_supply[_node_id[_orig_target]] = -_orig_flow_value;
|
|
957 |
_supply[_node_id[_psource]] = _pstflow;
|
|
958 |
_supply[_node_id[_ptarget]] = -_pstflow;
|
| 884 |
959 |
}
|
| 885 |
960 |
if (!valid_supply) return false;
|
| 886 |
961 |
|
| 887 |
962 |
// Set data for the artificial root node
|
| 888 |
963 |
_root = _node_num;
|
| 889 |
964 |
_parent[_root] = -1;
|
| 890 |
965 |
_pred[_root] = -1;
|
| 891 |
966 |
_thread[_root] = 0;
|
| 892 |
967 |
_rev_thread[0] = _root;
|
| 893 |
968 |
_succ_num[_root] = all_node_num;
|
| 894 |
969 |
_last_succ[_root] = _root - 1;
|
| 895 |
970 |
_supply[_root] = 0;
|
| 896 |
971 |
_pi[_root] = 0;
|
| 897 |
972 |
|
| 898 |
973 |
// Store the arcs in a mixed order
|
| 899 |
974 |
int k = std::max(int(sqrt(_arc_num)), 10);
|
| 900 |
975 |
int i = 0;
|
| 901 |
976 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
| 902 |
977 |
_arc_ref[i] = e;
|
| 903 |
978 |
if ((i += k) >= _arc_num) i = (i % k) + 1;
|
| 904 |
979 |
}
|
| 905 |
980 |
|
| 906 |
981 |
// Initialize arc maps
|
|
982 |
if (_pupper && _pcost) {
|
| 907 |
983 |
for (int i = 0; i != _arc_num; ++i) {
|
| 908 |
984 |
Arc e = _arc_ref[i];
|
| 909 |
985 |
_source[i] = _node_id[_graph.source(e)];
|
| 910 |
986 |
_target[i] = _node_id[_graph.target(e)];
|
| 911 |
|
_cost[i] = _orig_cost[e];
|
| 912 |
|
_cap[i] = _orig_cap[e];
|
|
987 |
_cap[i] = (*_pupper)[e];
|
|
988 |
_cost[i] = (*_pcost)[e];
|
|
989 |
}
|
|
990 |
} else {
|
|
991 |
for (int i = 0; i != _arc_num; ++i) {
|
|
992 |
Arc e = _arc_ref[i];
|
|
993 |
_source[i] = _node_id[_graph.source(e)];
|
|
994 |
_target[i] = _node_id[_graph.target(e)];
|
|
995 |
}
|
|
996 |
if (_pupper) {
|
|
997 |
for (int i = 0; i != _arc_num; ++i)
|
|
998 |
_cap[i] = (*_pupper)[_arc_ref[i]];
|
|
999 |
} else {
|
|
1000 |
Value val = std::numeric_limits<Value>::max();
|
|
1001 |
for (int i = 0; i != _arc_num; ++i)
|
|
1002 |
_cap[i] = val;
|
|
1003 |
}
|
|
1004 |
if (_pcost) {
|
|
1005 |
for (int i = 0; i != _arc_num; ++i)
|
|
1006 |
_cost[i] = (*_pcost)[_arc_ref[i]];
|
|
1007 |
} else {
|
|
1008 |
for (int i = 0; i != _arc_num; ++i)
|
|
1009 |
_cost[i] = 1;
|
|
1010 |
}
|
| 913 |
1011 |
}
|
| 914 |
1012 |
|
| 915 |
1013 |
// Remove non-zero lower bounds
|
| 916 |
|
if (_orig_lower) {
|
|
1014 |
if (_plower) {
|
| 917 |
1015 |
for (int i = 0; i != _arc_num; ++i) {
|
| 918 |
|
Capacity c = (*_orig_lower)[_arc_ref[i]];
|
|
1016 |
Value c = (*_plower)[_arc_ref[i]];
|
| 919 |
1017 |
if (c != 0) {
|
| 920 |
1018 |
_cap[i] -= c;
|
| 921 |
1019 |
_supply[_source[i]] -= c;
|
| 922 |
1020 |
_supply[_target[i]] += c;
|
| 923 |
1021 |
}
|
| 924 |
1022 |
}
|
| 925 |
1023 |
}
|
| 926 |
1024 |
|
| 927 |
1025 |
// Add artificial arcs and initialize the spanning tree data structure
|
| 928 |
|
Cost max_cost = std::numeric_limits<Cost>::max() / 4;
|
| 929 |
|
Capacity max_cap = std::numeric_limits<Capacity>::max();
|
|
1026 |
Value max_cap = std::numeric_limits<Value>::max();
|
|
1027 |
Value max_cost = std::numeric_limits<Value>::max() / 4;
|
| 930 |
1028 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
| 931 |
1029 |
_thread[u] = u + 1;
|
| 932 |
1030 |
_rev_thread[u + 1] = u;
|
| 933 |
1031 |
_succ_num[u] = 1;
|
| 934 |
1032 |
_last_succ[u] = u;
|
| 935 |
1033 |
_parent[u] = _root;
|
| 936 |
1034 |
_pred[u] = e;
|
| 937 |
1035 |
if (_supply[u] >= 0) {
|
| 938 |
1036 |
_flow[e] = _supply[u];
|
| 939 |
1037 |
_forward[u] = true;
|
| 940 |
1038 |
_pi[u] = -max_cost;
|
| 941 |
1039 |
} else {
|
| 942 |
1040 |
_flow[e] = -_supply[u];
|
| 943 |
1041 |
_forward[u] = false;
|
| 944 |
1042 |
_pi[u] = max_cost;
|
| 945 |
1043 |
}
|
| 946 |
1044 |
_cost[e] = max_cost;
|
| 947 |
1045 |
_cap[e] = max_cap;
|
| 948 |
1046 |
_state[e] = STATE_TREE;
|
| 949 |
1047 |
}
|
| 950 |
1048 |
|
| 951 |
1049 |
return true;
|
| 952 |
1050 |
}
|
| 953 |
1051 |
|
| ... |
... |
@@ -958,87 +1056,87 @@
|
| 958 |
1056 |
while (u != v) {
|
| 959 |
1057 |
if (_succ_num[u] < _succ_num[v]) {
|
| 960 |
1058 |
u = _parent[u];
|
| 961 |
1059 |
} else {
|
| 962 |
1060 |
v = _parent[v];
|
| 963 |
1061 |
}
|
| 964 |
1062 |
}
|
| 965 |
1063 |
join = u;
|
| 966 |
1064 |
}
|
| 967 |
1065 |
|
| 968 |
1066 |
// Find the leaving arc of the cycle and returns true if the
|
| 969 |
1067 |
// leaving arc is not the same as the entering arc
|
| 970 |
1068 |
bool findLeavingArc() {
|
| 971 |
1069 |
// Initialize first and second nodes according to the direction
|
| 972 |
1070 |
// of the cycle
|
| 973 |
1071 |
if (_state[in_arc] == STATE_LOWER) {
|
| 974 |
1072 |
first = _source[in_arc];
|
| 975 |
1073 |
second = _target[in_arc];
|
| 976 |
1074 |
} else {
|
| 977 |
1075 |
first = _target[in_arc];
|
| 978 |
1076 |
second = _source[in_arc];
|
| 979 |
1077 |
}
|
| 980 |
1078 |
delta = _cap[in_arc];
|
| 981 |
1079 |
int result = 0;
|
| 982 |
|
Capacity d;
|
|
1080 |
Value d;
|
| 983 |
1081 |
int e;
|
| 984 |
1082 |
|
| 985 |
1083 |
// Search the cycle along the path form the first node to the root
|
| 986 |
1084 |
for (int u = first; u != join; u = _parent[u]) {
|
| 987 |
1085 |
e = _pred[u];
|
| 988 |
1086 |
d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
|
| 989 |
1087 |
if (d < delta) {
|
| 990 |
1088 |
delta = d;
|
| 991 |
1089 |
u_out = u;
|
| 992 |
1090 |
result = 1;
|
| 993 |
1091 |
}
|
| 994 |
1092 |
}
|
| 995 |
1093 |
// Search the cycle along the path form the second node to the root
|
| 996 |
1094 |
for (int u = second; u != join; u = _parent[u]) {
|
| 997 |
1095 |
e = _pred[u];
|
| 998 |
1096 |
d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
|
| 999 |
1097 |
if (d <= delta) {
|
| 1000 |
1098 |
delta = d;
|
| 1001 |
1099 |
u_out = u;
|
| 1002 |
1100 |
result = 2;
|
| 1003 |
1101 |
}
|
| 1004 |
1102 |
}
|
| 1005 |
1103 |
|
| 1006 |
1104 |
if (result == 1) {
|
| 1007 |
1105 |
u_in = first;
|
| 1008 |
1106 |
v_in = second;
|
| 1009 |
1107 |
} else {
|
| 1010 |
1108 |
u_in = second;
|
| 1011 |
1109 |
v_in = first;
|
| 1012 |
1110 |
}
|
| 1013 |
1111 |
return result != 0;
|
| 1014 |
1112 |
}
|
| 1015 |
1113 |
|
| 1016 |
1114 |
// Change _flow and _state vectors
|
| 1017 |
1115 |
void changeFlow(bool change) {
|
| 1018 |
1116 |
// Augment along the cycle
|
| 1019 |
1117 |
if (delta > 0) {
|
| 1020 |
|
Capacity val = _state[in_arc] * delta;
|
|
1118 |
Value val = _state[in_arc] * delta;
|
| 1021 |
1119 |
_flow[in_arc] += val;
|
| 1022 |
1120 |
for (int u = _source[in_arc]; u != join; u = _parent[u]) {
|
| 1023 |
1121 |
_flow[_pred[u]] += _forward[u] ? -val : val;
|
| 1024 |
1122 |
}
|
| 1025 |
1123 |
for (int u = _target[in_arc]; u != join; u = _parent[u]) {
|
| 1026 |
1124 |
_flow[_pred[u]] += _forward[u] ? val : -val;
|
| 1027 |
1125 |
}
|
| 1028 |
1126 |
}
|
| 1029 |
1127 |
// Update the state of the entering and leaving arcs
|
| 1030 |
1128 |
if (change) {
|
| 1031 |
1129 |
_state[in_arc] = STATE_TREE;
|
| 1032 |
1130 |
_state[_pred[u_out]] =
|
| 1033 |
1131 |
(_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
|
| 1034 |
1132 |
} else {
|
| 1035 |
1133 |
_state[in_arc] = -_state[in_arc];
|
| 1036 |
1134 |
}
|
| 1037 |
1135 |
}
|
| 1038 |
1136 |
|
| 1039 |
1137 |
// Update the tree structure
|
| 1040 |
1138 |
void updateTreeStructure() {
|
| 1041 |
1139 |
int u, w;
|
| 1042 |
1140 |
int old_rev_thread = _rev_thread[u_out];
|
| 1043 |
1141 |
int old_succ_num = _succ_num[u_out];
|
| 1044 |
1142 |
int old_last_succ = _last_succ[u_out];
|
| ... |
... |
@@ -1137,110 +1235,110 @@
|
| 1137 |
1235 |
if (join != old_rev_thread && v_in != old_rev_thread) {
|
| 1138 |
1236 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
|
| 1139 |
1237 |
u = _parent[u]) {
|
| 1140 |
1238 |
_last_succ[u] = old_rev_thread;
|
| 1141 |
1239 |
}
|
| 1142 |
1240 |
} else {
|
| 1143 |
1241 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
|
| 1144 |
1242 |
u = _parent[u]) {
|
| 1145 |
1243 |
_last_succ[u] = _last_succ[u_out];
|
| 1146 |
1244 |
}
|
| 1147 |
1245 |
}
|
| 1148 |
1246 |
|
| 1149 |
1247 |
// Update _succ_num from v_in to join
|
| 1150 |
1248 |
for (u = v_in; u != join; u = _parent[u]) {
|
| 1151 |
1249 |
_succ_num[u] += old_succ_num;
|
| 1152 |
1250 |
}
|
| 1153 |
1251 |
// Update _succ_num from v_out to join
|
| 1154 |
1252 |
for (u = v_out; u != join; u = _parent[u]) {
|
| 1155 |
1253 |
_succ_num[u] -= old_succ_num;
|
| 1156 |
1254 |
}
|
| 1157 |
1255 |
}
|
| 1158 |
1256 |
|
| 1159 |
1257 |
// Update potentials
|
| 1160 |
1258 |
void updatePotential() {
|
| 1161 |
|
Cost sigma = _forward[u_in] ?
|
|
1259 |
Value sigma = _forward[u_in] ?
|
| 1162 |
1260 |
_pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
|
| 1163 |
1261 |
_pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
|
| 1164 |
1262 |
if (_succ_num[u_in] > _node_num / 2) {
|
| 1165 |
1263 |
// Update in the upper subtree (which contains the root)
|
| 1166 |
1264 |
int before = _rev_thread[u_in];
|
| 1167 |
1265 |
int after = _thread[_last_succ[u_in]];
|
| 1168 |
1266 |
_thread[before] = after;
|
| 1169 |
1267 |
_pi[_root] -= sigma;
|
| 1170 |
1268 |
for (int u = _thread[_root]; u != _root; u = _thread[u]) {
|
| 1171 |
1269 |
_pi[u] -= sigma;
|
| 1172 |
1270 |
}
|
| 1173 |
1271 |
_thread[before] = u_in;
|
| 1174 |
1272 |
} else {
|
| 1175 |
1273 |
// Update in the lower subtree (which has been moved)
|
| 1176 |
1274 |
int end = _thread[_last_succ[u_in]];
|
| 1177 |
1275 |
for (int u = u_in; u != end; u = _thread[u]) {
|
| 1178 |
1276 |
_pi[u] += sigma;
|
| 1179 |
1277 |
}
|
| 1180 |
1278 |
}
|
| 1181 |
1279 |
}
|
| 1182 |
1280 |
|
| 1183 |
1281 |
// Execute the algorithm
|
| 1184 |
|
bool start(PivotRuleEnum pivot_rule) {
|
|
1282 |
bool start(PivotRule pivot_rule) {
|
| 1185 |
1283 |
// Select the pivot rule implementation
|
| 1186 |
1284 |
switch (pivot_rule) {
|
| 1187 |
|
case FIRST_ELIGIBLE_PIVOT:
|
|
1285 |
case FIRST_ELIGIBLE:
|
| 1188 |
1286 |
return start<FirstEligiblePivotRule>();
|
| 1189 |
|
case BEST_ELIGIBLE_PIVOT:
|
|
1287 |
case BEST_ELIGIBLE:
|
| 1190 |
1288 |
return start<BestEligiblePivotRule>();
|
| 1191 |
|
case BLOCK_SEARCH_PIVOT:
|
|
1289 |
case BLOCK_SEARCH:
|
| 1192 |
1290 |
return start<BlockSearchPivotRule>();
|
| 1193 |
|
case CANDIDATE_LIST_PIVOT:
|
|
1291 |
case CANDIDATE_LIST:
|
| 1194 |
1292 |
return start<CandidateListPivotRule>();
|
| 1195 |
|
case ALTERING_LIST_PIVOT:
|
|
1293 |
case ALTERING_LIST:
|
| 1196 |
1294 |
return start<AlteringListPivotRule>();
|
| 1197 |
1295 |
}
|
| 1198 |
1296 |
return false;
|
| 1199 |
1297 |
}
|
| 1200 |
1298 |
|
| 1201 |
|
template<class PivotRuleImplementation>
|
|
1299 |
template <typename PivotRuleImpl>
|
| 1202 |
1300 |
bool start() {
|
| 1203 |
|
PivotRuleImplementation pivot(*this);
|
|
1301 |
PivotRuleImpl pivot(*this);
|
| 1204 |
1302 |
|
| 1205 |
|
// Execute the network simplex algorithm
|
|
1303 |
// Execute the Network Simplex algorithm
|
| 1206 |
1304 |
while (pivot.findEnteringArc()) {
|
| 1207 |
1305 |
findJoinNode();
|
| 1208 |
1306 |
bool change = findLeavingArc();
|
| 1209 |
1307 |
changeFlow(change);
|
| 1210 |
1308 |
if (change) {
|
| 1211 |
1309 |
updateTreeStructure();
|
| 1212 |
1310 |
updatePotential();
|
| 1213 |
1311 |
}
|
| 1214 |
1312 |
}
|
| 1215 |
1313 |
|
| 1216 |
1314 |
// Check if the flow amount equals zero on all the artificial arcs
|
| 1217 |
1315 |
for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
|
| 1218 |
1316 |
if (_flow[e] > 0) return false;
|
| 1219 |
1317 |
}
|
| 1220 |
1318 |
|
| 1221 |
1319 |
// Copy flow values to _flow_map
|
| 1222 |
|
if (_orig_lower) {
|
|
1320 |
if (_plower) {
|
| 1223 |
1321 |
for (int i = 0; i != _arc_num; ++i) {
|
| 1224 |
1322 |
Arc e = _arc_ref[i];
|
| 1225 |
|
_flow_map->set(e, (*_orig_lower)[e] + _flow[i]);
|
|
1323 |
_flow_map->set(e, (*_plower)[e] + _flow[i]);
|
| 1226 |
1324 |
}
|
| 1227 |
1325 |
} else {
|
| 1228 |
1326 |
for (int i = 0; i != _arc_num; ++i) {
|
| 1229 |
1327 |
_flow_map->set(_arc_ref[i], _flow[i]);
|
| 1230 |
1328 |
}
|
| 1231 |
1329 |
}
|
| 1232 |
1330 |
// Copy potential values to _potential_map
|
| 1233 |
1331 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1234 |
1332 |
_potential_map->set(n, _pi[_node_id[n]]);
|
| 1235 |
1333 |
}
|
| 1236 |
1334 |
|
| 1237 |
1335 |
return true;
|
| 1238 |
1336 |
}
|
| 1239 |
1337 |
|
| 1240 |
1338 |
}; //class NetworkSimplex
|
| 1241 |
1339 |
|
| 1242 |
1340 |
///@}
|
| 1243 |
1341 |
|
| 1244 |
1342 |
} //namespace lemon
|
| 1245 |
1343 |
|
| 1246 |
1344 |
#endif //LEMON_NETWORK_SIMPLEX_H
|