| ... |
... |
@@ -16,13 +16,13 @@
|
| 16 |
16 |
*
|
| 17 |
17 |
*/
|
| 18 |
18 |
|
| 19 |
19 |
#ifndef LEMON_NETWORK_SIMPLEX_H
|
| 20 |
20 |
#define LEMON_NETWORK_SIMPLEX_H
|
| 21 |
21 |
|
| 22 |
|
/// \ingroup min_cost_flow
|
|
22 |
/// \ingroup min_cost_flow_algs
|
| 23 |
23 |
///
|
| 24 |
24 |
/// \file
|
| 25 |
25 |
/// \brief Network Simplex algorithm for finding a minimum cost flow.
|
| 26 |
26 |
|
| 27 |
27 |
#include <vector>
|
| 28 |
28 |
#include <limits>
|
| ... |
... |
@@ -30,13 +30,13 @@
|
| 30 |
30 |
|
| 31 |
31 |
#include <lemon/core.h>
|
| 32 |
32 |
#include <lemon/math.h>
|
| 33 |
33 |
|
| 34 |
34 |
namespace lemon {
|
| 35 |
35 |
|
| 36 |
|
/// \addtogroup min_cost_flow
|
|
36 |
/// \addtogroup min_cost_flow_algs
|
| 37 |
37 |
/// @{
|
| 38 |
38 |
|
| 39 |
39 |
/// \brief Implementation of the primal Network Simplex algorithm
|
| 40 |
40 |
/// for finding a \ref min_cost_flow "minimum cost flow".
|
| 41 |
41 |
///
|
| 42 |
42 |
/// \ref NetworkSimplex implements the primal Network Simplex algorithm
|
| ... |
... |
@@ -99,56 +99,22 @@
|
| 99 |
99 |
/// \brief Constants for selecting the type of the supply constraints.
|
| 100 |
100 |
///
|
| 101 |
101 |
/// Enum type containing constants for selecting the supply type,
|
| 102 |
102 |
/// i.e. the direction of the inequalities in the supply/demand
|
| 103 |
103 |
/// constraints of the \ref min_cost_flow "minimum cost flow problem".
|
| 104 |
104 |
///
|
| 105 |
|
/// The default supply type is \c GEQ, since this form is supported
|
| 106 |
|
/// by other minimum cost flow algorithms and the \ref Circulation
|
| 107 |
|
/// algorithm, as well.
|
| 108 |
|
/// The \c LEQ problem type can be selected using the \ref supplyType()
|
| 109 |
|
/// function.
|
| 110 |
|
///
|
| 111 |
|
/// Note that the equality form is a special case of both supply types.
|
|
105 |
/// The default supply type is \c GEQ, the \c LEQ type can be
|
|
106 |
/// selected using \ref supplyType().
|
|
107 |
/// The equality form is a special case of both supply types.
|
| 112 |
108 |
enum SupplyType {
|
| 113 |
|
|
| 114 |
109 |
/// This option means that there are <em>"greater or equal"</em>
|
| 115 |
|
/// supply/demand constraints in the definition, i.e. the exact
|
| 116 |
|
/// formulation of the problem is the following.
|
| 117 |
|
/**
|
| 118 |
|
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
|
| 119 |
|
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
|
| 120 |
|
sup(u) \quad \forall u\in V \f]
|
| 121 |
|
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
|
| 122 |
|
*/
|
| 123 |
|
/// It means that the total demand must be greater or equal to the
|
| 124 |
|
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
|
| 125 |
|
/// negative) and all the supplies have to be carried out from
|
| 126 |
|
/// the supply nodes, but there could be demands that are not
|
| 127 |
|
/// satisfied.
|
|
110 |
/// supply/demand constraints in the definition of the problem.
|
| 128 |
111 |
GEQ,
|
| 129 |
|
/// It is just an alias for the \c GEQ option.
|
| 130 |
|
CARRY_SUPPLIES = GEQ,
|
| 131 |
|
|
| 132 |
112 |
/// This option means that there are <em>"less or equal"</em>
|
| 133 |
|
/// supply/demand constraints in the definition, i.e. the exact
|
| 134 |
|
/// formulation of the problem is the following.
|
| 135 |
|
/**
|
| 136 |
|
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
|
| 137 |
|
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
|
| 138 |
|
sup(u) \quad \forall u\in V \f]
|
| 139 |
|
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
|
| 140 |
|
*/
|
| 141 |
|
/// It means that the total demand must be less or equal to the
|
| 142 |
|
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
|
| 143 |
|
/// positive) and all the demands have to be satisfied, but there
|
| 144 |
|
/// could be supplies that are not carried out from the supply
|
| 145 |
|
/// nodes.
|
| 146 |
|
LEQ,
|
| 147 |
|
/// It is just an alias for the \c LEQ option.
|
| 148 |
|
SATISFY_DEMANDS = LEQ
|
|
113 |
/// supply/demand constraints in the definition of the problem.
|
|
114 |
LEQ
|
| 149 |
115 |
};
|
| 150 |
116 |
|
| 151 |
117 |
/// \brief Constants for selecting the pivot rule.
|
| 152 |
118 |
///
|
| 153 |
119 |
/// Enum type containing constants for selecting the pivot rule for
|
| 154 |
120 |
/// the \ref run() function.
|
| ... |
... |
@@ -212,12 +178,14 @@
|
| 212 |
178 |
private:
|
| 213 |
179 |
|
| 214 |
180 |
// Data related to the underlying digraph
|
| 215 |
181 |
const GR &_graph;
|
| 216 |
182 |
int _node_num;
|
| 217 |
183 |
int _arc_num;
|
|
184 |
int _all_arc_num;
|
|
185 |
int _search_arc_num;
|
| 218 |
186 |
|
| 219 |
187 |
// Parameters of the problem
|
| 220 |
188 |
bool _have_lower;
|
| 221 |
189 |
SupplyType _stype;
|
| 222 |
190 |
Value _sum_supply;
|
| 223 |
191 |
|
| ... |
... |
@@ -274,30 +242,31 @@
|
| 274 |
242 |
const IntVector &_source;
|
| 275 |
243 |
const IntVector &_target;
|
| 276 |
244 |
const CostVector &_cost;
|
| 277 |
245 |
const IntVector &_state;
|
| 278 |
246 |
const CostVector &_pi;
|
| 279 |
247 |
int &_in_arc;
|
| 280 |
|
int _arc_num;
|
|
248 |
int _search_arc_num;
|
| 281 |
249 |
|
| 282 |
250 |
// Pivot rule data
|
| 283 |
251 |
int _next_arc;
|
| 284 |
252 |
|
| 285 |
253 |
public:
|
| 286 |
254 |
|
| 287 |
255 |
// Constructor
|
| 288 |
256 |
FirstEligiblePivotRule(NetworkSimplex &ns) :
|
| 289 |
257 |
_source(ns._source), _target(ns._target),
|
| 290 |
258 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 291 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
|
259 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
|
|
260 |
_next_arc(0)
|
| 292 |
261 |
{}
|
| 293 |
262 |
|
| 294 |
263 |
// Find next entering arc
|
| 295 |
264 |
bool findEnteringArc() {
|
| 296 |
265 |
Cost c;
|
| 297 |
|
for (int e = _next_arc; e < _arc_num; ++e) {
|
|
266 |
for (int e = _next_arc; e < _search_arc_num; ++e) {
|
| 298 |
267 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 299 |
268 |
if (c < 0) {
|
| 300 |
269 |
_in_arc = e;
|
| 301 |
270 |
_next_arc = e + 1;
|
| 302 |
271 |
return true;
|
| 303 |
272 |
}
|
| ... |
... |
@@ -325,27 +294,27 @@
|
| 325 |
294 |
const IntVector &_source;
|
| 326 |
295 |
const IntVector &_target;
|
| 327 |
296 |
const CostVector &_cost;
|
| 328 |
297 |
const IntVector &_state;
|
| 329 |
298 |
const CostVector &_pi;
|
| 330 |
299 |
int &_in_arc;
|
| 331 |
|
int _arc_num;
|
|
300 |
int _search_arc_num;
|
| 332 |
301 |
|
| 333 |
302 |
public:
|
| 334 |
303 |
|
| 335 |
304 |
// Constructor
|
| 336 |
305 |
BestEligiblePivotRule(NetworkSimplex &ns) :
|
| 337 |
306 |
_source(ns._source), _target(ns._target),
|
| 338 |
307 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 339 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num)
|
|
308 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
|
| 340 |
309 |
{}
|
| 341 |
310 |
|
| 342 |
311 |
// Find next entering arc
|
| 343 |
312 |
bool findEnteringArc() {
|
| 344 |
313 |
Cost c, min = 0;
|
| 345 |
|
for (int e = 0; e < _arc_num; ++e) {
|
|
314 |
for (int e = 0; e < _search_arc_num; ++e) {
|
| 346 |
315 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 347 |
316 |
if (c < min) {
|
| 348 |
317 |
min = c;
|
| 349 |
318 |
_in_arc = e;
|
| 350 |
319 |
}
|
| 351 |
320 |
}
|
| ... |
... |
@@ -364,41 +333,42 @@
|
| 364 |
333 |
const IntVector &_source;
|
| 365 |
334 |
const IntVector &_target;
|
| 366 |
335 |
const CostVector &_cost;
|
| 367 |
336 |
const IntVector &_state;
|
| 368 |
337 |
const CostVector &_pi;
|
| 369 |
338 |
int &_in_arc;
|
| 370 |
|
int _arc_num;
|
|
339 |
int _search_arc_num;
|
| 371 |
340 |
|
| 372 |
341 |
// Pivot rule data
|
| 373 |
342 |
int _block_size;
|
| 374 |
343 |
int _next_arc;
|
| 375 |
344 |
|
| 376 |
345 |
public:
|
| 377 |
346 |
|
| 378 |
347 |
// Constructor
|
| 379 |
348 |
BlockSearchPivotRule(NetworkSimplex &ns) :
|
| 380 |
349 |
_source(ns._source), _target(ns._target),
|
| 381 |
350 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 382 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
|
351 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
|
|
352 |
_next_arc(0)
|
| 383 |
353 |
{
|
| 384 |
354 |
// The main parameters of the pivot rule
|
| 385 |
|
const double BLOCK_SIZE_FACTOR = 2.0;
|
|
355 |
const double BLOCK_SIZE_FACTOR = 0.5;
|
| 386 |
356 |
const int MIN_BLOCK_SIZE = 10;
|
| 387 |
357 |
|
| 388 |
358 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR *
|
| 389 |
|
std::sqrt(double(_arc_num))),
|
|
359 |
std::sqrt(double(_search_arc_num))),
|
| 390 |
360 |
MIN_BLOCK_SIZE );
|
| 391 |
361 |
}
|
| 392 |
362 |
|
| 393 |
363 |
// Find next entering arc
|
| 394 |
364 |
bool findEnteringArc() {
|
| 395 |
365 |
Cost c, min = 0;
|
| 396 |
366 |
int cnt = _block_size;
|
| 397 |
367 |
int e, min_arc = _next_arc;
|
| 398 |
|
for (e = _next_arc; e < _arc_num; ++e) {
|
|
368 |
for (e = _next_arc; e < _search_arc_num; ++e) {
|
| 399 |
369 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 400 |
370 |
if (c < min) {
|
| 401 |
371 |
min = c;
|
| 402 |
372 |
min_arc = e;
|
| 403 |
373 |
}
|
| 404 |
374 |
if (--cnt == 0) {
|
| ... |
... |
@@ -437,13 +407,13 @@
|
| 437 |
407 |
const IntVector &_source;
|
| 438 |
408 |
const IntVector &_target;
|
| 439 |
409 |
const CostVector &_cost;
|
| 440 |
410 |
const IntVector &_state;
|
| 441 |
411 |
const CostVector &_pi;
|
| 442 |
412 |
int &_in_arc;
|
| 443 |
|
int _arc_num;
|
|
413 |
int _search_arc_num;
|
| 444 |
414 |
|
| 445 |
415 |
// Pivot rule data
|
| 446 |
416 |
IntVector _candidates;
|
| 447 |
417 |
int _list_length, _minor_limit;
|
| 448 |
418 |
int _curr_length, _minor_count;
|
| 449 |
419 |
int _next_arc;
|
| ... |
... |
@@ -451,22 +421,23 @@
|
| 451 |
421 |
public:
|
| 452 |
422 |
|
| 453 |
423 |
/// Constructor
|
| 454 |
424 |
CandidateListPivotRule(NetworkSimplex &ns) :
|
| 455 |
425 |
_source(ns._source), _target(ns._target),
|
| 456 |
426 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 457 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
|
427 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
|
|
428 |
_next_arc(0)
|
| 458 |
429 |
{
|
| 459 |
430 |
// The main parameters of the pivot rule
|
| 460 |
431 |
const double LIST_LENGTH_FACTOR = 1.0;
|
| 461 |
432 |
const int MIN_LIST_LENGTH = 10;
|
| 462 |
433 |
const double MINOR_LIMIT_FACTOR = 0.1;
|
| 463 |
434 |
const int MIN_MINOR_LIMIT = 3;
|
| 464 |
435 |
|
| 465 |
436 |
_list_length = std::max( int(LIST_LENGTH_FACTOR *
|
| 466 |
|
std::sqrt(double(_arc_num))),
|
|
437 |
std::sqrt(double(_search_arc_num))),
|
| 467 |
438 |
MIN_LIST_LENGTH );
|
| 468 |
439 |
_minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
|
| 469 |
440 |
MIN_MINOR_LIMIT );
|
| 470 |
441 |
_curr_length = _minor_count = 0;
|
| 471 |
442 |
_candidates.resize(_list_length);
|
| 472 |
443 |
}
|
| ... |
... |
@@ -497,13 +468,13 @@
|
| 497 |
468 |
}
|
| 498 |
469 |
}
|
| 499 |
470 |
|
| 500 |
471 |
// Major iteration: build a new candidate list
|
| 501 |
472 |
min = 0;
|
| 502 |
473 |
_curr_length = 0;
|
| 503 |
|
for (e = _next_arc; e < _arc_num; ++e) {
|
|
474 |
for (e = _next_arc; e < _search_arc_num; ++e) {
|
| 504 |
475 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 505 |
476 |
if (c < 0) {
|
| 506 |
477 |
_candidates[_curr_length++] = e;
|
| 507 |
478 |
if (c < min) {
|
| 508 |
479 |
min = c;
|
| 509 |
480 |
min_arc = e;
|
| ... |
... |
@@ -543,13 +514,13 @@
|
| 543 |
514 |
const IntVector &_source;
|
| 544 |
515 |
const IntVector &_target;
|
| 545 |
516 |
const CostVector &_cost;
|
| 546 |
517 |
const IntVector &_state;
|
| 547 |
518 |
const CostVector &_pi;
|
| 548 |
519 |
int &_in_arc;
|
| 549 |
|
int _arc_num;
|
|
520 |
int _search_arc_num;
|
| 550 |
521 |
|
| 551 |
522 |
// Pivot rule data
|
| 552 |
523 |
int _block_size, _head_length, _curr_length;
|
| 553 |
524 |
int _next_arc;
|
| 554 |
525 |
IntVector _candidates;
|
| 555 |
526 |
CostVector _cand_cost;
|
| ... |
... |
@@ -571,23 +542,23 @@
|
| 571 |
542 |
public:
|
| 572 |
543 |
|
| 573 |
544 |
// Constructor
|
| 574 |
545 |
AlteringListPivotRule(NetworkSimplex &ns) :
|
| 575 |
546 |
_source(ns._source), _target(ns._target),
|
| 576 |
547 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
| 577 |
|
_in_arc(ns.in_arc), _arc_num(ns._arc_num),
|
| 578 |
|
_next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
|
|
548 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
|
|
549 |
_next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
|
| 579 |
550 |
{
|
| 580 |
551 |
// The main parameters of the pivot rule
|
| 581 |
552 |
const double BLOCK_SIZE_FACTOR = 1.5;
|
| 582 |
553 |
const int MIN_BLOCK_SIZE = 10;
|
| 583 |
554 |
const double HEAD_LENGTH_FACTOR = 0.1;
|
| 584 |
555 |
const int MIN_HEAD_LENGTH = 3;
|
| 585 |
556 |
|
| 586 |
557 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR *
|
| 587 |
|
std::sqrt(double(_arc_num))),
|
|
558 |
std::sqrt(double(_search_arc_num))),
|
| 588 |
559 |
MIN_BLOCK_SIZE );
|
| 589 |
560 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
|
| 590 |
561 |
MIN_HEAD_LENGTH );
|
| 591 |
562 |
_candidates.resize(_head_length + _block_size);
|
| 592 |
563 |
_curr_length = 0;
|
| 593 |
564 |
}
|
| ... |
... |
@@ -607,13 +578,13 @@
|
| 607 |
578 |
|
| 608 |
579 |
// Extend the list
|
| 609 |
580 |
int cnt = _block_size;
|
| 610 |
581 |
int last_arc = 0;
|
| 611 |
582 |
int limit = _head_length;
|
| 612 |
583 |
|
| 613 |
|
for (int e = _next_arc; e < _arc_num; ++e) {
|
|
584 |
for (int e = _next_arc; e < _search_arc_num; ++e) {
|
| 614 |
585 |
_cand_cost[e] = _state[e] *
|
| 615 |
586 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
| 616 |
587 |
if (_cand_cost[e] < 0) {
|
| 617 |
588 |
_candidates[_curr_length++] = e;
|
| 618 |
589 |
last_arc = e;
|
| 619 |
590 |
}
|
| ... |
... |
@@ -675,33 +646,33 @@
|
| 675 |
646 |
"The cost type of NetworkSimplex must be signed");
|
| 676 |
647 |
|
| 677 |
648 |
// Resize vectors
|
| 678 |
649 |
_node_num = countNodes(_graph);
|
| 679 |
650 |
_arc_num = countArcs(_graph);
|
| 680 |
651 |
int all_node_num = _node_num + 1;
|
| 681 |
|
int all_arc_num = _arc_num + _node_num;
|
|
652 |
int max_arc_num = _arc_num + 2 * _node_num;
|
| 682 |
653 |
|
| 683 |
|
_source.resize(all_arc_num);
|
| 684 |
|
_target.resize(all_arc_num);
|
|
654 |
_source.resize(max_arc_num);
|
|
655 |
_target.resize(max_arc_num);
|
| 685 |
656 |
|
| 686 |
|
_lower.resize(all_arc_num);
|
| 687 |
|
_upper.resize(all_arc_num);
|
| 688 |
|
_cap.resize(all_arc_num);
|
| 689 |
|
_cost.resize(all_arc_num);
|
|
657 |
_lower.resize(_arc_num);
|
|
658 |
_upper.resize(_arc_num);
|
|
659 |
_cap.resize(max_arc_num);
|
|
660 |
_cost.resize(max_arc_num);
|
| 690 |
661 |
_supply.resize(all_node_num);
|
| 691 |
|
_flow.resize(all_arc_num);
|
|
662 |
_flow.resize(max_arc_num);
|
| 692 |
663 |
_pi.resize(all_node_num);
|
| 693 |
664 |
|
| 694 |
665 |
_parent.resize(all_node_num);
|
| 695 |
666 |
_pred.resize(all_node_num);
|
| 696 |
667 |
_forward.resize(all_node_num);
|
| 697 |
668 |
_thread.resize(all_node_num);
|
| 698 |
669 |
_rev_thread.resize(all_node_num);
|
| 699 |
670 |
_succ_num.resize(all_node_num);
|
| 700 |
671 |
_last_succ.resize(all_node_num);
|
| 701 |
|
_state.resize(all_arc_num);
|
|
672 |
_state.resize(max_arc_num);
|
| 702 |
673 |
|
| 703 |
674 |
// Copy the graph (store the arcs in a mixed order)
|
| 704 |
675 |
int i = 0;
|
| 705 |
676 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
|
| 706 |
677 |
_node_id[n] = i;
|
| 707 |
678 |
}
|
| ... |
... |
@@ -1066,13 +1037,13 @@
|
| 1066 |
1037 |
}
|
| 1067 |
1038 |
}
|
| 1068 |
1039 |
|
| 1069 |
1040 |
// Initialize artifical cost
|
| 1070 |
1041 |
Cost ART_COST;
|
| 1071 |
1042 |
if (std::numeric_limits<Cost>::is_exact) {
|
| 1072 |
|
ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
|
|
1043 |
ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
|
| 1073 |
1044 |
} else {
|
| 1074 |
1045 |
ART_COST = std::numeric_limits<Cost>::min();
|
| 1075 |
1046 |
for (int i = 0; i != _arc_num; ++i) {
|
| 1076 |
1047 |
if (_cost[i] > ART_COST) ART_COST = _cost[i];
|
| 1077 |
1048 |
}
|
| 1078 |
1049 |
ART_COST = (ART_COST + 1) * _node_num;
|
| ... |
... |
@@ -1090,35 +1061,127 @@
|
| 1090 |
1061 |
_pred[_root] = -1;
|
| 1091 |
1062 |
_thread[_root] = 0;
|
| 1092 |
1063 |
_rev_thread[0] = _root;
|
| 1093 |
1064 |
_succ_num[_root] = _node_num + 1;
|
| 1094 |
1065 |
_last_succ[_root] = _root - 1;
|
| 1095 |
1066 |
_supply[_root] = -_sum_supply;
|
| 1096 |
|
_pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
|
|
1067 |
_pi[_root] = 0;
|
| 1097 |
1068 |
|
| 1098 |
1069 |
// Add artificial arcs and initialize the spanning tree data structure
|
| 1099 |
|
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
| 1100 |
|
_parent[u] = _root;
|
| 1101 |
|
_pred[u] = e;
|
| 1102 |
|
_thread[u] = u + 1;
|
| 1103 |
|
_rev_thread[u + 1] = u;
|
| 1104 |
|
_succ_num[u] = 1;
|
| 1105 |
|
_last_succ[u] = u;
|
| 1106 |
|
_cost[e] = ART_COST;
|
| 1107 |
|
_cap[e] = INF;
|
| 1108 |
|
_state[e] = STATE_TREE;
|
| 1109 |
|
if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
|
| 1110 |
|
_flow[e] = _supply[u];
|
| 1111 |
|
_forward[u] = true;
|
| 1112 |
|
_pi[u] = -ART_COST + _pi[_root];
|
| 1113 |
|
} else {
|
| 1114 |
|
_flow[e] = -_supply[u];
|
| 1115 |
|
_forward[u] = false;
|
| 1116 |
|
_pi[u] = ART_COST + _pi[_root];
|
|
1070 |
if (_sum_supply == 0) {
|
|
1071 |
// EQ supply constraints
|
|
1072 |
_search_arc_num = _arc_num;
|
|
1073 |
_all_arc_num = _arc_num + _node_num;
|
|
1074 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
|
1075 |
_parent[u] = _root;
|
|
1076 |
_pred[u] = e;
|
|
1077 |
_thread[u] = u + 1;
|
|
1078 |
_rev_thread[u + 1] = u;
|
|
1079 |
_succ_num[u] = 1;
|
|
1080 |
_last_succ[u] = u;
|
|
1081 |
_cap[e] = INF;
|
|
1082 |
_state[e] = STATE_TREE;
|
|
1083 |
if (_supply[u] >= 0) {
|
|
1084 |
_forward[u] = true;
|
|
1085 |
_pi[u] = 0;
|
|
1086 |
_source[e] = u;
|
|
1087 |
_target[e] = _root;
|
|
1088 |
_flow[e] = _supply[u];
|
|
1089 |
_cost[e] = 0;
|
|
1090 |
} else {
|
|
1091 |
_forward[u] = false;
|
|
1092 |
_pi[u] = ART_COST;
|
|
1093 |
_source[e] = _root;
|
|
1094 |
_target[e] = u;
|
|
1095 |
_flow[e] = -_supply[u];
|
|
1096 |
_cost[e] = ART_COST;
|
|
1097 |
}
|
| 1117 |
1098 |
}
|
| 1118 |
1099 |
}
|
|
1100 |
else if (_sum_supply > 0) {
|
|
1101 |
// LEQ supply constraints
|
|
1102 |
_search_arc_num = _arc_num + _node_num;
|
|
1103 |
int f = _arc_num + _node_num;
|
|
1104 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
|
1105 |
_parent[u] = _root;
|
|
1106 |
_thread[u] = u + 1;
|
|
1107 |
_rev_thread[u + 1] = u;
|
|
1108 |
_succ_num[u] = 1;
|
|
1109 |
_last_succ[u] = u;
|
|
1110 |
if (_supply[u] >= 0) {
|
|
1111 |
_forward[u] = true;
|
|
1112 |
_pi[u] = 0;
|
|
1113 |
_pred[u] = e;
|
|
1114 |
_source[e] = u;
|
|
1115 |
_target[e] = _root;
|
|
1116 |
_cap[e] = INF;
|
|
1117 |
_flow[e] = _supply[u];
|
|
1118 |
_cost[e] = 0;
|
|
1119 |
_state[e] = STATE_TREE;
|
|
1120 |
} else {
|
|
1121 |
_forward[u] = false;
|
|
1122 |
_pi[u] = ART_COST;
|
|
1123 |
_pred[u] = f;
|
|
1124 |
_source[f] = _root;
|
|
1125 |
_target[f] = u;
|
|
1126 |
_cap[f] = INF;
|
|
1127 |
_flow[f] = -_supply[u];
|
|
1128 |
_cost[f] = ART_COST;
|
|
1129 |
_state[f] = STATE_TREE;
|
|
1130 |
_source[e] = u;
|
|
1131 |
_target[e] = _root;
|
|
1132 |
_cap[e] = INF;
|
|
1133 |
_flow[e] = 0;
|
|
1134 |
_cost[e] = 0;
|
|
1135 |
_state[e] = STATE_LOWER;
|
|
1136 |
++f;
|
|
1137 |
}
|
|
1138 |
}
|
|
1139 |
_all_arc_num = f;
|
|
1140 |
}
|
|
1141 |
else {
|
|
1142 |
// GEQ supply constraints
|
|
1143 |
_search_arc_num = _arc_num + _node_num;
|
|
1144 |
int f = _arc_num + _node_num;
|
|
1145 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
|
1146 |
_parent[u] = _root;
|
|
1147 |
_thread[u] = u + 1;
|
|
1148 |
_rev_thread[u + 1] = u;
|
|
1149 |
_succ_num[u] = 1;
|
|
1150 |
_last_succ[u] = u;
|
|
1151 |
if (_supply[u] <= 0) {
|
|
1152 |
_forward[u] = false;
|
|
1153 |
_pi[u] = 0;
|
|
1154 |
_pred[u] = e;
|
|
1155 |
_source[e] = _root;
|
|
1156 |
_target[e] = u;
|
|
1157 |
_cap[e] = INF;
|
|
1158 |
_flow[e] = -_supply[u];
|
|
1159 |
_cost[e] = 0;
|
|
1160 |
_state[e] = STATE_TREE;
|
|
1161 |
} else {
|
|
1162 |
_forward[u] = true;
|
|
1163 |
_pi[u] = -ART_COST;
|
|
1164 |
_pred[u] = f;
|
|
1165 |
_source[f] = u;
|
|
1166 |
_target[f] = _root;
|
|
1167 |
_cap[f] = INF;
|
|
1168 |
_flow[f] = _supply[u];
|
|
1169 |
_state[f] = STATE_TREE;
|
|
1170 |
_cost[f] = ART_COST;
|
|
1171 |
_source[e] = _root;
|
|
1172 |
_target[e] = u;
|
|
1173 |
_cap[e] = INF;
|
|
1174 |
_flow[e] = 0;
|
|
1175 |
_cost[e] = 0;
|
|
1176 |
_state[e] = STATE_LOWER;
|
|
1177 |
++f;
|
|
1178 |
}
|
|
1179 |
}
|
|
1180 |
_all_arc_num = f;
|
|
1181 |
}
|
| 1119 |
1182 |
|
| 1120 |
1183 |
return true;
|
| 1121 |
1184 |
}
|
| 1122 |
1185 |
|
| 1123 |
1186 |
// Find the join node
|
| 1124 |
1187 |
void findJoinNode() {
|
| ... |
... |
@@ -1371,26 +1434,14 @@
|
| 1371 |
1434 |
updateTreeStructure();
|
| 1372 |
1435 |
updatePotential();
|
| 1373 |
1436 |
}
|
| 1374 |
1437 |
}
|
| 1375 |
1438 |
|
| 1376 |
1439 |
// Check feasibility
|
| 1377 |
|
if (_sum_supply < 0) {
|
| 1378 |
|
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
| 1379 |
|
if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
|
| 1380 |
|
}
|
| 1381 |
|
}
|
| 1382 |
|
else if (_sum_supply > 0) {
|
| 1383 |
|
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
| 1384 |
|
if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
|
| 1385 |
|
}
|
| 1386 |
|
}
|
| 1387 |
|
else {
|
| 1388 |
|
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
| 1389 |
|
if (_flow[e] != 0) return INFEASIBLE;
|
| 1390 |
|
}
|
|
1440 |
for (int e = _search_arc_num; e != _all_arc_num; ++e) {
|
|
1441 |
if (_flow[e] != 0) return INFEASIBLE;
|
| 1391 |
1442 |
}
|
| 1392 |
1443 |
|
| 1393 |
1444 |
// Transform the solution and the supply map to the original form
|
| 1394 |
1445 |
if (_have_lower) {
|
| 1395 |
1446 |
for (int i = 0; i != _arc_num; ++i) {
|
| 1396 |
1447 |
Value c = _lower[i];
|
| ... |
... |
@@ -1398,12 +1449,36 @@
|
| 1398 |
1449 |
_flow[i] += c;
|
| 1399 |
1450 |
_supply[_source[i]] += c;
|
| 1400 |
1451 |
_supply[_target[i]] -= c;
|
| 1401 |
1452 |
}
|
| 1402 |
1453 |
}
|
| 1403 |
1454 |
}
|
|
1455 |
|
|
1456 |
// Shift potentials to meet the requirements of the GEQ/LEQ type
|
|
1457 |
// optimality conditions
|
|
1458 |
if (_sum_supply == 0) {
|
|
1459 |
if (_stype == GEQ) {
|
|
1460 |
Cost max_pot = std::numeric_limits<Cost>::min();
|
|
1461 |
for (int i = 0; i != _node_num; ++i) {
|
|
1462 |
if (_pi[i] > max_pot) max_pot = _pi[i];
|
|
1463 |
}
|
|
1464 |
if (max_pot > 0) {
|
|
1465 |
for (int i = 0; i != _node_num; ++i)
|
|
1466 |
_pi[i] -= max_pot;
|
|
1467 |
}
|
|
1468 |
} else {
|
|
1469 |
Cost min_pot = std::numeric_limits<Cost>::max();
|
|
1470 |
for (int i = 0; i != _node_num; ++i) {
|
|
1471 |
if (_pi[i] < min_pot) min_pot = _pi[i];
|
|
1472 |
}
|
|
1473 |
if (min_pot < 0) {
|
|
1474 |
for (int i = 0; i != _node_num; ++i)
|
|
1475 |
_pi[i] -= min_pot;
|
|
1476 |
}
|
|
1477 |
}
|
|
1478 |
}
|
| 1404 |
1479 |
|
| 1405 |
1480 |
return OPTIMAL;
|
| 1406 |
1481 |
}
|
| 1407 |
1482 |
|
| 1408 |
1483 |
}; //class NetworkSimplex
|
| 1409 |
1484 |
|