Changeset 663:8b0df68370a4 in lemon-main
- Timestamp:
- 05/12/09 12:06:40 (16 years ago)
- Branch:
- default
- Phase:
- public
- Files:
-
- 1 added
- 3 edited
-
doc/Makefile.am (modified) (1 diff)
-
doc/groups.dox (modified) (3 diffs)
-
doc/min_cost_flow.dox (added)
-
lemon/network_simplex.h (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
doc/Makefile.am
r587 r663 9 9 doc/mainpage.dox \ 10 10 doc/migration.dox \ 11 doc/min_cost_flow.dox \ 11 12 doc/named-param.dox \ 12 13 doc/namespaces.dox \ -
doc/groups.dox
r660 r663 336 336 337 337 /** 338 @defgroup min_cost_flow Minimum Cost Flow Algorithms338 @defgroup min_cost_flow_algs Minimum Cost Flow Algorithms 339 339 @ingroup algs 340 340 … … 342 342 343 343 This group contains the algorithms for finding minimum cost flows and 344 circulations. 345 346 The \e minimum \e cost \e flow \e problem is to find a feasible flow of 347 minimum total cost from a set of supply nodes to a set of demand nodes 348 in a network with capacity constraints (lower and upper bounds) 349 and arc costs. 350 Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$, 351 \f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and 352 upper bounds for the flow values on the arcs, for which 353 \f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, 354 \f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow 355 on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the 356 signed supply values of the nodes. 357 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ 358 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with 359 \f$-sup(u)\f$ demand. 360 A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution 361 of the following optimization problem. 362 363 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 364 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq 365 sup(u) \quad \forall u\in V \f] 366 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 367 368 The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be 369 zero or negative in order to have a feasible solution (since the sum 370 of the expressions on the left-hand side of the inequalities is zero). 371 It means that the total demand must be greater or equal to the total 372 supply and all the supplies have to be carried out from the supply nodes, 373 but there could be demands that are not satisfied. 374 If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand 375 constraints have to be satisfied with equality, i.e. all demands 376 have to be satisfied and all supplies have to be used. 377 378 If you need the opposite inequalities in the supply/demand constraints 379 (i.e. the total demand is less than the total supply and all the demands 380 have to be satisfied while there could be supplies that are not used), 381 then you could easily transform the problem to the above form by reversing 382 the direction of the arcs and taking the negative of the supply values 383 (e.g. using \ref ReverseDigraph and \ref NegMap adaptors). 384 However \ref NetworkSimplex algorithm also supports this form directly 385 for the sake of convenience. 386 387 A feasible solution for this problem can be found using \ref Circulation. 388 389 Note that the above formulation is actually more general than the usual 390 definition of the minimum cost flow problem, in which strict equalities 391 are required in the supply/demand contraints, i.e. 392 393 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) = 394 sup(u) \quad \forall u\in V. \f] 395 396 However if the sum of the supply values is zero, then these two problems 397 are equivalent. So if you need the equality form, you have to ensure this 398 additional contraint for the algorithms. 399 400 The dual solution of the minimum cost flow problem is represented by node 401 potentials \f$\pi: V\rightarrow\mathbf{Z}\f$. 402 An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem 403 is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$ 404 node potentials the following \e complementary \e slackness optimality 405 conditions hold. 406 407 - For all \f$uv\in A\f$ arcs: 408 - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; 409 - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; 410 - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. 411 - For all \f$u\in V\f$ nodes: 412 - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$, 413 then \f$\pi(u)=0\f$. 414 415 Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc 416 \f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e. 417 \f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f] 418 419 All algorithms provide dual solution (node potentials) as well, 420 if an optimal flow is found. 421 422 LEMON contains several algorithms for solving minimum cost flow problems. 344 circulations. For more information about this problem and its dual 345 solution see \ref min_cost_flow "Minimum Cost Flow Problem". 346 347 LEMON contains several algorithms for this problem. 423 348 - \ref NetworkSimplex Primal Network Simplex algorithm with various 424 349 pivot strategies. … … 429 354 - \ref CancelAndTighten The Cancel and Tighten algorithm. 430 355 - \ref CycleCanceling Cycle-Canceling algorithms. 431 432 Most of these implementations support the general inequality form of the433 minimum cost flow problem, but CancelAndTighten and CycleCanceling434 only support the equality form due to the primal method they use.435 356 436 357 In general NetworkSimplex is the most efficient implementation, -
lemon/network_simplex.h
r643 r663 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 … … 34 34 namespace lemon { 35 35 36 /// \addtogroup min_cost_flow 36 /// \addtogroup min_cost_flow_algs 37 37 /// @{ 38 38 … … 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 … … 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 … … 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 … … 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 … … 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) { … … 329 298 const CostVector &_pi; 330 299 int &_in_arc; 331 int _ arc_num;300 int _search_arc_num; 332 301 333 302 public: … … 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 … … 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) { … … 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 … … 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 } … … 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) { … … 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 … … 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 … … 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), … … 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) { … … 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 … … 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 … … 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), … … 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]]); … … 679 650 _arc_num = countArcs(_graph); 680 651 int all_node_num = _node_num + 1; 681 int all_arc_num = _arc_num +_node_num;682 683 _source.resize( all_arc_num);684 _target.resize( all_arc_num);685 686 _lower.resize( all_arc_num);687 _upper.resize( all_arc_num);688 _cap.resize( all_arc_num);689 _cost.resize( all_arc_num);652 int max_arc_num = _arc_num + 2 * _node_num; 653 654 _source.resize(max_arc_num); 655 _target.resize(max_arc_num); 656 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 … … 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) … … 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(); … … 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]; 1117 } 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 } 1098 } 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; 1118 1181 } 1119 1182 … … 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 … … 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;
Note: See TracChangeset
for help on using the changeset viewer.
