38 //#define SORTED_LIST_PIVOT |
37 //#define SORTED_LIST_PIVOT |
39 |
38 |
40 //#define _DEBUG_ITER_ |
39 //#define _DEBUG_ITER_ |
41 |
40 |
42 |
41 |
43 /// \brief State constant for edges at their lower bounds. |
42 // State constant for edges at their lower bounds. |
44 #define LOWER 1 |
43 #define LOWER 1 |
45 /// \brief State constant for edges in the spanning tree. |
44 // State constant for edges in the spanning tree. |
46 #define TREE 0 |
45 #define TREE 0 |
47 /// \brief State constant for edges at their upper bounds. |
46 // State constant for edges at their upper bounds. |
48 #define UPPER -1 |
47 #define UPPER -1 |
49 |
48 |
50 #ifdef EDGE_BLOCK_PIVOT |
49 #ifdef EDGE_BLOCK_PIVOT |
51 #include <cmath> |
50 #include <cmath> |
52 /// \brief Lower bound for the size of blocks. |
51 #define MIN_BLOCK_SIZE 10 |
53 #define MIN_BLOCK_SIZE 10 |
|
54 #endif |
52 #endif |
55 |
53 |
56 #ifdef CANDIDATE_LIST_PIVOT |
54 #ifdef CANDIDATE_LIST_PIVOT |
57 #include <vector> |
55 #include <vector> |
58 #define LIST_LENGTH_DIV 50 |
56 #define LIST_LENGTH_DIV 50 |
59 #define MINOR_LIMIT_DIV 200 |
57 #define MINOR_LIMIT_DIV 200 |
60 #endif |
58 #endif |
61 |
59 |
62 #ifdef SORTED_LIST_PIVOT |
60 #ifdef SORTED_LIST_PIVOT |
63 #include <vector> |
61 #include <vector> |
64 #include <algorithm> |
62 #include <algorithm> |
65 #define LIST_LENGTH_DIV 100 |
63 #define LIST_LENGTH_DIV 100 |
66 #define LOWER_DIV 4 |
64 #define LOWER_DIV 4 |
67 #endif |
65 #endif |
68 |
66 |
69 namespace lemon { |
67 namespace lemon { |
70 |
68 |
71 /// \addtogroup min_cost_flow |
69 /// \addtogroup min_cost_flow |
72 /// @{ |
70 /// @{ |
73 |
71 |
74 /// \brief Implementation of the network simplex algorithm for |
72 /// \brief Implementation of the network simplex algorithm for |
75 /// finding a minimum cost flow. |
73 /// finding a minimum cost flow. |
76 /// |
74 /// |
77 /// \ref lemon::NetworkSimplex "NetworkSimplex" implements the |
75 /// \ref NetworkSimplex implements the network simplex algorithm for |
78 /// network simplex algorithm for finding a minimum cost flow. |
76 /// finding a minimum cost flow. |
79 /// |
77 /// |
80 /// \param Graph The directed graph type the algorithm runs on. |
78 /// \param Graph The directed graph type the algorithm runs on. |
81 /// \param LowerMap The type of the lower bound map. |
79 /// \param LowerMap The type of the lower bound map. |
82 /// \param CapacityMap The type of the capacity (upper bound) map. |
80 /// \param CapacityMap The type of the capacity (upper bound) map. |
83 /// \param CostMap The type of the cost (length) map. |
81 /// \param CostMap The type of the cost (length) map. |
84 /// \param SupplyMap The type of the supply map. |
82 /// \param SupplyMap The type of the supply map. |
85 /// |
83 /// |
86 /// \warning |
84 /// \warning |
87 /// - Edge capacities and costs should be nonnegative integers. |
85 /// - Edge capacities and costs should be non-negative integers. |
88 /// However \c CostMap::Value should be signed type. |
86 /// However \c CostMap::Value should be signed type. |
89 /// - Supply values should be signed integers. |
87 /// - Supply values should be signed integers. |
90 /// - \c LowerMap::Value must be convertible to |
88 /// - \c LowerMap::Value must be convertible to |
91 /// \c CapacityMap::Value and \c CapacityMap::Value must be |
89 /// \c CapacityMap::Value and \c CapacityMap::Value must be |
92 /// convertible to \c SupplyMap::Value. |
90 /// convertible to \c SupplyMap::Value. |
93 /// |
91 /// |
94 /// \author Peter Kovacs |
92 /// \author Peter Kovacs |
95 |
93 |
96 template < typename Graph, |
94 template < typename Graph, |
97 typename LowerMap = typename Graph::template EdgeMap<int>, |
95 typename LowerMap = typename Graph::template EdgeMap<int>, |
148 const SPotentialMap &pot_map; |
141 const SPotentialMap &pot_map; |
149 |
142 |
150 public: |
143 public: |
151 |
144 |
152 ReducedCostMap( const SGraph &_gr, |
145 ReducedCostMap( const SGraph &_gr, |
153 const SCostMap &_cm, |
146 const SCostMap &_cm, |
154 const SPotentialMap &_pm ) : |
147 const SPotentialMap &_pm ) : |
155 gr(_gr), cost_map(_cm), pot_map(_pm) {} |
148 gr(_gr), cost_map(_cm), pot_map(_pm) {} |
156 |
149 |
157 Cost operator[](const Edge &e) const { |
150 Cost operator[](const Edge &e) const { |
158 return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)]; |
151 return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)]; |
159 } |
152 } |
160 |
153 |
161 }; //class ReducedCostMap |
154 }; //class ReducedCostMap |
162 |
155 |
163 protected: |
156 protected: |
164 |
157 |
165 /// \brief The directed graph the algorithm runs on. |
158 /// The directed graph the algorithm runs on. |
166 SGraph graph; |
159 SGraph graph; |
167 /// \brief The original graph. |
160 /// The original graph. |
168 const Graph &graph_ref; |
161 const Graph &graph_ref; |
169 /// \brief The original lower bound map. |
162 /// The original lower bound map. |
170 const LowerMap *lower; |
163 const LowerMap *lower; |
171 /// \brief The capacity map. |
164 /// The capacity map. |
172 SCapacityMap capacity; |
165 SCapacityMap capacity; |
173 /// \brief The cost map. |
166 /// The cost map. |
174 SCostMap cost; |
167 SCostMap cost; |
175 /// \brief The supply map. |
168 /// The supply map. |
176 SSupplyMap supply; |
169 SSupplyMap supply; |
177 /// \brief The reduced cost map. |
170 /// The reduced cost map. |
178 ReducedCostMap red_cost; |
171 ReducedCostMap red_cost; |
179 /// \brief The sum of supply values equals zero. |
|
180 bool valid_supply; |
172 bool valid_supply; |
181 |
173 |
182 /// \brief The edge map of the current flow. |
174 /// The edge map of the current flow. |
183 SCapacityMap flow; |
175 SCapacityMap flow; |
184 /// \brief The edge map of the found flow on the original graph. |
176 /// The potential node map. |
|
177 SPotentialMap potential; |
185 FlowMap flow_result; |
178 FlowMap flow_result; |
186 /// \brief The potential node map. |
|
187 SPotentialMap potential; |
|
188 /// \brief The potential node map on the original graph. |
|
189 PotentialMap potential_result; |
179 PotentialMap potential_result; |
190 |
180 |
191 /// \brief Node reference for the original graph. |
181 /// Node reference for the original graph. |
192 NodeRefMap node_ref; |
182 NodeRefMap node_ref; |
193 /// \brief Edge reference for the original graph. |
183 /// Edge reference for the original graph. |
194 EdgeRefMap edge_ref; |
184 EdgeRefMap edge_ref; |
195 |
185 |
196 /// \brief The depth node map of the spanning tree structure. |
186 /// The \c depth node map of the spanning tree structure. |
197 IntNodeMap depth; |
187 IntNodeMap depth; |
198 /// \brief The parent node map of the spanning tree structure. |
188 /// The \c parent node map of the spanning tree structure. |
199 NodeNodeMap parent; |
189 NodeNodeMap parent; |
200 /// \brief The pred_edge node map of the spanning tree structure. |
190 /// The \c pred_edge node map of the spanning tree structure. |
201 EdgeNodeMap pred_edge; |
191 EdgeNodeMap pred_edge; |
202 /// \brief The thread node map of the spanning tree structure. |
192 /// The \c thread node map of the spanning tree structure. |
203 NodeNodeMap thread; |
193 NodeNodeMap thread; |
204 /// \brief The forward node map of the spanning tree structure. |
194 /// The \c forward node map of the spanning tree structure. |
205 BoolNodeMap forward; |
195 BoolNodeMap forward; |
206 /// \brief The state edge map. |
196 /// The state edge map. |
207 IntEdgeMap state; |
197 IntEdgeMap state; |
208 |
198 /// The root node of the starting spanning tree. |
|
199 Node root; |
|
200 |
|
201 // The entering edge of the current pivot iteration. |
|
202 Edge in_edge; |
|
203 // Temporary nodes used in the current pivot iteration. |
|
204 Node join, u_in, v_in, u_out, v_out; |
|
205 Node right, first, second, last; |
|
206 Node stem, par_stem, new_stem; |
|
207 // The maximum augment amount along the found cycle in the current |
|
208 // pivot iteration. |
|
209 Capacity delta; |
209 |
210 |
210 #ifdef EDGE_BLOCK_PIVOT |
211 #ifdef EDGE_BLOCK_PIVOT |
211 /// \brief The size of blocks for the "Edge Block" pivot rule. |
212 /// The size of blocks for the "Edge Block" pivot rule. |
212 int block_size; |
213 int block_size; |
213 #endif |
214 #endif |
214 #ifdef CANDIDATE_LIST_PIVOT |
215 #ifdef CANDIDATE_LIST_PIVOT |
215 /// \brief The list of candidate edges for the "Candidate List" |
216 /// \brief The list of candidate edges for the "Candidate List" |
216 /// pivot rule. |
217 /// pivot rule. |
256 /// \param _lower The lower bounds of the edges. |
245 /// \param _lower The lower bounds of the edges. |
257 /// \param _capacity The capacities (upper bounds) of the edges. |
246 /// \param _capacity The capacities (upper bounds) of the edges. |
258 /// \param _cost The cost (length) values of the edges. |
247 /// \param _cost The cost (length) values of the edges. |
259 /// \param _supply The supply values of the nodes (signed). |
248 /// \param _supply The supply values of the nodes (signed). |
260 NetworkSimplex( const Graph &_graph, |
249 NetworkSimplex( const Graph &_graph, |
261 const LowerMap &_lower, |
250 const LowerMap &_lower, |
262 const CapacityMap &_capacity, |
251 const CapacityMap &_capacity, |
263 const CostMap &_cost, |
252 const CostMap &_cost, |
264 const SupplyMap &_supply ) : |
253 const SupplyMap &_supply ) : |
265 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), |
254 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), |
266 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
255 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
267 potential_result(_graph), depth(graph), parent(graph), |
256 potential_result(_graph), depth(graph), parent(graph), |
268 pred_edge(graph), thread(graph), forward(graph), state(graph), |
257 pred_edge(graph), thread(graph), forward(graph), state(graph), |
269 node_ref(graph_ref), edge_ref(graph_ref), |
258 node_ref(graph_ref), edge_ref(graph_ref), |
270 red_cost(graph, cost, potential) |
259 red_cost(graph, cost, potential) |
271 { |
260 { |
272 // Checking the sum of supply values |
261 // Checking the sum of supply values |
273 Supply sum = 0; |
262 Supply sum = 0; |
274 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
263 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
275 sum += _supply[n]; |
264 sum += _supply[n]; |
276 if (!(valid_supply = sum == 0)) return; |
265 if (!(valid_supply = sum == 0)) return; |
277 |
266 |
278 // Copying graph_ref to graph |
267 // Copying graph_ref to graph |
279 graph.reserveNode(countNodes(graph_ref) + 1); |
268 graph.reserveNode(countNodes(graph_ref) + 1); |
280 graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref)); |
269 graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref)); |
281 copyGraph(graph, graph_ref) |
270 copyGraph(graph, graph_ref) |
282 .edgeMap(cost, _cost) |
271 .edgeMap(cost, _cost) |
283 .nodeRef(node_ref) |
272 .nodeRef(node_ref) |
284 .edgeRef(edge_ref) |
273 .edgeRef(edge_ref) |
285 .run(); |
274 .run(); |
286 |
275 |
287 // Removing nonzero lower bounds |
276 // Removing non-zero lower bounds |
288 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { |
277 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { |
289 capacity[edge_ref[e]] = _capacity[e] - _lower[e]; |
278 capacity[edge_ref[e]] = _capacity[e] - _lower[e]; |
290 } |
279 } |
291 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { |
280 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { |
292 Supply s = _supply[n]; |
281 Supply s = _supply[n]; |
293 for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) |
282 for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) |
294 s += _lower[e]; |
283 s += _lower[e]; |
295 for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) |
284 for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) |
296 s -= _lower[e]; |
285 s -= _lower[e]; |
297 supply[node_ref[n]] = s; |
286 supply[node_ref[n]] = s; |
298 } |
287 } |
299 } |
288 } |
300 |
289 |
301 /// \brief General constructor of the class (without lower bounds). |
290 /// \brief General constructor of the class (without lower bounds). |
302 /// |
291 /// |
305 /// \param _graph The directed graph the algorithm runs on. |
294 /// \param _graph The directed graph the algorithm runs on. |
306 /// \param _capacity The capacities (upper bounds) of the edges. |
295 /// \param _capacity The capacities (upper bounds) of the edges. |
307 /// \param _cost The cost (length) values of the edges. |
296 /// \param _cost The cost (length) values of the edges. |
308 /// \param _supply The supply values of the nodes (signed). |
297 /// \param _supply The supply values of the nodes (signed). |
309 NetworkSimplex( const Graph &_graph, |
298 NetworkSimplex( const Graph &_graph, |
310 const CapacityMap &_capacity, |
299 const CapacityMap &_capacity, |
311 const CostMap &_cost, |
300 const CostMap &_cost, |
312 const SupplyMap &_supply ) : |
301 const SupplyMap &_supply ) : |
313 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), |
302 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), |
314 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
303 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
315 potential_result(_graph), depth(graph), parent(graph), |
304 potential_result(_graph), depth(graph), parent(graph), |
316 pred_edge(graph), thread(graph), forward(graph), state(graph), |
305 pred_edge(graph), thread(graph), forward(graph), state(graph), |
317 node_ref(graph_ref), edge_ref(graph_ref), |
306 node_ref(graph_ref), edge_ref(graph_ref), |
318 red_cost(graph, cost, potential) |
307 red_cost(graph, cost, potential) |
319 { |
308 { |
320 // Checking the sum of supply values |
309 // Checking the sum of supply values |
321 Supply sum = 0; |
310 Supply sum = 0; |
322 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
311 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
323 sum += _supply[n]; |
312 sum += _supply[n]; |
324 if (!(valid_supply = sum == 0)) return; |
313 if (!(valid_supply = sum == 0)) return; |
325 |
314 |
326 // Copying graph_ref to graph |
315 // Copying graph_ref to graph |
327 copyGraph(graph, graph_ref) |
316 copyGraph(graph, graph_ref) |
328 .edgeMap(capacity, _capacity) |
317 .edgeMap(capacity, _capacity) |
329 .edgeMap(cost, _cost) |
318 .edgeMap(cost, _cost) |
330 .nodeMap(supply, _supply) |
319 .nodeMap(supply, _supply) |
331 .nodeRef(node_ref) |
320 .nodeRef(node_ref) |
332 .edgeRef(edge_ref) |
321 .edgeRef(edge_ref) |
333 .run(); |
322 .run(); |
334 } |
323 } |
335 |
324 |
336 /// \brief Simple constructor of the class (with lower bounds). |
325 /// \brief Simple constructor of the class (with lower bounds). |
337 /// |
326 /// |
338 /// Simple constructor of the class (with lower bounds). |
327 /// Simple constructor of the class (with lower bounds). |
342 /// \param _capacity The capacities (upper bounds) of the edges. |
331 /// \param _capacity The capacities (upper bounds) of the edges. |
343 /// \param _cost The cost (length) values of the edges. |
332 /// \param _cost The cost (length) values of the edges. |
344 /// \param _s The source node. |
333 /// \param _s The source node. |
345 /// \param _t The target node. |
334 /// \param _t The target node. |
346 /// \param _flow_value The required amount of flow from node \c _s |
335 /// \param _flow_value The required amount of flow from node \c _s |
347 /// to node \c _t (i.e. the supply of \c _s and the demand of |
336 /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t). |
348 /// \c _t). |
|
349 NetworkSimplex( const Graph &_graph, |
337 NetworkSimplex( const Graph &_graph, |
350 const LowerMap &_lower, |
338 const LowerMap &_lower, |
351 const CapacityMap &_capacity, |
339 const CapacityMap &_capacity, |
352 const CostMap &_cost, |
340 const CostMap &_cost, |
353 typename Graph::Node _s, |
341 typename Graph::Node _s, |
354 typename Graph::Node _t, |
342 typename Graph::Node _t, |
355 typename SupplyMap::Value _flow_value ) : |
343 typename SupplyMap::Value _flow_value ) : |
356 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), |
344 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph), |
357 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
345 supply(graph), flow(graph), flow_result(_graph), potential(graph), |
358 potential_result(_graph), depth(graph), parent(graph), |
346 potential_result(_graph), depth(graph), parent(graph), |
359 pred_edge(graph), thread(graph), forward(graph), state(graph), |
347 pred_edge(graph), thread(graph), forward(graph), state(graph), |
360 node_ref(graph_ref), edge_ref(graph_ref), |
348 node_ref(graph_ref), edge_ref(graph_ref), |
361 red_cost(graph, cost, potential) |
349 red_cost(graph, cost, potential) |
362 { |
350 { |
363 // Copying graph_ref to graph |
351 // Copying graph_ref to graph |
364 copyGraph(graph, graph_ref) |
352 copyGraph(graph, graph_ref) |
365 .edgeMap(cost, _cost) |
353 .edgeMap(cost, _cost) |
366 .nodeRef(node_ref) |
354 .nodeRef(node_ref) |
367 .edgeRef(edge_ref) |
355 .edgeRef(edge_ref) |
368 .run(); |
356 .run(); |
369 |
357 |
370 // Removing nonzero lower bounds |
358 // Removing non-zero lower bounds |
371 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { |
359 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) { |
372 capacity[edge_ref[e]] = _capacity[e] - _lower[e]; |
360 capacity[edge_ref[e]] = _capacity[e] - _lower[e]; |
373 } |
361 } |
374 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { |
362 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) { |
375 Supply s = 0; |
363 Supply s = 0; |
376 if (n == _s) s = _flow_value; |
364 if (n == _s) s = _flow_value; |
377 if (n == _t) s = -_flow_value; |
365 if (n == _t) s = -_flow_value; |
378 for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) |
366 for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e) |
379 s += _lower[e]; |
367 s += _lower[e]; |
380 for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) |
368 for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e) |
381 s -= _lower[e]; |
369 s -= _lower[e]; |
382 supply[node_ref[n]] = s; |
370 supply[node_ref[n]] = s; |
383 } |
371 } |
384 valid_supply = true; |
372 valid_supply = true; |
385 } |
373 } |
386 |
374 |
387 /// \brief Simple constructor of the class (without lower bounds). |
375 /// \brief Simple constructor of the class (without lower bounds). |
392 /// \param _capacity The capacities (upper bounds) of the edges. |
380 /// \param _capacity The capacities (upper bounds) of the edges. |
393 /// \param _cost The cost (length) values of the edges. |
381 /// \param _cost The cost (length) values of the edges. |
394 /// \param _s The source node. |
382 /// \param _s The source node. |
395 /// \param _t The target node. |
383 /// \param _t The target node. |
396 /// \param _flow_value The required amount of flow from node \c _s |
384 /// \param _flow_value The required amount of flow from node \c _s |
397 /// to node \c _t (i.e. the supply of \c _s and the demand of |
385 /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t). |
398 /// \c _t). |
|
399 NetworkSimplex( const Graph &_graph, |
386 NetworkSimplex( const Graph &_graph, |
400 const CapacityMap &_capacity, |
387 const CapacityMap &_capacity, |
401 const CostMap &_cost, |
388 const CostMap &_cost, |
402 typename Graph::Node _s, |
389 typename Graph::Node _s, |
403 typename Graph::Node _t, |
390 typename Graph::Node _t, |
404 typename SupplyMap::Value _flow_value ) : |
391 typename SupplyMap::Value _flow_value ) : |
405 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), |
392 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph), |
406 supply(graph, 0), flow(graph), flow_result(_graph), potential(graph), |
393 supply(graph, 0), flow(graph), flow_result(_graph), potential(graph), |
407 potential_result(_graph), depth(graph), parent(graph), |
394 potential_result(_graph), depth(graph), parent(graph), |
408 pred_edge(graph), thread(graph), forward(graph), state(graph), |
395 pred_edge(graph), thread(graph), forward(graph), state(graph), |
409 node_ref(graph_ref), edge_ref(graph_ref), |
396 node_ref(graph_ref), edge_ref(graph_ref), |
410 red_cost(graph, cost, potential) |
397 red_cost(graph, cost, potential) |
411 { |
398 { |
412 // Copying graph_ref to graph |
399 // Copying graph_ref to graph |
413 copyGraph(graph, graph_ref) |
400 copyGraph(graph, graph_ref) |
414 .edgeMap(capacity, _capacity) |
401 .edgeMap(capacity, _capacity) |
415 .edgeMap(cost, _cost) |
402 .edgeMap(cost, _cost) |
416 .nodeRef(node_ref) |
403 .nodeRef(node_ref) |
417 .edgeRef(edge_ref) |
404 .edgeRef(edge_ref) |
418 .run(); |
405 .run(); |
419 supply[node_ref[_s]] = _flow_value; |
406 supply[node_ref[_s]] = _flow_value; |
420 supply[node_ref[_t]] = -_flow_value; |
407 supply[node_ref[_t]] = -_flow_value; |
421 valid_supply = true; |
408 valid_supply = true; |
|
409 } |
|
410 |
|
411 /// \brief Runs the algorithm. |
|
412 /// |
|
413 /// Runs the algorithm. |
|
414 /// |
|
415 /// \return \c true if a feasible flow can be found. |
|
416 bool run() { |
|
417 return init() && start(); |
422 } |
418 } |
423 |
419 |
424 /// \brief Returns a const reference to the flow map. |
420 /// \brief Returns a const reference to the flow map. |
425 /// |
421 /// |
426 /// Returns a const reference to the flow map. |
422 /// Returns a const reference to the flow map. |
489 Supply sum = 0; |
476 Supply sum = 0; |
490 Node last = root; |
477 Node last = root; |
491 Edge e; |
478 Edge e; |
492 Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
479 Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
493 for (NodeIt u(graph); u != INVALID; ++u) { |
480 for (NodeIt u(graph); u != INVALID; ++u) { |
494 if (u == root) continue; |
481 if (u == root) continue; |
495 thread[last] = u; |
482 thread[last] = u; |
496 last = u; |
483 last = u; |
497 parent[u] = root; |
484 parent[u] = root; |
498 depth[u] = 1; |
485 depth[u] = 1; |
499 sum += supply[u]; |
486 sum += supply[u]; |
500 if (supply[u] >= 0) { |
487 if (supply[u] >= 0) { |
501 e = graph.addEdge(u, root); |
488 e = graph.addEdge(u, root); |
502 flow[e] = supply[u]; |
489 flow[e] = supply[u]; |
503 forward[u] = true; |
490 forward[u] = true; |
504 potential[u] = max_cost; |
491 potential[u] = max_cost; |
505 } else { |
492 } else { |
506 e = graph.addEdge(root, u); |
493 e = graph.addEdge(root, u); |
507 flow[e] = -supply[u]; |
494 flow[e] = -supply[u]; |
508 forward[u] = false; |
495 forward[u] = false; |
509 potential[u] = -max_cost; |
496 potential[u] = -max_cost; |
510 } |
497 } |
511 cost[e] = max_cost; |
498 cost[e] = max_cost; |
512 capacity[e] = std::numeric_limits<Capacity>::max(); |
499 capacity[e] = std::numeric_limits<Capacity>::max(); |
513 state[e] = TREE; |
500 state[e] = TREE; |
514 pred_edge[u] = e; |
501 pred_edge[u] = e; |
515 } |
502 } |
516 thread[last] = root; |
503 thread[last] = root; |
517 |
504 |
518 #ifdef EDGE_BLOCK_PIVOT |
505 #ifdef EDGE_BLOCK_PIVOT |
519 // Initializing block_size for the edge block pivot rule |
506 // Initializing block_size for the edge block pivot rule |
520 int edge_num = countEdges(graph); |
507 int edge_num = countEdges(graph); |
521 block_size = 2 * int(sqrt(countEdges(graph))); |
508 block_size = 2 * int(sqrt(countEdges(graph))); |
522 if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE; |
509 if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE; |
523 // block_size = edge_num >= BLOCK_NUM * MIN_BLOCK_SIZE ? |
|
524 // edge_num / BLOCK_NUM : MIN_BLOCK_SIZE; |
|
525 #endif |
510 #endif |
526 #ifdef CANDIDATE_LIST_PIVOT |
511 #ifdef CANDIDATE_LIST_PIVOT |
527 int edge_num = countEdges(graph); |
512 int edge_num = countEdges(graph); |
528 minor_count = 0; |
513 minor_count = 0; |
529 list_length = edge_num / LIST_LENGTH_DIV; |
514 list_length = edge_num / LIST_LENGTH_DIV; |
653 private: |
638 private: |
654 const IntEdgeMap &st; |
639 const IntEdgeMap &st; |
655 const ReducedCostMap &rc; |
640 const ReducedCostMap &rc; |
656 public: |
641 public: |
657 SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) : |
642 SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) : |
658 st(_st), rc(_rc) {} |
643 st(_st), rc(_rc) {} |
659 bool operator()(const Edge &e1, const Edge &e2) { |
644 bool operator()(const Edge &e1, const Edge &e2) { |
660 return st[e1] * rc[e1] < st[e2] * rc[e2]; |
645 return st[e1] * rc[e1] < st[e2] * rc[e2]; |
661 } |
646 } |
662 }; |
647 }; |
663 |
648 |
664 /// \brief Finds entering edge according to the "Sorted List" |
649 /// \brief Finds entering edge according to the "Sorted List" |
665 /// pivot rule. |
650 /// pivot rule. |
666 bool findEnteringEdge() { |
651 bool findEnteringEdge() { |
667 static SortFunc sort_func(state, red_cost); |
652 static SortFunc sort_func(state, red_cost); |
668 |
653 |
669 // Minor iteration |
654 // Minor iteration |
670 while (list_index < candidates.size()) { |
655 while (list_index < candidates.size()) { |
671 in_edge = candidates[list_index++]; |
656 in_edge = candidates[list_index++]; |
672 if (state[in_edge] * red_cost[in_edge] < 0) return true; |
657 if (state[in_edge] * red_cost[in_edge] < 0) return true; |
673 } |
658 } |
674 |
659 |
675 // Major iteration |
660 // Major iteration |
676 candidates.clear(); |
661 candidates.clear(); |
677 Cost curr, min = 0; |
662 Cost curr, min = 0; |
678 for (EdgeIt e(graph); e != INVALID; ++e) { |
663 for (EdgeIt e(graph); e != INVALID; ++e) { |
679 if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) { |
664 if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) { |
680 candidates.push_back(e); |
665 candidates.push_back(e); |
681 if (curr < min) min = curr; |
666 if (curr < min) min = curr; |
682 if (candidates.size() == list_length) break; |
667 if (candidates.size() == list_length) break; |
683 } |
668 } |
684 } |
669 } |
685 if (candidates.size() == 0) return false; |
670 if (candidates.size() == 0) return false; |
686 sort(candidates.begin(), candidates.end(), sort_func); |
671 sort(candidates.begin(), candidates.end(), sort_func); |
687 in_edge = candidates[0]; |
672 in_edge = candidates[0]; |
688 list_index = 1; |
673 list_index = 1; |
694 Node findJoinNode() { |
679 Node findJoinNode() { |
695 // Finding the join node |
680 // Finding the join node |
696 Node u = graph.source(in_edge); |
681 Node u = graph.source(in_edge); |
697 Node v = graph.target(in_edge); |
682 Node v = graph.target(in_edge); |
698 while (u != v) { |
683 while (u != v) { |
699 if (depth[u] == depth[v]) { |
684 if (depth[u] == depth[v]) { |
700 u = parent[u]; |
685 u = parent[u]; |
701 v = parent[v]; |
686 v = parent[v]; |
702 } |
687 } |
703 else if (depth[u] > depth[v]) u = parent[u]; |
688 else if (depth[u] > depth[v]) u = parent[u]; |
704 else v = parent[v]; |
689 else v = parent[v]; |
705 } |
690 } |
706 return u; |
691 return u; |
707 } |
692 } |
708 |
693 |
709 /// \brief Finds the leaving edge of the cycle. Returns \c true if |
694 /// \brief Finds the leaving edge of the cycle. Returns \c true if |
710 /// the leaving edge is not the same as the entering edge. |
695 /// the leaving edge is not the same as the entering edge. |
711 bool findLeavingEdge() { |
696 bool findLeavingEdge() { |
712 // Initializing first and second nodes according to the direction |
697 // Initializing first and second nodes according to the direction |
713 // of the cycle |
698 // of the cycle |
714 if (state[in_edge] == LOWER) { |
699 if (state[in_edge] == LOWER) { |
715 first = graph.source(in_edge); |
700 first = graph.source(in_edge); |
716 second = graph.target(in_edge); |
701 second = graph.target(in_edge); |
717 } else { |
702 } else { |
718 first = graph.target(in_edge); |
703 first = graph.target(in_edge); |
719 second = graph.source(in_edge); |
704 second = graph.source(in_edge); |
720 } |
705 } |
721 delta = capacity[in_edge]; |
706 delta = capacity[in_edge]; |
722 bool result = false; |
707 bool result = false; |
723 Capacity d; |
708 Capacity d; |
724 Edge e; |
709 Edge e; |
725 |
710 |
726 // Searching the cycle along the path form the first node to the |
711 // Searching the cycle along the path form the first node to the |
727 // root node |
712 // root node |
728 for (Node u = first; u != join; u = parent[u]) { |
713 for (Node u = first; u != join; u = parent[u]) { |
729 e = pred_edge[u]; |
714 e = pred_edge[u]; |
730 d = forward[u] ? flow[e] : capacity[e] - flow[e]; |
715 d = forward[u] ? flow[e] : capacity[e] - flow[e]; |
731 if (d < delta) { |
716 if (d < delta) { |
732 delta = d; |
717 delta = d; |
733 u_out = u; |
718 u_out = u; |
734 u_in = first; |
719 u_in = first; |
735 v_in = second; |
720 v_in = second; |
736 result = true; |
721 result = true; |
737 } |
722 } |
738 } |
723 } |
739 // Searching the cycle along the path form the second node to the |
724 // Searching the cycle along the path form the second node to the |
740 // root node |
725 // root node |
741 for (Node u = second; u != join; u = parent[u]) { |
726 for (Node u = second; u != join; u = parent[u]) { |
742 e = pred_edge[u]; |
727 e = pred_edge[u]; |
743 d = forward[u] ? capacity[e] - flow[e] : flow[e]; |
728 d = forward[u] ? capacity[e] - flow[e] : flow[e]; |
744 if (d <= delta) { |
729 if (d <= delta) { |
745 delta = d; |
730 delta = d; |
746 u_out = u; |
731 u_out = u; |
747 u_in = second; |
732 u_in = second; |
748 v_in = first; |
733 v_in = first; |
749 result = true; |
734 result = true; |
750 } |
735 } |
751 } |
736 } |
752 return result; |
737 return result; |
753 } |
738 } |
754 |
739 |
755 /// \brief Changes flow and state edge maps. |
740 /// \brief Changes \c flow and \c state edge maps. |
756 void changeFlows(bool change) { |
741 void changeFlows(bool change) { |
757 // Augmenting along the cycle |
742 // Augmenting along the cycle |
758 if (delta > 0) { |
743 if (delta > 0) { |
759 Capacity val = state[in_edge] * delta; |
744 Capacity val = state[in_edge] * delta; |
760 flow[in_edge] += val; |
745 flow[in_edge] += val; |
761 for (Node u = graph.source(in_edge); u != join; u = parent[u]) { |
746 for (Node u = graph.source(in_edge); u != join; u = parent[u]) { |
762 flow[pred_edge[u]] += forward[u] ? -val : val; |
747 flow[pred_edge[u]] += forward[u] ? -val : val; |
763 } |
748 } |
764 for (Node u = graph.target(in_edge); u != join; u = parent[u]) { |
749 for (Node u = graph.target(in_edge); u != join; u = parent[u]) { |
765 flow[pred_edge[u]] += forward[u] ? val : -val; |
750 flow[pred_edge[u]] += forward[u] ? val : -val; |
766 } |
751 } |
767 } |
752 } |
768 // Updating the state of the entering and leaving edges |
753 // Updating the state of the entering and leaving edges |
769 if (change) { |
754 if (change) { |
770 state[in_edge] = TREE; |
755 state[in_edge] = TREE; |
771 state[pred_edge[u_out]] = |
756 state[pred_edge[u_out]] = |
772 (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER; |
757 (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER; |
773 } else { |
758 } else { |
774 state[in_edge] = -state[in_edge]; |
759 state[in_edge] = -state[in_edge]; |
775 } |
760 } |
776 } |
761 } |
777 |
762 |
778 /// \brief Updates thread and parent node maps. |
763 /// \brief Updates \c thread and \c parent node maps. |
779 void updateThreadParent() { |
764 void updateThreadParent() { |
780 Node u; |
765 Node u; |
781 v_out = parent[u_out]; |
766 v_out = parent[u_out]; |
782 |
767 |
783 // Handling the case when join and v_out coincide |
768 // Handling the case when join and v_out coincide |
784 bool par_first = false; |
769 bool par_first = false; |
785 if (join == v_out) { |
770 if (join == v_out) { |
786 for (u = join; u != u_in && u != v_in; u = thread[u]) ; |
771 for (u = join; u != u_in && u != v_in; u = thread[u]) ; |
787 if (u == v_in) { |
772 if (u == v_in) { |
788 par_first = true; |
773 par_first = true; |
789 while (thread[u] != u_out) u = thread[u]; |
774 while (thread[u] != u_out) u = thread[u]; |
790 first = u; |
775 first = u; |
791 } |
776 } |
792 } |
777 } |
793 |
778 |
794 // Finding the last successor of u_in (u) and the node after it |
779 // Finding the last successor of u_in (u) and the node after it |
795 // (right) according to the thread index |
780 // (right) according to the thread index |
796 for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ; |
781 for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ; |
797 right = thread[u]; |
782 right = thread[u]; |
798 if (thread[v_in] == u_out) { |
783 if (thread[v_in] == u_out) { |
799 for (last = u; depth[last] > depth[u_out]; last = thread[last]) ; |
784 for (last = u; depth[last] > depth[u_out]; last = thread[last]) ; |
800 if (last == u_out) last = thread[last]; |
785 if (last == u_out) last = thread[last]; |
801 } |
786 } |
802 else last = thread[v_in]; |
787 else last = thread[v_in]; |
803 |
788 |
804 // Updating stem nodes |
789 // Updating stem nodes |
805 thread[v_in] = stem = u_in; |
790 thread[v_in] = stem = u_in; |
806 par_stem = v_in; |
791 par_stem = v_in; |
807 while (stem != u_out) { |
792 while (stem != u_out) { |
808 thread[u] = new_stem = parent[stem]; |
793 thread[u] = new_stem = parent[stem]; |
809 |
794 |
810 // Finding the node just before the stem node (u) according to |
795 // Finding the node just before the stem node (u) according to |
811 // the original thread index |
796 // the original thread index |
812 for (u = new_stem; thread[u] != stem; u = thread[u]) ; |
797 for (u = new_stem; thread[u] != stem; u = thread[u]) ; |
813 thread[u] = right; |
798 thread[u] = right; |
814 |
799 |
815 // Changing the parent node of stem and shifting stem and |
800 // Changing the parent node of stem and shifting stem and |
816 // par_stem nodes |
801 // par_stem nodes |
817 parent[stem] = par_stem; |
802 parent[stem] = par_stem; |
818 par_stem = stem; |
803 par_stem = stem; |
819 stem = new_stem; |
804 stem = new_stem; |
820 |
805 |
821 // Finding the last successor of stem (u) and the node after it |
806 // Finding the last successor of stem (u) and the node after it |
822 // (right) according to the thread index |
807 // (right) according to the thread index |
823 for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ; |
808 for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ; |
824 right = thread[u]; |
809 right = thread[u]; |
825 } |
810 } |
826 parent[u_out] = par_stem; |
811 parent[u_out] = par_stem; |
827 thread[u] = last; |
812 thread[u] = last; |
828 |
813 |
829 if (join == v_out && par_first) { |
814 if (join == v_out && par_first) { |
830 if (first != v_in) thread[first] = right; |
815 if (first != v_in) thread[first] = right; |
831 } else { |
816 } else { |
832 for (u = v_out; thread[u] != u_out; u = thread[u]) ; |
817 for (u = v_out; thread[u] != u_out; u = thread[u]) ; |
833 thread[u] = right; |
818 thread[u] = right; |
834 } |
819 } |
835 } |
820 } |
836 |
821 |
837 /// \brief Updates pred_edge and forward node maps. |
822 /// \brief Updates \c pred_edge and \c forward node maps. |
838 void updatePredEdge() { |
823 void updatePredEdge() { |
839 Node u = u_out, v; |
824 Node u = u_out, v; |
840 while (u != u_in) { |
825 while (u != u_in) { |
841 v = parent[u]; |
826 v = parent[u]; |
842 pred_edge[u] = pred_edge[v]; |
827 pred_edge[u] = pred_edge[v]; |
843 forward[u] = !forward[v]; |
828 forward[u] = !forward[v]; |
844 u = v; |
829 u = v; |
845 } |
830 } |
846 pred_edge[u_in] = in_edge; |
831 pred_edge[u_in] = in_edge; |
847 forward[u_in] = (u_in == graph.source(in_edge)); |
832 forward[u_in] = (u_in == graph.source(in_edge)); |
848 } |
833 } |
849 |
834 |
850 /// \brief Updates depth and potential node maps. |
835 /// \brief Updates \c depth and \c potential node maps. |
851 void updateDepthPotential() { |
836 void updateDepthPotential() { |
852 depth[u_in] = depth[v_in] + 1; |
837 depth[u_in] = depth[v_in] + 1; |
853 potential[u_in] = forward[u_in] ? |
838 potential[u_in] = forward[u_in] ? |
854 potential[v_in] + cost[pred_edge[u_in]] : |
839 potential[v_in] + cost[pred_edge[u_in]] : |
855 potential[v_in] - cost[pred_edge[u_in]]; |
840 potential[v_in] - cost[pred_edge[u_in]]; |
856 |
841 |
857 Node u = thread[u_in], v; |
842 Node u = thread[u_in], v; |
858 while (true) { |
843 while (true) { |
859 v = parent[u]; |
844 v = parent[u]; |
860 if (v == INVALID) break; |
845 if (v == INVALID) break; |
861 depth[u] = depth[v] + 1; |
846 depth[u] = depth[v] + 1; |
862 potential[u] = forward[u] ? |
847 potential[u] = forward[u] ? |
863 potential[v] + cost[pred_edge[u]] : |
848 potential[v] + cost[pred_edge[u]] : |
864 potential[v] - cost[pred_edge[u]]; |
849 potential[v] - cost[pred_edge[u]]; |
865 if (depth[u] <= depth[v_in]) break; |
850 if (depth[u] <= depth[v_in]) break; |
866 u = thread[u]; |
851 u = thread[u]; |
867 } |
852 } |
868 } |
853 } |
869 |
854 |
870 /// \brief Executes the algorithm. |
855 /// \brief Executes the algorithm. |
871 bool start() { |
856 bool start() { |
878 while (findEnteringEdge(next_edge)) |
863 while (findEnteringEdge(next_edge)) |
879 #else |
864 #else |
880 while (findEnteringEdge()) |
865 while (findEnteringEdge()) |
881 #endif |
866 #endif |
882 { |
867 { |
883 join = findJoinNode(); |
868 join = findJoinNode(); |
884 bool change = findLeavingEdge(); |
869 bool change = findLeavingEdge(); |
885 changeFlows(change); |
870 changeFlows(change); |
886 if (change) { |
871 if (change) { |
887 updateThreadParent(); |
872 updateThreadParent(); |
888 updatePredEdge(); |
873 updatePredEdge(); |
889 updateDepthPotential(); |
874 updateDepthPotential(); |
890 } |
875 } |
891 #ifdef _DEBUG_ITER_ |
876 #ifdef _DEBUG_ITER_ |
892 ++iter_num; |
877 ++iter_num; |
893 #endif |
878 #endif |
894 } |
879 } |
895 |
880 |
896 #ifdef _DEBUG_ITER_ |
881 #ifdef _DEBUG_ITER_ |
897 std::cout << "Network Simplex algorithm finished. " << iter_num |
882 std::cout << "Network Simplex algorithm finished. " << iter_num |
898 << " pivot iterations performed." << std::endl; |
883 << " pivot iterations performed." << std::endl; |
899 #endif |
884 #endif |
900 |
885 |
901 // Checking if the flow amount equals zero on all the |
886 // Checking if the flow amount equals zero on all the |
902 // artificial edges |
887 // artificial edges |
903 for (InEdgeIt e(graph, root); e != INVALID; ++e) |
888 for (InEdgeIt e(graph, root); e != INVALID; ++e) |
904 if (flow[e] > 0) return false; |
889 if (flow[e] > 0) return false; |
905 for (OutEdgeIt e(graph, root); e != INVALID; ++e) |
890 for (OutEdgeIt e(graph, root); e != INVALID; ++e) |
906 if (flow[e] > 0) return false; |
891 if (flow[e] > 0) return false; |
907 |
892 |
908 // Copying flow values to flow_result |
893 // Copying flow values to flow_result |
909 if (lower) { |
894 if (lower) { |
910 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) |
895 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) |
911 flow_result[e] = (*lower)[e] + flow[edge_ref[e]]; |
896 flow_result[e] = (*lower)[e] + flow[edge_ref[e]]; |
912 } else { |
897 } else { |
913 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) |
898 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) |
914 flow_result[e] = flow[edge_ref[e]]; |
899 flow_result[e] = flow[edge_ref[e]]; |
915 } |
900 } |
916 // Copying potential values to potential_result |
901 // Copying potential values to potential_result |
917 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
902 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) |
918 potential_result[n] = potential[node_ref[n]]; |
903 potential_result[n] = potential[node_ref[n]]; |
919 |
904 |
920 return true; |
905 return true; |
921 } |
906 } |
922 |
907 |
923 }; //class NetworkSimplex |
908 }; //class NetworkSimplex |